# A Mortality Model for Multi-populations:
A Semi-Parametric Approach††thanks: The authors gratefully acknowledge
financial support from the Deutsche Forschungsgemeinschaft through the
International Research Training Group IRTG 1792 "High Dimensional Non
Stationary Time Series" and the Collaborative Research Center CRC 649
"Economic Risk".

Lei Fang1, Wolfgang K. Härdle1, Juhyun Park2
1Humboldt-Universität zu Berlin, 2Lancaster University
Corresponding author: Department of Mathematics and Statistics, Lancaster
University, Lancaster LA1 4YF, UK., email<EMAIL_ADDRESS>
###### Abstract

Mortality is different across countries, states and regions. Several empirical
research works however reveal that mortality trends exhibit a common pattern
and show similar structures across populations. The key element in analyzing
mortality rate is a time-varying indicator curve. Our main interest lies in
validating the existence of the common trends among these curves, the similar
gender differences and their variability in location among the curves at the
national level. Motivated by the empirical findings, we make the study of
estimating and forecasting mortality rates based on a semi-parametric
approach, which is applied to multiple curves with the shape-related nonlinear
variation. This approach allows us to capture the common features contained in
the curve functions and meanwhile provides the possibility to characterize the
nonlinear variation via a few deviation parameters. These parameters carry an
instructive summary of the time-varying curve functions and can be further
used to make a suggestive forecast analysis for countries with barren data
sets. In this research the model is illustrated with mortality rates of Japan
and China, and extended to incorporate more countries.

JEL classification: C14, C32, C38, J11, J13

Keywords: Mortality forecasting; Common trend; Lee-Carter method; Multi-
populations; Semi-parametric modeling

##  1 Introduction

In recent years, global population trend has received widespread attention
with regard to the growth of the aging populations on the globe, which raises
demographic risk in most developed and even some developing countries.
Demographic risk is understood to be an imbalance of the age distribution of a
society with the obvious implications in economic growth, social stability,
political decisions and resource allocation. The factor demography is
particularly important in understanding the challenges of developing countries
and their global impacts, yet, due to limited access and available resources,
most studies have focused on the cases from developed countries. This
motivated us to develop a novel approach to forecasting a population trend
with limited data, using the example from China.

China, as a large developing Asian country, is experiencing the transformation
to an aging society at an even faster rate, and is therefore a good example
with which to study demographic risk. However, due to political reasons and
delays in construction of a national system for systematically collecting
statistics, the problem of insufficient and unsatisfying demographic data sets
remains unsolved, and this situation poses a unique challenge in statistical
analysis. Japan, China’s neighbour, has undergone a dramatic demographic
change during the last several decades. In addition, the Japanese government
has set up a complete national statistical system in the middle of the last
century, which provides fine demographic data sets in longer time horizons to
help researchers explore Japan’s demographic transition. A typical example of
demographic data in terms of mortality and fertility is shown in Figure 1 from
Japan.

Figure 1: Japan’s descriptive demography - male (left) and female (middle)
mortality dynamics from 1947 to 2012, fertility (right) dynamics from 1947 to
2009, rates in different years are plotted in rainbow palette order.

MuPoMo

As we shall demonstrate later with our new formulation of the problem, those
Japanese data sets could offer a good reference point to studying the pattern
in China. However, due to the time-lag in the phase of the development of each
country, a direct comparison is not feasible. In this work, we propose to
incorporate a nonlinear transformation across different periods by formulating
a semi-parametric regression model. The advantages and shortcomings of our
approach will be discussed in detail. To overcome some of the shortcomings,
our approach will be extended by taking into account a trend in multi-
populations using global Human Mortality Database (2013).

Our approach could be understood as an extension of the _coherent forecasting_
advocated by (Li & Lee, 2005) with the idea that a local trend estimation
could be improved by taking into account the patterns in a larger group it
belongs to. This requires a serious consideration of the notion of similarity
between the members. Until now, the idea of coherent forecasting was mostly
discussed in the context of concurrent comparison with a small well-defined
neighboring countries. In this article we attempt to demonstrate its
generality in a wider context, with a new formulation under the shape
invariant models (Lawton et al., 1972). We first introduce the datasets at our
disposal to demonstrate our approach before presenting a brief literature
review.

###  1.1 Data Sets

The demographic data sets are collected from different sources: China data
sets are extracted from the China Statistical Year Book, the other 35
countries are obtained from the Human Mortality Database (HMD). For data sets
on China the mortality rates are gender-age specific starting from 0 to 90+
years old, while the mortality rates of other countries are gender-age
specific from 0 to 110+ since the Human Mortality Database makes estimates and
adjustments on the raw data to extend into wider age group. As for the sample
size, they are all different as well: China mortality data spans from 1994 to
2010 with missing data of years 1996, 1997, 2001 and 2006, while the sample
size from other countries range from 14 years (Chile) to 261 years (Sweden).
Referring to the missing values, we use moving average of neighboring five
years to compute these.

Recall that the definition of the mortality rate is the number of deaths per
1000 living individuals per single calender year. To fit the mortality trend
more precisely and for visual convenience, we present the log mortality.

Data and sources codes are available online
(https://github.com/QuantLet/MuPoMo) and linked to the figures.

###  1.2 Literature Review

Since 1980, one of the challenges in demography has been to analyze and
forecast mortality in a purely statistical way without involving the
subjective opinions of experts. Lee & Carter (1992) (LC) firstly proposed a
stochastic method based on a Singular Value Decomposition technique to explore
the unobserved demographic information. This proved insightful and gained a
good reputation. Since then, several methods based on stochastic population
modelling and forecasting have been developed, see e.g. Cairn et al. (2008),
Booth & Tickle (2008) and Booth (2006) for review. Among all the stochastic
models, the most popular one is the LC model, which was used to analyze the
U.S. mortality rates from 1933 to 1987. Based on the idea of the LC model,
comparisons of different methods and some variants or extensions have been
developed. Lee & Miller (2001) compared the forecasts of LC model with the
U.S. social security system forecasts. Li et al. (2004) proposed another
method when there are fewer observations at uneven intervals, and applied it
to China and South Korea. Hyndman & Ullah (2007) developed a more general
method by treating the underlying demographic process as functional data,
employing the functional principal components analysis to extract more than
one explanatory components and providing robust estimation and forecast. This
idea is further extended to multilevel functional data in Shang (2016); Gao &
Shang (2017). These linear extensions are successful in capturing complex
variation at the aggregate level at the expanse of ease of interpretability of
the LC model.

In light of limited data access combined with fragile quality, less
technically refined methods are available for Asian countries compared to, for
example, developed western countries. An exception is the stochastic
population approach on Asian data by Li et al. (2004), who implemented the LC
model to sparse data. In their work, they generated central forecast with just
the first and last observations along the time horizon, improved the estimates
by additional observations and evaluated its performance with other existing
methods. Raftery et al. (2012) proposed a Bayesian method for probabilistic
population projections for all countries, where the Bayesian hierarchical
models, estimated via Markov chain Monte Carlo, are applied to the United
Nations population data. In cases of limited data and similar demographic
trends between two populations (regional or national level), the Bayesian
stochastic modelling for two populations is proposed by Cairns et al. (2011).
This motivates us to analyze Chinese demography via taking Japan into
reference.

One might wonder why Japan, not Taiwan or any other neighboring countries. An
interesting finding from Fang & Härdle (2015) is that China has a demographic
trend closer to Japan than to Taiwan, particularly visible in the mortality
trend. On the economic level, China and Japan have been both important
economies in last several decades and the development pattern is also quite
similar. Hanewald (2011) found that the LC mortality index $k_{t}$ correlates
significantly with macroeconomic fluctuations in some periods, which provides
a good reference with which to connect the mortality trends between China and
Japan. Nevertheless, due to their differences in the developmental phase, it
is much desirable to incorporate nonlinear transformation of the trend to
capture the effect of time trend in a flexible manner.

###  1.3 Goal and Outline

Due to the sensitivity of fertility to social policies and induced
unpredictability, our research is restricted to mortality analysis. The
purpose of our analysis is to build a modeling framework that takes into
account similarities in trend across different countries. We focus on the
popular time-varying mortality indicator $k_{t}$ extracted from the standard
LC method. Following the usual interpretation of the mortality indicator as a
time trend, we seek to compare the trends in terms of their shape variation in
time or in phase to measure similarities among different countries. Most
linear approaches such as those based on principal component analysis are
inefficient in capturing these nonlinear variations and thus difficult to
interpret the results. Instead we directly model the shape variation of the
curves. To increase flexibility and interpretability of the shape of the trend
function, we adopt the framework of the semiparametric comparison approach
(Härdle & Marron, 1990; Kneip & Engel, 1995) and accordingly demonstrate
potential improvements in mortality forecasting. We analyse similarities in
mortality between China and Japan, and then extend our approach to a global
common mortality trend and study a sub-group pattern.

A typical setting for the semiparametric approach assumes that (i)
measurements are available on a common interval, often at common and dense
grid points, (ii) the population has a well defined trend in terms of
identifiable features and (iii) the errors are independent, all of which are
violated in our demographic application. We demonstrate its ramifications and
additional considerations throughout the analysis.

The general methodology is presented in Section 2. Section 3 focuses on
mortality forecast with two country comparison, using the example of China and
Japan. This serves a motivation for further development in Section 4 with
multiple populations using the Human Mortality Database. More discussions on
economic insights, global aging trend influences and suggestions are given in
the last section.

##  2 Methodology

In this section, we firstly introduce the parameters of interests and then
outline the LC method, semi-parametric comparison of nonlinear curves and
common trend modeling in details.

###  2.1 Notations and Parameters of Interest

The parameters of interests are age-specific mortality rates from multiple
countries. We use the symbols $m$ to denote the mortality rate. All the
parameters are indexed by a one-year age group, denoted by $x$, and in
addition indexed by time, denoted by $t$. For instance, $m(x,t)$ is the
mortality rate for age $x$ in year $t$. Based on the LC model, $m(x,t)$ is
assumed to be decomposed into the average age pattern and the time-varying
index $k_{t}$ (for single country, state or region). Characterizing $k_{t}$ is
essential in understanding the mortality trend, which is the main goal of this
paper. When it comes to multiple countries, we use $k_{i}(t)$ to denote the
derived time-varying mortality indicator for country $i$, with
$i\in\\{1,...,N\\}$.

###  2.2 Lee-Carter (LC) Method

The benchmark LC model employs the Singular Value Decomposition (SVD) to
analyse the time series on the log of the age-specific mortality. The method
relies on the standard statistical analysis of the time series. Nonetheless,
the LC model does not fit well in some cases where missing data is common or
the horizon of time series is not sufficient, the reason being the assumption
of long-term stationarity.

The basic idea for demography dynamics analysis is to regress mortality
$m(x,t)$ on non-observable regressors for prediction. The regressors are
obtained via SVD of the demographic indicators. It separates the age pattern
from the time-dependent components, takes time series analysis on the time-
dependent components only and hence forecasts the future trend.

The mortality rate $m(x,t)$ is hence calibrated via the following model:

$\displaystyle\log\\{m(x,t)\\}$ $\displaystyle=$ $\displaystyle
a_{x}+b_{x}k_{t}+\varepsilon_{x,t},$ (1)

or

$\displaystyle m(x,t)$ $\displaystyle=$
$\displaystyle\mbox{exp}(a_{x}+b_{x}k_{t}+\varepsilon_{x,t}),$

where $a_{x}$ is the derived age pattern averaged across years and $k_{t}$
represents the only time-varying index of mortality level. Thus $b_{x}$
measures the sensitivity of the mortality rates to the change of $k_{t}$,
reflecting how fast the mortality rate changes over ages and
$\varepsilon_{x,t}$ is the residual term at age $x$ in year $t$ with
$\mathop{\mbox{\sf E}}(\varepsilon_{x,t})=0$ and $\mathop{\mbox{\sf
Var}}(\varepsilon_{x,t})=\sigma_{\varepsilon}^{2}$.

Three unobserved parameters $a_{x}$, $b_{x}$ and $k_{t}$ in the single
equation (1) mean that the LC model is over-parameterized and therefore two
normalisation constraints are imposed:

$\displaystyle\sum k_{t}=0,\ \sum b_{x}=1\,.$

By SVD, one obtains $k_{t}$ and $b_{x}$. The evolution of $k_{t}$ can be
further fitted by standard time series models such as using ARIMA techniques.
Lee & Carter (1992) found that a random walk with drift describes $k_{t}$
quite well:

$\displaystyle k_{t}=k_{t-1}+d+e_{t}$

where $d$ is the drift parameter reflecting the average annual change and
$e_{t}$ is an uncorrelated error. Others such as Chan et al. (2008) pointed
out that the mortality index $k_{t}$ may be better fitted with a trend-
stationary model for Canada, England and United States when accounting for a
break in mortality rate decline during the 1970s.

Given an $h$-step ahead forecasting $k_{t+h}$ from the time series models, the
forecast of the mortality rates in future period $t+h$ can be made via the
following formula:

$\displaystyle m(x,t+h)$ $\displaystyle=$
$\displaystyle\mbox{exp}(a_{x}+b_{x}k_{t+h}).$

###  2.3 Semi-Parametric Comparison of Nonlinear Curves

When the observable curves are noisy versions of similar regression curves,
comparison of regression curves from related samples is not trivial. The
problem aggravates when the data are sparse, partially available for some and
the domain of the curves do not match, as demonstrated in Figure 2 where we
wish to compare the short segment of black dots on the right and the long blue
curve below. The black dots on the right with the red solid curve underneath
represent the mortality rates from China, while the elongated blue curve below
represents the mortality rates from Sweden. The curve above in cyan is a
transformed Swedish curve that gives a best match to the partial observations
of China, which will be defined precisely later.

Figure 2: Semi-parametric Comparison of Nonlinear Curves: the dark blue curve
and the red curve have similar pattern, while the light blue curve is semi-
parametrically shifted to represent the red curve via the dark blue curve.

MuPoMo

It is not unreasonable to assume that developing countries will catch up with
the developed countries, and due to globalization the trend in mortality rates
is expected to follow a similar pattern. Then the main differences among those
curves could be explained by shifted time axis and vertical re-scaling, which
could further be parametrized for parsimonious representation and ease of
interpretation. This is our motivation to take a view of semiparametric
comparison in taking into account mortality variations of other countries.

To handle the noisy data, nonparametric smoothing techniques could be
incorporated to estimate the underlying curves when the solid theory is
unavailable in modeling them. Hence, we seek general semi-parametric models
that allow for nonlinear transformation on the domain as well as the image of
the curves.

Simply denoting the underlying curves by $f_{1}$ and $f_{2}$, the semi-
parametric comparison of the nonlinear curves can be expressed as

$\displaystyle
f_{2}(t)=\theta_{1}f_{1}\left(\frac{t-\theta_{2}}{\theta_{3}}\right)+\theta_{4},$
(2)

where we assume that $f_{2}$ has a similar pattern to $f_{1}$ and
$\theta=(\theta_{1},\theta_{2},\theta_{3},\theta_{4})^{\top}$ are shape
deviation parameters. Since normally only noisy measurements of the curves are
available, $f_{1}$ and $f_{2}$ are not directly available, there will be
measurements models given by

$Y_{i}(t_{ij})\equiv Y_{ij}=f_{i}(t_{ij})+\varepsilon_{ij}\qquad
i=1,2;j=1,\ldots,n_{i}\,.$

More detailed discussions on this method can be referred to Härdle & Marron
(1990) on semiparametric comparison of regression curves.

###  2.4 Common Trend Modeling

When more than two regression curves share similar pattern or trend, a common
trend model can be built based on the technique of semi-parametric comparison
of nonlinear curves we discussed previously. Suppose we are given $N$ noisy
curves $\mbox{Y}_{i},i=1,\ldots,N$ that exhibit some similar patterns. A
general regression model can be expressed as,

$\mbox{Y}_{i}=f_{i}+\varepsilon_{i}\,,$ (3)

where $f_{i}$ denote unknown smoothing regression functions while
$\varepsilon_{i}$ represent independent errors with mean 0 and variance
$\sigma_{i}^{2}$.

The relationship among these similar curves can be described as

$\displaystyle
f_{i}(t)=\theta_{i1}g\left(\frac{t-\theta_{i2}}{\theta_{i3}}\right)+\theta_{i4}.$
(4)

Here $\theta_{i}=\left(\theta_{i1},\theta_{i2},\theta_{i3},\theta_{i4}\right)$
are unknown parameters describing shape deviations, and $g$ is a unknown
function specifying the common shape of these curves, which can be interpreted
as a reference curve. The model in (4) is commonly known as shape invariant
model (SIM), firstly proposed by Lawton et al. (1972) and further studied by
Kneip & Engel (1995) and provides an extension of the model in (2) to multiple
curves. A detailed investigation of this model to estimate the mortality trend
will be given in the following section.

##  3 Mortality forecasting: two-country case

(a)  (b)

Figure 3: China mortality trend vs. Japan mortality trend: female, male
(left) and Japan trend, Japan smoothed trend, China trend  and China smoothed
trends (right).The dotted lines (right) from left to right are shifted Japan’s
smoothed trends of 20-, 23- and 25- years forward respectively.

MuPoMo

With reasons mentioned in Section 1, we will first focus on analysing
mortality similarities between China and Japan. Before presenting our
formulation, we first graphically show the empirical trends in the left of
Figure 3, which suggests that the mortality trends from both gender groups of
China correlate with those of Japan respectively. However, due to limited
sample they all seemingly reflect linear mortality trend over time, which
would contradict to the theoretical and conceptual views on mortality trends.
To provide intuitive comparison between these two countries, the horizontal
shift of Japan’s mortality curve over time axis is plotted in the right part
of Figure 3. The dotted lines from left to right are shifted Japan’s smoothed
trends of 20-, 23- and 25- years forward respectively. Graphically we see that
Japan’s mortality trend is 23 years earlier than China.

###  3.1 Model formulation

To parameterize the potential relationship between China and Japan mortality
trend, we specify the model as following, using $k_{t}$ derived from LC model
in (1):

$\displaystyle\log\\{m(x,t)\\}$ $\displaystyle=$ $\displaystyle
a_{x}+b_{x}k_{t}+\varepsilon_{x,t}\,.$

Then we infer China’s mortality trend via Japan’s trend through the technique
of semi-parametric comparison of regression curves defined in (2)

$\displaystyle
k_{c}(t)=\theta_{1}k_{j}\left(\frac{t-\theta_{2}}{\theta_{3}}\right)+\theta_{4},$
(5)

where $k_{c}(t)$ is the time-varying indicator for China, $k_{j}(t)$ is the
time-varying indicator for Japan, and
$\theta=(\theta_{1},\theta_{2},\theta_{3},\theta_{4})^{\top}$ are shape
deviation parameters.

##### Understanding $\theta$

It is probably easiest to interpret the parameters (5) by starting with
$\theta=(\theta_{1},\theta_{2},\theta_{3},\theta_{4})^{\top}=(1,\theta_{2},1,\theta_{4})^{\top}$.

  * •

$\theta_{1}$ is the general trend adjustment, here selected as 1.

  * •

$\theta_{2}$ is the time-delay parameter

  * •

$\theta_{3}$ is the time acceleration parameter, here selected as 1.

  * •

$\theta_{4}$ is the vertical shift parameter

In Figure 4, it demonstrates how $\theta$ influences shift of these two
curves. In the left of Figure 4, China’s female $k_{t}$ can reach a similar
behavioral area by shifting horizontally $\theta_{2}=-23$, while in the right
it shows that another acceptable area could be obtained via approximately
vertical shift of $\theta_{4}=85$.

(a)  (b)

Figure 4: Time delay $\theta_{2}=-23$ (left) and vertical shift
$\theta_{4}=85$ (right).

MuPoMo

###  3.2 Estimation

In order to find the optimal solution for shape deviation parameters, we
minimize the following loss function.

$\displaystyle\min_{\theta}\int_{t_{c}}\left\\{\hat{k}_{c}(u)-\theta_{1}\hat{k}_{j}\left(\frac{u-\theta_{2}}{\theta_{3}}\right)-\theta_{4}\right\\}^{2}w(u)du,$
(6)

where $\hat{k}_{c}(t)$ and $\hat{k}_{j}(t)$ are the nonparametric estimates of
the original time-varying indicators $k_{c}(t)$ and $k_{j}(t)$, and $t_{c}$ is
the time interval of China mortality data. Here we have used local linear
smoother for our implementation but any other smoothing methods could be used
(Simonoff, 1996).

The comparison region needs to satisfy the following condition, in order to
make sure the parameter estimation is compared only in the common region
defined by $w(u)$.

$\displaystyle w(u)=\prod_{t_{j}}1_{[a,b]}\\{(u-\theta_{2})/\theta_{3}\\},$

where $t_{j}$ is the time interval of Japan’s mortality data, $a\geq
inf(t_{j})$ and $b\leq sup(t_{j})$. To consider the importance of more recent
data’s impact on future trend, one could impose different weights on different
support intervals.

##### Algorithm

To estimate the parameters by the nonlinear least squares estimation criterion
given in (6), we first obtain the estimates of $k_{c}$ and $k_{t}$ by
nonparametric local linear smoothing, denoted by $\hat{k}_{c}$ and
$\hat{k}_{t}$ respectively. Then we set up the initial estimates
$\theta^{0}=(\theta_{1}^{0},\theta_{2}^{0},\theta_{3}^{0},\theta_{4}^{0})$ and
solve the nonlinear least squares estimation problem by iteratively updating
the estimates until convergence.

(a)

(b)

Figure 5: Loss surface of $\theta_{2}$ and $\theta_{4}$ (left) and Contour of
$\theta_{2}$ and $\theta_{4}$ (right) with $\theta_{1}=\theta_{3}=1$.

MuPoMo

##### Initial choice of $\theta_{2}$ and $\theta_{4}$

Figure 6: Loss function of $\theta_{2}$ with
$(\theta_{1},\theta_{3},\theta_{4})^{\top}=(1,1,0)^{\top}$.

MuPoMo

From previous Figure 3 and in our initial analysis, we see that there is a
potential ambiguity between $\theta_{2}$ and $\theta_{4}$. It is also clear
that we can find the replacement relationship from Figure 4. To investigate
further, we show the criterion function in Figure 5 as a function of
$(\theta_{2},\theta_{4})$. There is a valley area in the loss surface function
of $\theta_{2}$ and $\theta_{4}$ in the left plot and also in the contour of
$\theta_{2}$ and $\theta_{4}$ in the right one, which suggests that there
exists an approximate linear combination of $\theta_{2}$ and $\theta_{4}$ in
searching $\theta$ for an optimal solution. It will bring about difficulties
in finding the optimal parameters $\theta$ in numerical optimization. This
ambiguity in identifiability of the parameters is one of the ramifications of
the partial observation in our setting.

In order to find the optimal value, we should be very careful with selecting
initial values of $\theta$ and in the consideration of identifiability issues
we need to decide whether the analysis concentrates on time delay or vertical
shift. In our analysis we stick with time delay influence $\theta_{2}$ since
it is more valuable in prediction perspective, and thus the initial value of
$\theta$ is determined as $(\theta_{1},\theta_{3})^{\top}=(1,1)^{\top}$ and
set $\theta_{4}=0$. Hence the optimal initial $\theta_{2}$ is obtained around
-23, see Figure 6.

##### Goodness of Fit

Based on the initial
$\theta^{(0)}=(\theta_{1},\theta_{2},\theta_{3})^{\top}=(1,-23,1)^{\top}$ and
algorithm, the optimal parameter $\theta$ is reached at
$\hat{\theta}=(1.205,-22.621,1.000)^{\top}$. This is used in Figure 7, showing
that after the optimal transformation of the curve of Japan (based on the
optimal value of $\theta$) the $k_{t}$ of China fits quite well in the $k_{t}$
of Japan of years around from 1970 to 1990.

Figure 7: Goodness of Fit and forecast of China’s mortality trend from 2011 to
2030 via Japan’s historical data: black dots represent the original $k_{t}$
from Japan and China, and Japan smoothed trend, China smoothed trend  are
displayed as well; the fitted trend is plotted as light blue dashed line,
while the overlapping part is colored in purple; the light blue dashed line
after year 2011 is the forecast part.

MuPoMo

###  3.3 Forecast

Afterwards we can forecast $k_{t}$ for China via the data from Japan and the
optimal estimated shape deviation parameter $\hat{\theta}$ to extend the
forecasting horizon.

$\displaystyle
k_{c}(t+h)=\hat{\theta}_{1}k_{j}\left\\{\frac{(t+h)-\hat{\theta}_{2}}{\hat{\theta}_{3}}\right\\},$
(7)

where $\hat{\theta}=(1.205,-22.621,1.000)^{\top}$ and
$t=1994,1995,...,2010;h=1,2,...,20$.

Compared with the traditional forecasting method with time series analysis,
our proposed method can extend the forecasting horizon from 5 years to 23
years, which is a big advantage from this semi-parametric comparison technique
of regression curves, see Figure 7. However, searching for numerical optimal
solutions for seemingly linear regression curves is still a challenge in this
context. This has motivated us to develop an extension to the multi-
populations mortality models, because the identifiability issues can be better
solved by borrowing information across multiple curves.

##  4 Mortality forecasting: Multi-population models

Based on the previous study of two countries under sparse data, investigation
of multi-populations becomes promising and necessary since more information on
nonlinear trend will be provided in case of multi-countries. At the same time
it enables us to study global mortality trend in the past century and its
future.

(a)  (b)

Figure 8: Similarities of mortality trend among countries: different colors
represent different countries (left), and the red thick curve in the right
plot stands for reference curve while grey ones are smoothed curves.

MuPoMo

Denote by $k_{i}$ the time-varying mortality indicator for country $i$, with
$i\in\\{1,...,N\\}$. Figure 8 displays the estimates of $k_{i}$ for 36
countries in the database. The left panel shows the noisy curves that are
originally estimated $k_{t}$ from the LC model without any further
nonparametric smoothing, while on the right it shows the smoothed $k_{t}$ with
an initial estimate of the reference curve overlayed. Later on, we will
discuss why 31 countries are selected for analysis. By design, the available
time measurements vary among countries. Nevertheless, we notice remarkable
similarities in the trend across the countries, subject to individual
variability.

In order to investigate this structural similarity and to borrow the
information across the countries, we consider the shape invariant model
introduced in Section 2.4 of methodology part.

###  4.1 Model formulation

Specifically we assume the additive noise model defined in (3) for the derived
$k_{i}$s from LC model to account for noise and suppose that the underlying
curves share some common trend and can be represented in the form

$\displaystyle
k_{i}(t)=\theta_{i1}k_{0}\left(\frac{t-\theta_{i2}}{\theta_{i3}}\right)+\theta_{i4},$
(8)

where $k_{0}(t)$ is a reference curve, understood as common trend and
$\theta_{i}=(\theta_{i1},\theta_{i2},\theta_{i3},\theta_{i4})^{\top}$ are
shape deviation parameters. In order to be able to interpret the reference
curve $k_{0}$ as a mean trend, we can use the normalizing constraints on the
parameter $\theta_{i}$ as

$\displaystyle
N^{-1}\sum_{i=1}^{N}\theta_{i1}=N^{-1}\sum_{i=1}^{N}\theta_{i3}=1,$ (9)
$\displaystyle
N^{-1}\sum_{i=1}^{N}\theta_{i2}=N^{-1}\sum_{i=1}^{N}\theta_{i4}=0\,.$ (10)

Alternatively, we can use any country as a reference curve, for example,
Sweden as the longest record holder could be a reasonable choice, in which
case the reference curve is set to be $k_{t}$ of Sweden with
$\theta_{0}=(1,0,1,0)$ and $\theta_{i}$ will measure the deviation with
respect to the reference curve. In this work we will consider the mean curve
as a reference curve and use the above normalization constraints.

###  4.2 Estimation

##### Estimation of parameters

Suppose that $k_{0}$ and $k_{i}$ are given. Then for each country $i$, the
parameter $\theta_{i}$ can be determined by minimizing the least squares
criterion as

$\int\Big{\\{}k_{i}(t)-\theta_{i1}k_{0}\Big{(}\frac{t-\theta_{i2}}{\theta_{i3}}\Big{)}-\theta_{i4}\Big{\\}}^{2}w_{i}(t)\,dt$
(11)

where $w_{i}$ is chosen to ensure that the two functions are evaluated over
the common domain as in the case of two countries. In practice, $k_{0}$ and
$k_{i}$ are replaced by its nonparametric estimate.

##### Estimation of Common Trend

For given parameters $\theta_{i},i=1,\ldots,n$, the functional relationship in
(8) implies that

$\displaystyle
k_{i}(\theta_{i3}t+\theta_{i2})=\theta_{i1}k_{0}(t)+\theta_{i4},$ (12)

Thanks to the normalizing conditions on $\theta_{i1}$ and $\theta_{i4}$, this
implies that

$k_{0}(t)=N^{-1}\sum_{i=1}^{N}k_{i}(\theta_{i3}t+\theta_{i2})\,.$ (13)

That is, if $k_{i}$ is appropriately transformed with respect to the
individual parameters $\theta_{i}$, then $k_{0}$ is simply the average. In
practice, we have noisy version of $k_{i}$ available at different number of
time points. Then the functional mean can be estimated more efficiently with
nonparametric smoothing, which essentially gives rises to a weighted average
estimate.

##### Estimation algorithm

Combining the above two steps leads to the following iterative algorithm for
estimation of the parameters.

  1. (a)

Given $k_{i}$, obtain an initial estimate of $k_{0}$ based on all country-
level mortality rates

  2. (b)

Given $k_{0}$, update $\theta_{i}$ by minimizing the nonlinear least squares
criterion in (11) for each $i=1,\ldots,N$.

  3. (c)

Normalize the parameters to satisfy the constraints.

  4. (d)

Given $\theta_{i},i=1,\ldots,N$, update $k_{0}$ by (13).

  5. (e)

Iterate (b)-(d) until convergence.

The general algorithm was proposed and studied by Kneip & Engel (1995). We
adapted accordingly to account for incomplete observations in our sample.

##### Initial values of $k_{0}$ and $k_{i}$

To initialize $k_{0}$, we choose the trimmed mean of the sample estimates,
based on the middle 50% of the countries in terms of the length of the
recording period. The estimation of $k_{0}$ and $k_{i}$ is done with local
linear kernel smoothing method to account for measurement error. The smoothing
parameter for $k_{i}$ was selected to maintain comparable smoothness across
the samples using common degrees of freedom (Bowman & Azzalini, 1997).

##### Computational Issues with the Parameterization

The shape invariant model implicitly assumes that there are identifiable
features that are common across the sample. It is easy to check for densely
observed curves (with non-monotone functions) by means of the derivative
estimation, but for sparsely observed curves, there could be an ambiguity in
identifying the parameters. In the case of the mortality curves, due to the
limited measurements available, the ambiguity occurs in distinguishing the
differential effect of vertical shift ($\theta_{4})$ and horizontal shift
$(\theta_{2})$ in time. In this case, we choose to attribute the effect as
horizontal shift and set $\theta_{4}=0$, as this is more amenable to
interpretation and meanwhile promising to extend forecast horizon. Comparison
of these two parameterization cases will be illustrated in Section 4.3.2.

##### Bootstrap for Prediction Interval

The forecasting for Chinese mortality (8) is now updated based on the common
reference curve $k_{0}$. In order to construct a prediction interval of
China’s mortality trend, we also need a reliable estimate of the measure of
variability. Unlike the standard setting studied by Kneip & Engel (1995) for
relatively densely observed data on common intervals, it is difficult to
derive an asymptotic distribution of the estimators for sparse and incomplete
data as ours. Here, we opt for a bootstrap method to approximate the
uncertainties in estimation and prediction.

The standard bootstrap techniques relying on identically independent
distributed (i.i.d.) observations are not appropriate here. Recently re-sample
methods for dependent data have considered several options: bootstrap with
i.i.d. innovations, bootstrap with block segments and model-based bootstrap.
Due to limited sample of China’s mortality time series, bootstrap with block
segments do not create the ideal re-sampled time series. Alternatively we
bootstrap the mortality data based on i.i.d. innovations obtained from fitting
time series model, and afterwards we carry out estimation on the re-sampled
data and generate prediction interval at different levels.

Suppose that for the time series $k_{1},...,k_{n}$ and some fixed $p\in N$,
there exists a parametric estimator of the conditional expectation
$E(k_{t}|k_{t-1},...,k_{t-p})$ denoted by
$\widehat{m}_{n}(k_{t-1},...,k_{t-p})$. This estimator leads to residuals

$\widehat{e_{t}}:=k_{t}-\widehat{m}_{n}(k_{t-1},...,k_{t-p}),t=p+1,...,n.$
(14)

Resampling from these residuals leads to a bootstrap sample of time series

$k_{t}^{*}=\widehat{m}_{n}(k_{t-1}^{*},...,k_{t-p}^{*})+e_{t}^{*},t=p+1,...,n.$
(15)

The idea of parametric fit to the conditional expectation can be executed by
ARIMA models. Kreiss & Lahiri (2012) discuss different situations with
parametric and nonparametric modeling of predictor $k_{t}$ and respective
asymptotic consistence properties.

###  4.3 Global Trend

####  4.3.1 Outlying countries

In light of similarities across countries, seeking global mortality trend
becomes very reasonable and natural. Some research also find out mortality
trend has connection with economic development level, GDP for instance, see
Hanewald (2011). In the remaining section, we are going through the empirical
analysis of common trend in different groups via slightly different models.

Figure 9: Five outlying countries: geographically neighbors in east Europe.
Source: Google Map.

Evidence tells us that the mortality rate is decreasing as time evolves, due
to medical improvement, economic development and social stability. However, we
notice that there are some remarkable outliers like Russia, Estonia, Latvia,
Lithuania and Belarus in the database.

As shown in Figure 9, all these five countries locate themselves in east
Europe and are used to be members of the Soviet Union. They share similar
geographic characteristics and meanwhile experienced parallel economic and
social progress, and within our expectation they reflect comparable mortality
moving path as well, as displayed in Figure 10 and Appendix 1. Note that solid
(blue) curves represent global mortality trend, dashed (yan) curves are
representing estimated individual country-level mortality trends based on
global trend, the short solid (red) curves are smoothed original individual
country-level mortality trends and (black) circles are original individual
country-level mortality trends.

Surprisingly, they exhibit a quite opposite tendency in contrast with other 31
countries. The mortality rates go through a short period of decrease, then
stay stable or slightly increase for several years and afterwards go back to
declining path again. One possible reason on this different phenomenon perhaps
is connected with political event of dissolution of the soviet union.

(a)  (b)

Figure 10: Different mortality movements of Lithuania and Latvia: blue curves
represent global mortality trend, light blue curves are representing estimated
individual country-level mortality trends based on global trend, red curves
are smoothed original individual country-level mortality trends and black dots
are original individual country-level mortality trends.

MuPoMo

To ideally demonstrate global mortality movements in majority of countries, we
remove these five countries for remaining study to reduce influences from
minor outlying ones.

####  4.3.2 Mortality Trend among the Majority

As discussed previously, for sparsely observed curves there could be an
ambiguity in identifying all of the four parameters. Therefore, we choose to
compare the original 4-parameters model with a simplified 3-parameter model of
setting $\theta_{4}=0$.

(a)  (b)

Figure 11: Common mortality trends estimated by 3-parameters model (left) vs.
4-parameters model (right).

MuPoMo

Under two different parameterization, we estimate the parameters and reference
curve respectively. In Figure 11, the reference curves and common trend for
these two cases are plotted: the left one is calculated from the 3-parameters
model and the right one is generated from the 4-parameters one. The red curve
is the initial reference curve, blue ones are one-step ahead updated reference
curves. From these two plots, no clear and obvious difference can be seen. But
from the viewpoint of analytic thinking, we choose to attribute the effect as
horizontal shift and set $\theta_{4}=0$, as this is more amenable to
interpretation and meanwhile promising to extend forecast horizon.

Figure 12 explains the common mortality trend generated from 3-parameters
model compared with individual nation-level mortality trend. In this graph,
the black solid curve is the initial reference curve, cyan, green, blue and
red ones represent the updated ones at different iteration stage while the
grey ones are the non-smoothed mortality trend from each country. It is
obvious that after one step optimization, reference curves are already showing
a quite similar pattern. The common mortality trend is adjusted to an upper
level, mostly because more developing countries (Czech Republic, Hungary and
China, for example) started collecting demographic data at a later time period
in contrast with more developed countries (such as Sweden, Norway and France),
and also developing countries have higher mortality rates generally. The
figures on illustrating individual case are provided in Appendix 2.

Figure 12: Common mortality trend compared with individual nation-level
mortality trends: red curve is the initial reference curve, blue and black
ones are updated reference curves convergent to common trend, while the grey
lines represent individual country.

MuPoMo

###  4.4 China and Global Mortality Trend

Since a common mortality trend is available, it could be applied to help
improve estimation and forecasting of individual case. Especially when sample
size from individual country is relatively smaller than anticipated, semi-
parametric comparison of common mortality trend with each individual nation-
level one will be a promising way with respect to forecasting. At the same
time, it helps to reduce ambiguity in identifying the parameters in case of
seemingly linear co-movement between regression curves, like the case of
comparing China and Japan.

Figure 13 displays newly estimated China mortality trend via semi-parametric
comparison with common trend. The (blue) solid line is the common trend or
updated reference curve and the (cyan) dashed line is the estimated China
mortality trend based on the common trend. In comparison, the raw China
mortality curve estimate is marked by black circles, with the individual
smoothing estimate overlayed in short (red) line. Thanks to parameters
deviation on time axis $\frac{t-\theta_{2}}{\theta_{3}}$, we could extend
forecasting horizon of China approximately 40 years through the information
from common trend, see Figure 13. Referring to estimation and forecasting of
the other 30 countries, they are arranged in Appendix 2.

Figure 13: Common mortality trend and estimated China mortality trend based on
common trend.

MuPoMo

With model-based bootstrap approach, we simulate 500 re-sampled China’s
mortality time series from 1994 to 2010 based on ARIMA model. From each
simulation, we estimate the optimal shape deviation parameters $\theta$ and
accordingly calculate the estimated China mortality with longer time horizon
based on the common trend.

(a)  (b)

(c)

Figure 14: Histograms of $\hat{\theta}_{1}$,$\hat{\theta}_{2}$ and
$\hat{\theta}_{3}$.

MuPoMo

In Figure 14, we display the variation of
$\hat{\theta}_{1}$,$\hat{\theta}_{2}$ and $\hat{\theta}_{3}$ across the
countries. From histogram of $\hat{\theta}_{1}$, $50\%$ of $\hat{\theta}_{1}$
lies between 1.5 and 2.5, which indicates the overall accelerating declining
mortality trend of China compared with global trend. More than $95\%$ of
$\hat{\theta}_{2}$, the parameter describing time delay, falls into the
interval of $(0,10)$, which further confirms that there exists a time delay of
China’s mortality trend around 10 years later than global situation. Majority
of $\hat{\theta}_{3}$ ranges from 0.98 to 1, which reveals a little time
acceleration in China’s mortality trend.

In Figure 15, confidence intervals at different levels are displayed. On the
left part, confidence intervals at $80\%$ and $90\%$ are plotted in grey zone
and blue zone respectively, while yellow zone highlights the central area of
possible forecast path. On the right one, only $90\%$ confidence interval is
presented. Black line stands for global mortality trend, red one is original
China’s mortality trend and light blue curve shows the estimated China’s
mortality trend based common trend and original China’s data.

(a)  (b)

Figure 15: Confidence intervals at different levels: 80% (left) vs. 90%
confidence interval (right).

MuPoMo

##  5 Discussion

The global mortality, as expected, is undergoing a shift toward exhibiting
declining tendency, and with a dramatic decreasing movement in the last
several decades in contrast with hundreds years ago. It also depicts that most
of countries are converging with a similar mortality pattern of decreasing
over time. The improvement possibly results from economic development and
medical improvement, and it is also not difficult to imagine that there will
be gradually declining mortality rate in the near future due to technology
progress.

In addition to the global mortality trend, each country still behaves
differently from others to some extent. From this perspective, it might be
possible for life insurance companies to design insurance products among
different countries to hedge global longevity risk. That is, if wider range of
countries is covered by a particular insurance company, it is possible to
redistribute longevity risk among them.

Another advantage from this research is to establish a better forecasting
regime to foresee mortality change in longer time horizon, particularly for
countries with limited historical mortality data, such as China and Chile.

## References

  * Booth (2006) Booth, H. (2006).  Demographic forecasting: 1980 to 2005 in review.  International Journal of Forecasting, 22, 547–581.
  * Booth & Tickle (2008) Booth, H., & Tickle, L. (2008).  Mortality modelling and forecasting: a review of methods.  Annals of Actuarial Science, 3, 3–43.
  * Bowman & Azzalini (1997) Bowman, A., & Azzalini, A. (1997).  Applied Smoothing Techniques for Data Analysis: the Kernel Approach with S-Plus Illustrations.  Oxford University Press, Oxford.
  * Cairn et al. (2008) Cairn, A., Blake, D., & Dowd, K. (2008).  Modelling and management of mortality risk: a review.  Scandinavian Actuarial Journal, 108, 79–113.
  * Cairns et al. (2011) Cairns, A., Blake, D., Dowd, K., Coughlan, G., & Khalaf-Allah, M. (2011).  Bayesian stochastic mortality modelling for two populations.  ASTIN Bulletin, 41, 29–55.
  * Chan et al. (2008) Chan, W., Li, S., & Cheung, S. (2008).  Testing deterministic versus stochastic trends in the lee-carter mortality indexes and its implications for projecting mortality improvements at advanced ages.  Living to 100, .
  * Fang & Härdle (2015) Fang, L., & Härdle, W. (2015).  Stochastic population analysis: a functional data analysis.  Technical Report SFB 649 Discussion Paper Series 2015\.
  * Gao & Shang (2017) Gao, Y., & Shang, H. L. (2017).  Multivariate functional time series forecasting: An application to age-specific mortality rates.  Risks, 5(2), 21\.
  * Hanewald (2011) Hanewald, K. (2011).  Explaining mortality dynamics: the role of macroeconomic fluctuations and cause of death trends.  North American Actuarial Journal, 15.
  * Härdle & Marron (1990) Härdle, W., & Marron, J. (1990).  Semiparametric comparison of regression curves.  Annals of Statistics, 18, 63–89.
  * Human Mortality Database (2013) Human Mortality Database (2013).  University of California, Berkeley (USA), and Max Planck Institute for Demographic Research (Germany).  Available at www.mortality.org or www.humanmortality.de (data downloaded on 13 June 2013).
  * Hyndman & Ullah (2007) Hyndman, R., & Ullah (2007).  Robust forecasting of mortality and fertility rates: A functional data approach.  Computational Statistics and Data Analysis, 51, 4942–4956.
  * Kneip & Engel (1995) Kneip, A., & Engel, J. (1995).  Model estimation in nonlinear regression under shape invariance.  Annuals of Statistics, 23, 551–570.
  * Kreiss & Lahiri (2012) Kreiss, J., & Lahiri, S. (2012).  Bootstrap methods for time series.  Handbook of Statistics, 30, 3–26.
  * Lawton et al. (1972) Lawton, W. H., Sylvestre, E. A., & Maggio, M. S. (1972).  Self modeling nonlinear regression.  Technometrics, 14, 513–532.
  * Lee & Carter (1992) Lee, R. D., & Carter, L. R. (1992).  Modeling and forecasting u.s. mortality.  Journal of the American Statistical Association, 87, 659–671.
  * Lee & Miller (2001) Lee, R. D., & Miller, T. (2001).  Evaluating the performance of the lee-carter method for forecasting mortality.  Demography, 38, 537–549.
  * Li & Lee (2005) Li, N., & Lee, R. (2005).  Coherent mortality forecasts for a group of populations: An extension of the lee-cartermethod.  Demography, 42, 575–594.
  * Li et al. (2004) Li, N., Lee, R., & Tuljapurkar, S. (2004).  Using the lee-carter method to forecast mortality for populations with limited data.  International Statistical Review, 72, 19–36.
  * Raftery et al. (2012) Raftery, A., Li, N., Ševčíková, H., Gerland, P., & Heilig, G. (2012).  Bayesian probabilistic population projections for all countries.  Proceedings of the National Academy of Sciences, 109, 13915–13921.
  * Shang (2016) Shang, H. L. (2016).  Mortality and life expectancy forecasting for a group of populations in developed countries: A multilevel functional data method.  The Annals of Applied Statistics, 10(3), 1639–1672.
  * Simonoff (1996) Simonoff, J. S. (1996).  Smoothing methods in statistics.  Springer.

##  6 Appendices

###  6.1 Appendix 1

(a)  (b)

(c)

MuPoMo

###  6.2 Appendix 2

(d)  (e)

(f)  (g)

(h)  (i)

MuPoMo

(j)  (k)

(l)  (m)

(n)  (o)

MuPoMo

(p)  (q)

(r)  (s)

(t)  (u)

MuPoMo

(v)  (w)

(x)  (y)

(z)  (aa)

MuPoMo

(ab)  (ac)

(ad)  (ae)

(af)  (ag)

MuPoMo


# Julia as a Unifying End-to-End Workflow Language
on the Frontier Exascale System

William F. Godoy, Pedro Valero-Lara, Caira Anderson, Katrina W. Lee, Ana
Gainaru, Rafael Ferreira da Silva, and Jeffrey S. Vetter  Oak Ridge National
LaboratoryOak RidgeTNUSA
<EMAIL_ADDRESS>
(2023)

###### Abstract.

We evaluate Julia as a single language and ecosystem paradigm powered by LLVM
to develop workflow components for high-performance computing. We run a Gray-
Scott, 2-variable diffusion-reaction application using a memory-bound, 7-point
stencil kernel on Frontier, the US Department of Energy’s first exascale
supercomputer. We evaluate the performance, scaling, and trade-offs of (i) the
computational kernel on AMD’s MI250x GPUs, (ii) weak scaling up to 4,096 MPI
processes/GPUs or 512 nodes, (iii) parallel I/O writes using the ADIOS2
library bindings, and (iv) Jupyter Notebooks for interactive analysis. Results
suggest that although Julia generates a reasonable LLVM-IR, a nearly 50%
performance difference exists vs. native AMD HIP stencil codes when running on
the GPUs. As expected, we observed near-zero overhead when using MPI and
parallel I/O bindings for system-wide installed implementations. Consequently,
Julia emerges as a compelling high-performance and high-productivity workflow
composition language, as measured on the fastest supercomputer in the world.

Julia, end-to-end workflows, High-Performance Computing, HPC, data analysis,
exascale, Frontier supercomputer, Jupyter notebooks

This manuscript has been authored by UT-Battelle LLC under contract DE-
AC05-00OR22725 with the US Department of Energy (DOE). The publisher
acknowledges the US government license to provide public access under the DOE
Public Access Plan (https://energy.gov/downloads/doe-public-access-plan).

††journalyear: 2023††copyright: licensedusgovmixed††conference: Workshops of
The International Conference on High Performance Computing, Network, Storage,
and Analysis; November 12–17, 2023; Denver, CO, USA††booktitle: Workshops of
The International Conference on High Performance Computing, Network, Storage,
and Analysis (SC-W 2023), November 12–17, 2023, Denver, CO, USA††price:
15.00††doi: 10.1145/3624062.3624278††isbn: 979-8-4007-0785-8/23/11

##  1\. Introduction

Recent emphasis on the end-to-end workflow development process for high-
performance computing (HPC) applications acknowledges the increasing
complexity required for achieving performance, portability, and productivity
(Ben-Nun et al., 2020). This complexity is primarily driven by two factors:
(i) the evolving application requirements for experimental, observational, and
computational science and (ii) the extreme heterogeneity of our computing and
data generation and processing systems (Ferreira da Silva et al., 2017;
Deelman et al., 2018; Vetter et al., 2018). Consequently, the community vision
for a research and development roadmap (Ferreira da Silva et al., 2021, 2023)
has identified key challenges posed by integrating HPC, AI, and FAIR
(findable, accessible, interoperable, and reusable) (Sirvent et al., 2022;
Goble et al., 2020) workflows at exascale.

The Julia programming language (Bezanson et al., 2017) was designed to provide
a powerful unifying strategy to close the gaps between scientific computing
and data science. Julia unifies aspects of the scientific workflow development
process: simulation, communication, visualization, parallel data I/O, AI, and
interactive computing. This is achieved by providing (i) a dynamic just-in-
time (JIT) compiled front end to LLVM (Lattner and Adve, 2004), (ii) a
lightweight interoperability layer for existing C and Fortran HPC codes, and
(iii) a unified community ecosystem (e.g., packaging and testing). The status
quo in HPC software is to write simulation code by using a compiled language
(Fortran, C, or C++) and the Message Passing Interface (MPI) (Gropp et al.,
1999) alongside on-node standard (OpenMP (OpenMP Architecture Review Board,
2021)), vendor (CUDA (NVIDIA, 2022), HIP (AMD, 2022)) or third-party (Kokkos
(Carter Edwards et al., 2014), RAJA (Beckingsale et al., 2019)) parallel
programming models, also known as MPI+X, while data analysis and workflow
orchestration tasks are programmed by using high-level interfaces and
languages (e.g., Python)—not without caveats and trade-offs (Straßel et al.,
2020). As a result, Julia provides a valuable alternative in the convergence
of high-productivity and high-performance that must be tested on exascale
hardware.

In this work, we measure and analyze the computational performance aspects of
a Gray-Scott diffusion-reaction HPC workflow application (Pearson, 1993)
written in Julia and running on Frontier, the recently deployed exascale
system at Oak Ridge National Laboratory.111https://www.olcf.ornl.gov/frontier
We focus on the computational trade-offs of solving a representative
application of a 7-point stencil kernel on Frontier’s AMD’s MI250x GPUs at
different scales. Simultaneously, we focus on the production and consumption
of data in Jupyter Notebooks through the parallel file system. More
specifically, we present the following:

  * •

Compute measurements on AMD MI250x GPUs and comparison with a similar HIP
7-point stencil kernel

  * •

Weak scaling on up to 4,092 GPU (or 512 nodes) using MPI showing the
variability on local process time-to-solution

  * •

Parallel I/O scaling results using the ADIOS-2 Julia bindings

  * •

The overall feasibility of Julia as a single language and ecosystem for all
above components and its use on JupyterHub for data analysis of the generated
data.

The paper is organized as follows: Section 2 presents related work in the use
of Julia in HPC environments. Section 3 describes the components of the tested
workflow, including parallel simulations and data analysis. Section 4
describes the Julia implementation as a single unifying paradigm for GPU, MPI,
parallel I/O, and data analysis using Frontier’s compute nodes and the Lustre
file system with the ADIOS2 (Godoy et al., 2020) library. Section 5 shows
performance and weak scaling results when using up to 4,092 GPUs or 512 nodes
of the Frontier supercomputer, including achieved GPU bandwidths, weak
scaling, and write performance for large simulations. Conclusions and future
directions are presented in Section 6. Appendix A presents the Artifact
Description for the reproducibility of this study.

##  2\. Related Work

Recent work has reported the use of Julia in different aspects of HPC, and its
use is still an area of active exploration and community engagement (Churavy
et al., 2022).

On the compute side, Lin and McIntosh-Smith Lin and McIntosh-Smith (2021)
compared memory-bound (BabelStream) and
compute-bound (miniBUDE) benchmark mini-app ports in Julia against several
implementations on different CPU and GPU vendor systems. They concluded that
Julia’s on-node performance is either on par or slightly behind the vendor
compiled implementations for their use cases. Godoy et al. Godoy et al. (2023)
showed that Julia is competitive for a simple matrix multiplication kernel on
modern exascale nodes when compared alongside the high-level Python’s Numba
(Lam et al., 2015) JIT layer running on top of LLVM and the C++ Kokkos
programming model against vendor OpenMP on CPU and vendor CUDA and HIP
implementations on GPU. Faingnaert et al. Faingnaert et al. (2022) provide
optimized GEMM kernels in Julia that are competitive with cuBLAS and CUTLASS
implementations. Ranocha et al. Ranocha et al. (2022) describe similar levels
of performance for the Trixi.jl partial differential equation (PDE) solver at
scale when the Julia implementation is compared against compiled HPC
languages.

On the network communication side, Giordano et al. Giordano et al. (2022)
measured comparable system-level scalability performance for Julia’s MPI.jl
(Byrne et al., 2021) bindings on the Arm-based Fujitsu A64FX Fugaku system
running point-to-point communication benchmarks for different floating point
representations. Hunold et al.Hunold and Steiner (2020) used a Julia port of
the STREAM benchmark (McCalpin et al., 1995) to show negligible overheads in
the collective Bcast, Alltoall, and Allreduce MPI operations when compared
against equivalent C implementations on up to 1,152 processes. Shang et al.
Shang et al. (2022) used the Julia language ecosystem on the many-cores Sunway
supercomputer to conduct quantum computational chemistry simulations. In their
experiments they achieved up to 91% efficiency when measuring weak scaling on
up to 21M cores. Regier et al. Regier et al. (2019) provided Celeste, an early
application of Julia on multi-threaded CPU HPC systems. Celeste achieved 1.54
petaflops using 1.3 million threads on 9,300 Knights Landing (KNL) nodes of
the Cori supercomputer, while exploiting burst buffer for data I/O. More
recently, Lu et al. Lu et al. (2022) explored the integration of Julia in
their framework for moving computation as data in binary or LLVM bit code
representations. They found that their implementation can leverage JIT
optimizations, but further exploration is needed.

To the best of our knowledge, our work represents the pioneering effort to
harness the potential of Julia in modern HPC end-to-end workflows targeting a
GPU-powered exascale system.

##  3\. Workflow Description

This section describes the characteristics of the computation, communication,
parallel data I/O, and data analysis of the representative Gray-Scott
diffusion-reaction simulation code running on Frontier.

Figure 1 shows a schematic representation of a typical parallel simulation and
data analysis workflow components that apply to the current Gray-Scott use
case. Julia’s value proposition lies in its capacity to utilize a single
language, along with an integrated ecosystem, to write these components. This
unified approach allows the same language to not only operate at scale but
also to conduct data analysis in an interactive environment, thereby fostering
a seamless connection between computational execution and analytical
exploration.

Figure 1. Gray-Scott workflow components schematic illustrating the unifying
nature of Julia.

###  3.1. Simulation

Gray-Scott is a two-variable, diffusion-reaction 3D model described by the
PDEs shown in Equations (1a) and (1b):

(1a) $\displaystyle\frac{\partial U}{\partial
t}=D_{U}\nabla^{2}U-UV^{2}+F\left(1-U\right)+nr$ (1b)
$\displaystyle\frac{\partial V}{\partial
t}=D_{V}\nabla^{2}V+UV^{2}+-\left(F+k\right)V$

where $U$ and $V$ are the output concentrations of two reacting and diffusing
chemicals, while the inputs are listed as follows:

  * •

$D_{u}$ and $D_{v}$ are the diffusion rates for $U$ and $V$;

  * •

$F$ is the feed rate of $U$ into the system;

  * •

$k$ is the kill rate of $V$ from the system;

  * •

$n$ is the magnitude of the noise to be added to the system; and

  * •

$r$ is a uniformly distributed random number between $-$1 and 1 for each time
and spatial coordinate.

As illustrated in Equations (2a) and (2b), the set of governing equations is
discretized in time, $t$, and space, $i,j,k$, on a regular normalized mesh by
using simple forward and central differences, respectively:

(2a) $\displaystyle U^{t+1}_{i,j,k}$ $\displaystyle=$ $\displaystyle
U^{t}_{i,j,k}+\Delta t\left[D_{U}\nabla^{2}U^{t}_{i,j,k}+S^{t}_{U}\right]$
(2b) $\displaystyle V^{t+1}_{i,j,k}$ $\displaystyle=$ $\displaystyle
V^{t}_{i,j,k}+\Delta t\left[D_{V}\nabla^{2}V^{t}_{i,j,k}+S^{t}_{V}\right]$

where $\Delta t$ is an input time step variable, $S$ represents the local
source terms for $U$ and $V$ defined in Equations (1a) and (1b), and the
Laplacian operator, $\nabla^{2}$, is defined in Equation (3) for the 3D
nearest-neighbor Jacobi 7-point stencil in normalized spatial units:

$\displaystyle\nabla^{2}U^{t}_{i,j,k}=-\,U^{t}_{i,j,k}+\frac{1}{6}\left[U^{t}_{i-1,j,k}+U^{t}_{i+1,j,k}+\right.$
$\displaystyle U^{t}_{i,j-1,k}+U^{t}_{i,j+1,k}+$ (3)
$\displaystyle\left.U^{t}_{i,j,k-1}+U^{t}_{i,j,k+1}\right]$

The random nature of the source term caused by $r$ and the coupled nature of
$U$ and $V$ provide slightly more complexity than a typical diffusion
calculation (e.g., heat transfer equation), as shown in Figure 2 for the
visualization of the solution of the field variable $U$ at a center plane over
a few time iterations.

Figure 2. Gray-Scott snapshot of the solution of the 3D field $V$ variable
showing the random nature of the source terms in the PDE.

###  3.2. On-Node Parallelization

The stencil calculation is the main computational kernel in Gray-Scott.
Stencil calculations are an important class of PDE solvers on structured grids
(Datta et al., 2008), and they have been successfully optimized for different
hardware architectures, including GPUs (Holewinski et al., 2012; Krotkiewski
and Dabrowski, 2013; Godoy and Liu, 2011).

Because each new time step value for $U$ and $V$ only depends on neighboring
values, the problem is parallelizable on GPUs—and CPUs—thereby providing an
extra memory allocation for temporary $U_{t+1}$ and $V_{t+1}$ field values for
the next time iteration. Figure 3 illustrates the underlying contiguous memory
allocated for each variable, which is accessed from near and far away points
in memory; hence, some cache misses are expected in this typical kernel. Each
variable requires 7 reads (fetch) $+$ 1 write $=$ 8 floating point access
operations per variable for a total of 16 per cell in solving Equations (2a)
and (2b).

Figure 3. Seven-point stencil memory access for one variable in the parallel
Gray-Scott solver.

###  3.3. MPI Communication

The 3D regular domain shown in Figure 2 is decomposed using an MPI Cartesian
communicator. Each subdomain source shares face cell nodes with their
neighboring subdomain destination at every computational iteration. This is
illustrated in Figure 4, which shows a typical MPI_Send/MPI_Recv pattern with
ghost cell surfaces. Because the exchanged surfaces are not memory contiguous,
we define a new strided vector type by using MPI_Datatypes and MPI_Type_vector
(Wu et al., 2016). We keep the MPI exchange from CPU-allocated memory to
populate the strided vector contents coming from the GPU. We did not
experiment with GPU-aware MPI for this work.

Figure 4. MPI send/receive memory patterns in a Cartesian communicator
decomposition used in Gray-Scott for each variable $U$ and $V$.

###  3.4. Parallel I/O and Data Analysis

Data production and consumption is performed to enable post-processing (e.g.,
visualization) of the solution of the $U$ and $V$ field variables over time.
As such, a simple visualization schema and the corresponding metadata
accounting for the input provenance are stored in the resulting dataset. Due
to the volume of data produced from the simulation, drastically reducing the
frequency of writes to the parallel file system is often required. Listing 1
presents a sample provenance record of the data produced from the Gray-Scott
simulation.

Listing 1: Provenance of the data generated from the Gray-Scott simulation.

⬇

double Du attr = 0.2

double Dv attr = 0.1

double F attr = 0.02

double U 1000*{1024, 1024, 1024}

Min/Max -0.120795 / 1.46671

double V 1000*{1024, 1024, 1024}

Min/Max 0 / 0.959875

double dt attr = 1

double k attr = 0.048

double noise attr = 0.1

int32_t step 50*scalar = 20 / 1000

Attribute visualization schemas: FIDES, VTX

##  4\. Implementation

We developed GrayScott.jl as a comprehensive HPC workflow application packaged
as a Julia module. The resulting Julia implementation using the AMDGPU.jl
(Samaroo et al., 2022) package is shown in Listing 2. This implementation is a
straightforward application of a typical GPU kernel and maps each component of
the kernel onto a 3D computational grid. When launching this kernel and
allocating the number of threads and grid size, one must consider that Julia
arrays are column-major. Therefore, the fastest index, being the first one,
should be structured to avoid splitting across threads on the GPU.

Although the goal is not to increase the efficiency of the present stencil
calculation, we aim to evaluate how the performance of the JIT generated Julia
LLVM-IR stencil code compares with theoretical bandwidths and a simple AMD-
provided Laplacian kernel. The metrics for this implementation are described
and provided in Section 5.

Listing 2: Julia AMDGPU.jl Gray-Scott kernel.

⬇

using AMDGPU

using Distributions

function _laplacian(i, j, k, var)

l = var[i - 1, j, k] + var[i + 1, j, k]

+ var[i, j - 1, k] + var[i, j + 1, k]

+ var[i, j, k - 1] + var[i, j, k + 1]

- 6.0 * var[i, j, k]

return l / 6.0

end

function _kernel_amdgpu!(u, v, u_temp, v_temp,

sizes, Du, Dv, F, K,

noise, dt)

k = (workgroupIdx().x - 1) * workgroupDim().x

+ workitemIdx().x

j = (workgroupIdx().y - 1) * workgroupDim().y

+ workitemIdx().y

i = (workgroupIdx().z - 1) * workgroupDim().z

+ workitemIdx().z

if k == 1 || k >= sizes[3] ||

j == 1 || j >= sizes[2] ||

i == 1 || i >= sizes[1]

return

end

@inbounds begin

u_ijk = u[i, j, k]

v_ijk = v[i, j, k]

du = Du * _laplacian(i, j, k, u)

- u_ijk * v_ijk^2 + F * (1.0 - u_ijk)

+ noise * rand(Uniform(-1, 1))

dv = Dv * _laplacian(i, j, k, v)

+ u_ijk * v_ijk^2 - (F + K) * v_ijk

u_temp[i, j, k] = u_ijk + du * dt

v_temp[i, j, k] = v_ijk + dv * dt

end

return nothing

end

For managing communication, we use the MPI.jl package, which provides a high-
level binding interface with the underlying MPI back end. The vector data
types and MPI send/receive communication code for a single variable ghost cell
surface are described in Listing 3. MPI.jl uniquely enables referencing the
internal memory through @view to form a buffer structure prior to an actual
send/receive operation. We purposely introduce variations in the amount of
data shared across the $x$,$y$, and $z$ dimensions to capture distinct
communication patterns.

Listing 3: GrayScott.jl MPI.jl implementation

⬇

# if in GPU memory move first to CPU

u = typeof(fields.u) <: Array ? fields.u :

Array(fields.u)

xy_face_t = MPI.Types.create_vector(size_y + 2,

size_x, size_x + 2, MPI.Datatype(T))

MPI.Types.commit!(xy_face_t)

...

send = MPI.Buffer(@view(u[2, 1, size_z + 1]),

1, data_type)

recv = MPI.Buffer(@view(u[2, 1, 1]),

1, data_type)

MPI.Sendrecv!(send, recv, comm, dest = rank1,

source = rank2)

We use the ADIOS2 library via the Julia ADIOS2.jl bindings. ADIOS2 is an MPI-
based library, and the variables are described by accounting for their
partitioning in a parallel I/O in which each MPI process produces a block of
data in the global $U$ and $V$ 3D array variables. A typical ADIOS2 dataset is
shown in Listing 1 in which the combination of all data producer blocks from
the MPI subdomain form the global $U$ and $V$ arrays from the PDE solution. We
add schema attributes to enable later visualization in ParaView (Ayachit,
2015) by using FIDES (Pugmire et al., 2021) and VTX readers implemented in VTK
(Schroeder et al., 2000).

###  4.1. The Frontier Supercomputer

Frontier is equipped with AMD CPUs and GPUs and achieved 1.194 EFLOPS in the
HPL benchmark (Dongarra et al., 2003) according to the June 2023 TOP500
list.222https://www.top500.org/lists/top500/2023/06/ Table 1 lists Frontier
system characteristics and the software stack used for this study. In
structuring our numerical experiments, we have designed them to grow by a
factor of 8, so each dimension will double ($2^{3}=8$), thereby aligning with
the GPU’s maximum limits of 1,024 threads per dimension in a 3D computational
kernel placement (Listing 2).

Table 1. Summary of Frontier hardware and software used in this study.
Frontier Characteristics
---
Nodes | 9,408
—CPU Architecture—
CPU | AMD EPYC 7A53
Cores | 64
Memory | DDR4 512 GB
Bandwidth | 205 GB/s
—GPU Architecture—
GPU |  4$\times$AMD MI250X
|  8$\times$GCDs
Memory | HBM2E 64 GB
Bandwidth | 1,600 GB/s per GCD
—Connectivity—
GPU-to-GPU | Infinity Fabric
Bandwidth | 50–100 GB/s
GPU-to-CPU | Infinity Fabric
Bandwidth | 36 GB/s
—File System—
Type | Lustre Orion
Capacity | 679 PB
Nodes | 40 metadata
| 450 object storage service (OSS)
Write speed | 5.5 TBps
Read speed | 4.5 TBps
—Software versions—
Julia | 1.9.2
AMDGPU.jl | 0.4.15
ROCm | 5.4.0
MPI.jl | 0.20.12
Cray-MPICH | 8.1.23
ADIOS2.jl | 1.2.1
ADIOS2 | 2.8.3

##  5\. Results and Discussion

In this section, we present the results from our evaluation of the Gray-Scott
HPC workflow application regarding (i) compute capabilities on Frontier’s
MI250x AMD GPU using AMDGPU.jl, (ii) weak scalability of the simulation wall-
clock times using Julia’s MPI.jl, and (iii) simulation data production and
consumption rates in JupyterHub using Julia’s bindings to ADIOS2,
ADIOS2.jl.333https://github.com/eschnett/ADIOS2.jl

Figure 5. Gray-Scott simulation trace obtained with rocprof showing
computational load on GPU and memory transfer to CPU for communication.

###  5.1. Computation on MI250x GPUs

As explained in Section 3.1, the metrics for the memory-bound Gray-Scott
kernel presented in Listing 2 are based on the performance of accessing the
GPU allocated arrays $u$ and $v$ for 7 read (fetch/load) operations and 1
write operation. We collect two types of bandwidth information: (i) effective
bandwidth, meaning the expected minimal data movement (fetch/load and
write/store) according to the number of operations in a single domain, and
(ii) the total measured bandwidth obtained from AMD’s rocprof profiler for
FETCH_SIZE and WRITE_SIZE, which represent the total number of KB moved to and
from the GPU memory (Figure 5).

The effective fetch and write sizes can be estimated from the number of
floating point operations per cell by considering that edge cells execute
fewer operations due to reduced stencil size. A regular grid with $L$ number
of cells per direction is provided in Equations (4a) and (4b), as reported in
AMD’s training material.444https://github.com/AMD/amd-lab-notes

(4a) $\displaystyle
fetch\\_size_{effective}=\left[L^{3}-8-12\left(L-2\right)\right]\cdot
sizeof(T)$ (4b) $\displaystyle write\\_size_{effective}=(L-2)^{3}\cdot
sizeof(T)$

Hence, the effective and total bandwidth metrics are shown in Equations (5a)
and (5b) for the effective (theoretical) and the actual (measured with
rocprof) GPU data access bandwidths.

(5a) $\displaystyle
bandwidth_{effective}=\frac{\left(fetch\\_size+write\\_size\right)_{effective}}{kernel\\_time}$
(5b) $\displaystyle
bandwidth_{total}=\frac{\left(FETCH\\_SIZE+WRITE\\_SIZE\right)_{rocprof}}{kernel\\_time}$

Table 2. Average bandwidth comparison of different stencil implementations on a single GPU. Kernel | Bandwidth (GB/s)
---|---
| Effective | Total
Julia GrayScott.jl |  |
\- 2-variable (application) | 312 | 570
\- 1-variable no random | 312 | 625
HIP single variable | 599 | 1,163
Theoretical peak MI250x | 1,600
Table 3. The rocprof outputs for HIP 1-variable and Julia Gray-Scott.jl implementations. Kernel | HIP | GrayScott.jl
---|---|---
metric | 1-var | 1-variable | 2-variable
|  | no random | (application)
wgr | 256 | 512 | 512
lds | 0 | 29,184 | 29,184
scr | 0 | 8,192 | 8,192
FETCH_SIZE (GB) | 25.08 | 25.40 | 50.80
WRITE_SIZE (GB) | 8.35 | 8.38 | 16.78
TCC_HIT (M) | 9.14 | 10.80 | 24.60
TCC_MISS (M) | 8.36 | 8.69 | 17.19
Avg Duration (ms) | 28.74 | 54.03 | 111.07

Listing 4 presents the generated LLVM-IR memory accesses in the Julia
implementation. These numbers are consistent with the algorithm patterns in
Listing 2 of 16 loads and 2 stores. Consequently, we did not observe
additional overhead typical of high-level implementations. The reasons for the
bandwidth performance difference in Table 3 are beyond the IR level and might
present opportunities for lower-level optimizations at the LLVM instruction
set architecture level, thereby providing an opportunity to close gaps from
vendor-specific optimizations as Julia’s AMDGPU.jl and the AMD HIP/ROCm stack
continues to rapidly evolve.

Listing 4: GrayScott.jl application kernel unique memory loads (14) and store
(2) in LLVM-IR.

⬇

%94 = load double, double addrspace(1)* %93, align 8

%103 = load double, double addrspace(1)* %102, align 8

%107 = load double, double addrspace(1)* %106, align 8

%110 = load double, double addrspace(1)* %109, align 8

%114 = load double, double addrspace(1)* %113, align 8

%117 = load double, double addrspace(1)* %116, align 8

%122 = load double, double addrspace(1)* %121, align 8

%126 = load double, double addrspace(1)* %125, align 8

%312 = load double, double addrspace(1)* %311, align 8

%315 = load double, double addrspace(1)* %314, align 8

%318 = load double, double addrspace(1)* %317, align 8

%321 = load double, double addrspace(1)* %320, align 8

%325 = load double, double addrspace(1)* %324, align 8

%329 = load double, double addrspace(1)* %328, align 8

...

store double %345, double addrspace(1)* %353, align 8

store double %355, double addrspace(1)* %363, align 8

###  5.2. Weak scaling

We ran GrayScott.jl with different cell sizes using 1 GCD (GPU) per MPI
process. Weak scaling results are illustrated in Figure 6 and show the total
wall-clock times per MPI process for a constant workload of $nx=ny=nz=1,024$
per dimension, or $1,024^{3}=1,073\times 10^{9}$ cells per GPU.

Figure 6. Weak scaling with single MPI process wall-clock time variability
obtained with Gray-Scott.jl on Frontier.

Per-process wall-clock time variability in HPC workflows is proportional to
the number of MPI processes. As expected, the overall communication overhead
is dictated by the slowest time-to-solution processes (dashed blue line). The
Julia code could run successfully on 4,096 GPUs or 512 nodes of Frontier
(5.44% of the system). A clear trend of small variability of a few percent
(2%–3%) is apparent up to 512 MPI processes, whereas large variability
(12%–15%) appears for the largest case. Although we do not have a baseline
comparison, this initial measurement can help with further optimizations
because Julia can produce realistic scenarios for MPI communication on a
system such as Frontier. When attempting runs for the next size for a factor 8
at 32,768 GPUs (4,096 nodes, 50% of the Frontier system), unpredictable
failures occurred at the underlying MPI layers during the ghost cell exchange
stage. Nevertheless, all 32,768 GPUs showed initial runs keeping the bandwidth
close to the expected value of 312 GB/s, as reported in Table 2.

In the context of HPC workflows, understanding the interplay between scale and
performance is critical. To this end, we delve deeper into GPU bandwidth
metrics from large-scale runs. Figure 7 shows the obtained bandwidth
distributions for the initial JIT compiled run, which represents the initial
overhead in running Julia kernels, and the optimized code when running on
4,096 GPUs. The initial JIT compilation represents the early overhead
associated with running Julia kernels in an HPC workflow—an
$\approx$12.5$\times$ cost, or 8% on average, of the bandwidth in the
optimized kernel. Although JIT compilation represents an amortized cost, it
should be incorporated in the estimation of total wall-clock run times. While
this overhead can be addressed in HPC workflow planning, Julia’s ahead-of-time
mechanism was not explored in this study. These insights not only demonstrate
Julia’s capacity for HPC workflows but also reveal areas for further research
and optimization, which could enable more efficient utilization of resources
in complex computational environments.

Figure 7. Bandwidth distribution for JIT and optimized kernel code running on
4,096 Frontier GPUs for 20 simulation steps.

###  5.3. Parallel I/O

In HPC workflows, the efficient use of parallel file systems is crucial for
both saving computational results and ensuring smooth data transfer between
different stages of a workflow. To this end, we evaluate using the parallel
file system within the GrayScott.jl package, focusing on saving one output
step of the cases presented in Figure 6. Results from using the ADIOS2.jl
Julia bindings package are presented in Figure 8, revealing the performance
characteristics of weak scaling in the writing of the dataset described in
Section 4. The underlying ADIOS2 library uses the default BP5 format, which
generates a single sub-file per node. At 512 nodes, write times are fairly
flat, which represents an increase of write bandwidth up to 434 GB/s for the
largest case. This represents an 8% peak file system write performance of 5.5
TB/s when using 5% of the machine (Table 1). In the context of HPC workflows,
such a performance metric is significant, as parallel I/O performance is
dictated by low-level POSIX write operations, the ADIOS2 buffering management,
and real-time file system usage. With the Julia binding, we observed
negligible overhead, indicating that Julia can serve as an effective tool in
HPC workflow optimization.

Figure 8. Weak scaling on parallel I/O showing wall-clock times and bandwidth
performance when using ADIOS2.jl on Frontier.

Finally, we look at Julia in an end-to-end workflow development process for
data analysis. We tested a Julia implementation by reading a plotting of
simple 2D slices of the stored $U$ and $V$ 3D arrays with the Makie.jl library
(Danisch and Krumbiegel, 2021). Figure 9 shows a snapshot of the resulting
JupyterHub service consuming the ADIOS2-generated data shared on Frontier’s
file system. This instance demonstrates Julia’s potential as a unified
language that can seamlessly connect different stages of an HPC workflow, from
computational simulation to data visualization. Future work includes further
performance tuning and evaluation of the trade-offs for in-memory streaming
data pipelines (Poeschel et al., 2022).

Figure 9. Snapshot of Julia data analysis implementation on the OLCF’s
JupyterHub server.

By bridging various components in a typical HPC workflow, Julia could emerge
as a compelling option for scientists and engineers working on complex, large-
scale computational problems. It also sets the stage for more integrated,
responsive, and adaptive HPC workflows that can fully leverage the
capabilities of modern supercomputing infrastructures.

##  6\. Conclusions

On up to 4,096 GPUs and MPI processes (512 nodes of the Frontier
supercomputer), we used a 2-variable diffusion-reaction code, Gray-Scott, to
test the performance of the Julia HPC ecosystem in the development of tightly
coupled workflow components. Overall, the Julia stencil solver achieves close
to a 50% bandwidth on Frontier’s MI250x AMD GPUs vs. a HIP implementation;
hence, performance gaps still exist and must be closed as we look forward to
future versions of the actively developed AMDGPU.jl. Meanwhile, the measured
weak scaling from the MPI communication and parallel I/O components suggest
that bindings available in Julia are lightweight and perfectly suitable
running on top of the underlying MPI and ADIOS2 library implementations. The
Julia implementation shows similar patterns in overhead and variability
typical of network and file system communication in HPC systems. We also
showcased the data analysis feasibility on available JupyterHub systems.
Overall, the LLVM-based Julia HPC ecosystem presents an attractive alternative
for developing co-design components given the high-performance and high-
productivity requirements for the end-to-end workflows that power scientific
discovery (e.g., AI, FAIR).

###### Acknowledgements.

This research was supported by the Exascale Computing Project (17-SC-20-SC), a
collaborative effort of the US Department of Energy Office of Science and the
National Nuclear Security Administration. This research used resources of the
Oak Ridge Leadership Computing Facility and the Experimental Computing
Laboratory (ExCL) at the Oak Ridge National Laboratory, which is supported by
the Office of Science of the US Department of Energy under Contract No. DE-
AC05-00OR22725. This work is funded, in part, by Bluestone, a X-Stack project
in the DOE Advanced Scientific Computing Office with program manager Hal
Finkel.

## References

  * (1)
  * AMD (2022) AMD. 2022.  AMD ROCm v5.2 Release.  https://rocmdocs.amd.com/en/latest/Current_Release_Notes/Current-Release-Notes.html#amd-rocm-v5-2-release
  * Ayachit (2015) Utkarsh Ayachit. 2015\.  _The ParaView guide: a parallel visualization application_.  Kitware, Inc.
  * Beckingsale et al. (2019) David A. Beckingsale, Jason Burmark, Rich Hornung, Holger Jones, William Killian, Adam J. Kunen, Olga Pearce, Peter Robinson, Brian S. Ryujin, and Thomas RW Scogland. 2019.  RAJA: Portable Performance for Large-Scale Scientific Applications. In _2019 IEEE/ACM International Workshop on Performance, Portability and Productivity in HPC (P3HPC)_. 71–81.  https://doi.org/10.1109/P3HPC49587.2019.00012
  * Ben-Nun et al. (2020) Tal Ben-Nun, Todd Gamblin, D. S. Hollman, Hari Krishnan, and Chris J. Newburn. 2020.  Workflows are the New Applications: Challenges in Performance, Portability, and Productivity. In _2020 IEEE/ACM International Workshop on Performance, Portability and Productivity in HPC (P3HPC)_. 57–69.  https://doi.org/10.1109/P3HPC51967.2020.00011
  * Bezanson et al. (2017) Jeff Bezanson, Alan Edelman, Stefan Karpinski, and Viral B Shah. 2017.  Julia: A Fresh Approach to Numerical Computing.  _SIAM Rev._ 59, 1 (Jan. 2017), 65–98.  https://doi.org/10.1137/141000671
  * Byrne et al. (2021) Simon Byrne, Lucas C. Wilcox, and Valentin Churavy. 2021\.  MPI.jl: Julia bindings for the Message Passing Interface.  _Proceedings of the JuliaCon Conferences_ 1, 1 (2021), 68\.  https://doi.org/10.21105/jcon.00068
  * Carter Edwards et al. (2014) H. Carter Edwards, Christian R. Trott, and Daniel Sunderland. 2014.  Kokkos: Enabling manycore performance portability through polymorphic memory access patterns.  _J. Parallel and Distrib. Comput._ 74, 12 (2014), 3202–3216.  https://doi.org/10.1016/j.jpdc.2014.07.003 Domain-Specific Languages and High-Level Frameworks for High-Performance Computing.
  * Churavy et al. (2022) Valentin Churavy, William F Godoy, Carsten Bauer, Hendrik Ranocha, Michael Schlottke-Lakemper, Ludovic Räss, Johannes Blaschke, Mosè Giordano, Erik Schnetter, Samuel Omlin, Jeffrey S. Vetter, and Alan Edelman. 2022\.  Bridging HPC Communities through the Julia Programming Language.  arXiv:2211.02740 [cs.DC]
  * Danisch and Krumbiegel (2021) Simon Danisch and Julius Krumbiegel. 2021.  Makie.jl: Flexible high-performance data visualization for Julia.  _Journal of Open Source Software_ 6, 65 (2021), 3349\.  https://doi.org/10.21105/joss.03349
  * Datta et al. (2008) Kaushik Datta, Mark Murphy, Vasily Volkov, Samuel Williams, Jonathan Carter, Leonid Oliker, David Patterson, John Shalf, and Katherine Yelick. 2008.  Stencil computation optimization and auto-tuning on state-of-the-art multicore architectures. In _SC ’08: Proceedings of the 2008 ACM/IEEE Conference on Supercomputing_. 1–12.  https://doi.org/10.1109/SC.2008.5222004
  * Deelman et al. (2018) Ewa Deelman, Tom Peterka, Ilkay Altintas, Christopher D Carothers, Kerstin Kleese van Dam, Kenneth Moreland, Manish Parashar, Lavanya Ramakrishnan, Michela Taufer, and Jeffrey Vetter. 2018\.  The future of scientific workflows.  _The International Journal of High Performance Computing Applications_ 32, 1 (2018), 159–175.  https://doi.org/10.1177/1094342017704893 arXiv:https://doi.org/10.1177/1094342017704893
  * Dongarra et al. (2003) Jack J Dongarra, Piotr Luszczek, and Antoine Petitet. 2003\.  The LINPACK benchmark: past, present and future.  _Concurrency and Computation: practice and experience_ 15, 9 (2003), 803–820.
  * Faingnaert et al. (2022) Thomas Faingnaert, Tim Besard, and Bjorn De Sutter. 2022\.  Flexible Performant GEMM Kernels on GPUs.  _IEEE Transactions on Parallel and Distributed Systems_ 33, 9 (2022), 2230–2248.  https://doi.org/10.1109/TPDS.2021.3136457
  * Ferreira da Silva et al. (2023) Rafael Ferreira da Silva, Rosa M. Badia, Venkat Bala, Debbie Bard, Timo Bremer, Ian Buckley, Silvina Caino-Lores, Kyle Chard, Carole Goble, Shantenu Jha, Daniel S. Katz, Daniel Laney, Manish Parashar, Frederic Suter, Nick Tyler, Thomas Uram, Ilkay Altintas, et al. 2023\.  _Workflows Community Summit 2022: A Roadmap Revolution_.  Technical Report ORNL/TM-2023/2885. Oak Ridge National Laboratory.  https://doi.org/10.5281/zenodo.7750670
  * Ferreira da Silva et al. (2021) Rafael Ferreira da Silva, Henri Casanova, Kyle Chard, Ilkay Altintas, Rosa M Badia, Bartosz Balis, Tainã Coleman, Frederik Coppens, Frank Di Natale, Bjoern Enders, Thomas Fahringer, Rosa Filgueira, Grigori Fursin, Daniel Garijo, Carole Goble, Dorran Howell, Shantenu Jha, Daniel S. Katz, Daniel Laney, Ulf Leser, Maciej Malawski, Kshitij Mehta, Loïc Pottier, Jonathan Ozik, J. Luc Peterson, Lavanya Ramakrishnan, Stian Soiland-Reyes, Douglas Thain, and Matthew Wolf. 2021.  A Community Roadmap for Scientific Workflows Research and Development. In _2021 IEEE Workshop on Workflows in Support of Large-Scale Science (WORKS)_. 81–90.  https://doi.org/10.1109/WORKS54523.2021.00016
  * Ferreira da Silva et al. (2017) Rafael Ferreira da Silva, Rosa Filgueira, Ilia Pietri, Ming Jiang, Rizos Sakellariou, and Ewa Deelman. 2017.  A Characterization of Workflow Management Systems for Extreme-Scale Applications.  _Future Generation Computer Systems_ 75 (2017), 228–238.  https://doi.org/10.1016/j.future.2017.02.026
  * Giordano et al. (2022) Mosè Giordano, Milan Klöwer, and Valentin Churavy. 2022\.  Productivity meets Performance: Julia on A64FX. In _2022 IEEE International Conference on Cluster Computing (CLUSTER)_. 549–555.  https://doi.org/10.1109/CLUSTER51413.2022.00072
  * Goble et al. (2020) Carole Goble, Sarah Cohen-Boulakia, Stian Soiland-Reyes, Daniel Garijo, Yolanda Gil, Michael R Crusoe, Kristian Peters, and Daniel Schober. 2020.  FAIR computational workflows.  _Data Intelligence_ 2, 1-2 (2020), 108–121.
  * Godoy and Liu (2011) William F. Godoy and Xu Liu. 2011.  Introduction of Parallel GPGPU Acceleration Algorithms for the Solution of Radiative Transfer.  _Numerical Heat Transfer, Part B: Fundamentals_ 59, 1 (2011), 1–25.  https://doi.org/10.1080/10407790.2010.541359
  * Godoy et al. (2020) William F. Godoy, Norbert Podhorszki, Ruonan Wang, Chuck Atkins, Greg Eisenhauer, Junmin Gu, Philip Davis, Jong Choi, Kai Germaschewski, Kevin Huck, Axel Huebl, Mark Kim, James Kress, Tahsin Kurc, Qing Liu, Jeremy Logan, Kshitij Mehta, George Ostrouchov, Manish Parashar, Franz Poeschel, David Pugmire, Eric Suchyta, Keichi Takahashi, Nick Thompson, Seiji Tsutsumi, Lipeng Wan, Matthew Wolf, Kesheng Wu, and Scott Klasky. 2020.  Adios 2: The Adaptable Input Output System. A framework for high-performance data management.  _SoftwareX_ 12 (2020), 100561.  https://doi.org/10.1016/j.softx.2020.100561
  * Godoy et al. (2023) William F. Godoy, Pedro Valero-Lara, T. Elise Dettling, Christian Trefftz, Ian Jorquera, Thomas Sheehy, Ross G. Miller, Marc Gonzalez-Tallada, Jeffrey S. Vetter, and Valentin Churavy. 2023.  Evaluating performance and portability of high-level programming models: Julia, Python/Numba, and Kokkos on exascale nodes. In _2023 IEEE International Parallel and Distributed Processing Symposium Workshops (IPDPSW)_. 373–382.  https://doi.org/10.1109/IPDPSW59300.2023.00068
  * Gropp et al. (1999) William Gropp, Rajeev Thakur, and Ewing Lusk. 1999\.  _Using MPI-2: Advanced features of the message passing interface_.  MIT press.
  * Holewinski et al. (2012) Justin Holewinski, Louis-Noël Pouchet, and P. Sadayappan. 2012.  High-Performance Code Generation for Stencil Computations on GPU Architectures. In _Proceedings of the 26th ACM International Conference on Supercomputing_ (San Servolo Island, Venice, Italy) _(ICS ’12)_. Association for Computing Machinery, New York, NY, USA, 311–320.  https://doi.org/10.1145/2304576.2304619
  * Hunold and Steiner (2020) Sascha Hunold and Sebastian Steiner. 2020.  Benchmarking Julia’s Communication Performance: Is Julia HPC ready or Full HPC?. In _2020 IEEE/ACM Performance Modeling, Benchmarking and Simulation of High Performance Computer Systems (PMBS)_. IEEE.  https://doi.org/10.1109/pmbs51919.2020.00008
  * Krotkiewski and Dabrowski (2013) Marcin Krotkiewski and Marcin Dabrowski. 2013.  Efficient 3D stencil computations using CUDA.  _Parallel Comput._ 39, 10 (2013), 533–548.  https://doi.org/10.1016/j.parco.2013.08.002
  * Lam et al. (2015) Siu Kwan Lam, Antoine Pitrou, and Stanley Seibert. 2015\.  Numba: A LLVM-based Python JIT compiler. In _Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC_. 1–6.
  * Lattner and Adve (2004) Chris Lattner and Vikram Adve. 2004.  LLVM: A compilation framework for lifelong program analysis & transformation. In _International Symposium on Code Generation and Optimization, 2004. CGO 2004._ IEEE, 75–86.
  * Lin and McIntosh-Smith (2021) Wei-Chen Lin and Simon McIntosh-Smith. 2021.  Comparing Julia to Performance Portable Parallel Programming Models for HPC. In _2021 International Workshop on Performance Modeling, Benchmarking and Simulation of High Performance Computer Systems (PMBS)_. 94–105.  https://doi.org/10.1109/PMBS54543.2021.00016
  * Lu et al. (2022) Wenbin Lu, Luis E. Peña, Pavel Shamis, Valentin Churavy, Barbara Chapman, and Steve Poole. 2022\.  Bring the BitCODE-Moving Compute and Data in Distributed Heterogeneous Systems. In _2022 IEEE International Conference on Cluster Computing (CLUSTER)_. 12–22.  https://doi.org/10.1109/CLUSTER51413.2022.00017
  * McCalpin et al. (1995) John D McCalpin et al. 1995\.  Memory bandwidth and machine balance in current high performance computers.  _IEEE computer society technical committee on computer architecture (TCCA) newsletter_ 2, 19-25 (1995).
  * NVIDIA (2022) NVIDIA. 2022.  CUDA Toolkit Documentation - v11.7.0.  https://developer.nvidia.com/cuda-toolkit
  * OpenMP Architecture Review Board (2021) OpenMP Architecture Review Board. 2021\.  OpenMP Application Program Interface Version 5.2.  https://www.openmp.org/wp-content/uploads/OpenMP-API-Specification-5-2.pdf
  * Pearson (1993) John E. Pearson. 1993\.  Complex Patterns in a Simple System.  _Science_ 261, 5118 (1993), 189–192.  https://doi.org/10.1126/science.261.5118.189 arXiv:https://www.science.org/doi/pdf/10.1126/science.261.5118.189
  * Poeschel et al. (2022) Franz Poeschel, Juncheng E, William F. Godoy, Norbert Podhorszki, Scott Klasky, Greg Eisenhauer, Philip E. Davis, Lipeng Wan, Ana Gainaru, Junmin Gu, Fabian Koller, René Widera, Michael Bussmann, and Axel Huebl. 2022\.  Transitioning from File-Based HPC Workflows to Streaming Data Pipelines with openPMD and ADIOS2. In _Driving Scientific and Engineering Discoveries Through the Integration of Experiment, Big Data, and Modeling and Simulation_ , Jeffrey Nichols, Arthur ‘Barney’ Maccabe, James Nutaro, Swaroop Pophale, Pravallika Devineni, Theresa Ahearn, and Becky Verastegui (Eds.). Springer International Publishing, Cham, 99–118.
  * Pugmire et al. (2021) David Pugmire, Caitlin Ross, Nicholas Thompson, James Kress, Chuck Atkins, Scott Klasky, and Berk Geveci. 2021.  Fides: a general purpose data model library for streaming data. In _High Performance Computing: ISC High Performance Digital 2021 International Workshops, Frankfurt am Main, Germany, June 24–July 2, 2021, Revised Selected Papers 36_. Springer, 495–507.
  * Ranocha et al. (2022) Hendrik Ranocha, Michael Schlottke-Lakemper, Andrew R. Winters, Erik Faulhaber, Jesse Chan, and Gregor J. Gassner. 2022.  Adaptive numerical simulations with Trixi.jl: A case study of Julia for scientific computing.  _Proceedings of the JuliaCon Conferences_ 1, 1 (2022), 77\.  https://doi.org/10.21105/jcon.00077
  * Regier et al. (2019) Jeffrey Regier, Keno Fischer, Kiran Pamnany, Andreas Noack, Jarrett Revels, Maximilian Lam, Steve Howard, Ryan Giordano, David Schlegel, Jon McAuliffe, Rollin Thomas, and Prabhat. 2019\.  Cataloging the visible universe through Bayesian inference in Julia at petascale.  _J. Parallel and Distrib. Comput._ 127 (2019), 89–104.  https://doi.org/10.1016/j.jpdc.2018.12.008
  * Samaroo et al. (2022) Julian Samaroo, Valentin Churavy, Wiktor Phillips, Ali Ramadhan, Jason Barmparesos, Julia TagBot, Ludovic Räss, Michel Schanen, Tim Besard, Anton Smirnov, Takafumi Arakaki, Stephan Antholzer, Alessandro, Chris Elrod, Matin Raayai, and Tom Hu. 2022\.  _JuliaGPU/AMDGPU.jl: v0.4.1_.  https://doi.org/10.5281/zenodo.6949520
  * Schroeder et al. (2000) William J Schroeder, Lisa Sobierajski Avila, and William Hoffman. 2000.  Visualizing with VTK: a tutorial.  _IEEE Computer graphics and applications_ 20, 5 (2000), 20–27.
  * Shang et al. (2022) Honghui Shang, Li Shen, Yi Fan, Zhiqian Xu, Chu Guo, Jie Liu, Wenhao Zhou, Huan Ma, Rongfen Lin, Yuling Yang, Fang Li, Zhuoya Wang, Yunquan Zhang, and Zhenyu Li. 2022.  Large-Scale Simulation of Quantum Computational Chemistry on a New Sunway Supercomputer. In _SC22: International Conference for High Performance Computing, Networking, Storage and Analysis_. 1–14.  https://doi.org/10.1109/SC41404.2022.00019
  * Sirvent et al. (2022) Raül Sirvent, Javier Conejero, Francesc Lordan, Jorge Ejarque, Laura Rodríguez-Navas, José M. Fernández, Salvador Capella-Gutiérrez, and Rosa M. Badia. 2022.  Automatic, Efficient and Scalable Provenance Registration for FAIR HPC Workflows. In _2022 IEEE/ACM Workshop on Workflows in Support of Large-Scale Science (WORKS)_. 1–9.  https://doi.org/10.1109/WORKS56498.2022.00006
  * Straßel et al. (2020) Dominik Straßel, Philipp Reusch, and Janis Keuper. 2020\.  Python Workflows on HPC Systems. In _2020 IEEE/ACM 9th Workshop on Python for High-Performance and Scientific Computing (PyHPC)_. 32–40.  https://doi.org/10.1109/PyHPC51966.2020.00009
  * Vetter et al. (2018) Jeffrey S. Vetter, Ron Brightwell, Maya Gokhale, Pat McCormick, Rob Ross, John Shalf, Katie Antypas, David Donofrio, Travis Humble, Catherine Schuman, Brian Van Essen, Shinjae Yoo, Alex Aiken, David Bernholdt, Suren Byna, Kirk Cameron, Frank Cappello, Barbara Chapman, Andrew Chien, Mary Hall, Rebecca Hartman-Baker, Zhiling Lan, Michael Lang, John Leidel, Sherry Li, Robert Lucas, John Mellor-Crummey, Paul Peltz Jr., Thomas Peterka, Michelle Strout, and Jeremiah Wilke. 2018.  Extreme Heterogeneity 2018 - Productive Computational Science in the Era of Extreme Heterogeneity: Report for DOE ASCR Workshop on Extreme Heterogeneity.  (12 2018).  https://doi.org/10.2172/1473756
  * Wu et al. (2016) Wei Wu, George Bosilca, Rolf Vandevaart, Sylvain Jeaugey, and Jack Dongarra. 2016.  GPU-Aware Non-contiguous Data Movement In Open MPI. In _Proceedings of the 25th ACM International Symposium on High-Performance Parallel and Distributed Computing_. 231–242.

##  Appendix A Artifact Description

The GrayScott.jl implementation is available on GitHub:
https://github.com/JuliaORNL/GrayScott.jl.

The repository contains instructions for the code and how to run it locally,
including setting up the JSON configuration files used in this study, which
are also available:
https://github.com/JuliaORNL/GrayScott.jl/blob/main/examples/settings-
files.json.

To run the codes on Frontier, we created a specific branch to avoid issues
related to CUDA.jl and random number generation on AMDGPU.jl until fixes are
available:
https://github.com/JuliaORNL/GrayScott.jl/tree/amdgpu-frontier.

We provide a configuration file to set up the modules environment as of
June/July 2023: https://github.com/JuliaORNL/GrayScott.jl/blob/amdgpu-
frontier/scripts/config_frontier.sh.

We used the latest available release of Julia, version 1.9.2.