Recalibrating probabilistic forecasts of epidemics

Aaron Rumack1*, Ryan J. Tibshirani1,2, Roni Rosenfeld1

1 Machine Learning Department, Carnegie Mellon University, Pittsburgh, PA, USA

2 Department of Statistics & Data Science, Carnegie Mellon University,
Pittsburgh, PA, USA

*Corresponding author<EMAIL_ADDRESS>
## Abstract

Distributional forecasts are important for a wide variety of applications,
including forecasting epidemics. Often, forecasts are miscalibrated, or
unreliable in assigning uncertainty to future events. We present a
recalibration method that can be applied to a black-box forecaster given
retrospective forecasts and observations, as well as an extension to make this
method more effective in recalibrating epidemic forecasts. This method is
guaranteed to improve calibration and log score performance when trained and
measured in-sample. We also prove that the increase in expected log score of a
recalibrated forecaster is equal to the entropy of the PIT distribution. We
apply this recalibration method to the 27 influenza forecasters in the
FluSight Network and show that recalibration reliably improves forecast
accuracy and calibration. This method is effective, robust, and easy to use as
a post-processing tool to improve epidemic forecasts.

## Author summary

Epidemics of infectious disease cause millions of deaths worldwide each year,
and reliable epidemic forecasts can allow public health officials to respond
to mitigate the effects of epidemics. However, because epidemic forecasting is
a relatively new field and there is often insufficient training data, many
epidemic forecasts are not calibrated. Calibration is a desired property of
any forecast, and we provide a post-processing method that recalisrates
forecasts. We demonstrate the effectiveness of this method in improving
accuracy and calibration on a wide variety of influenza forecasters. We also
show a quantitative relationship between calibration and a forecaster’s
expected score. Our recalibration method is a tool that any forecaster can
use, regardless of model choice, to improve forecast accuracy and reliability.
This work provides a bridge between forecasting theory, which rarely deals
with applications in domains that are new or have little data, and some recent
applications of epidemic forecasting, where forecast calibration is rarely
analyzed systematically.

##  1 Introduction

Epidemic forecasting is an important tool to inform the public health response
to outbreaks of infectious diseases. Often, decision makers can take more
effective action with an estimate of the uncertainty in a forecasted target.
For this reason, distributional forecasts are more desirable than point
forecasts. A distributional forecast is a probability distribution over the
target variable and measures the uncertainty in the prediction, as opposed to
a point forecast, which is just a scalar value for each target and has no
measure of uncertainty. A desired property of distributional forecasts is
calibration, or reliability between forecasts and the true distribution of the
variable forecasted (a mathematical definition is given in Section 2).
Calibration is one of three components of a forecaster’s accuracy as measured
by any proper score [1], with better calibration resulting in a better score.
It is therefore important for a forecaster to produce calibrated forecasts.

Previous work has described general forecasting theory and calibration and
evaluated the calibration of certain forecasts [2, 3, 4, 5]. Later work has
gone beyond just describing calibration, presenting post-processing algorithms
to recalibrate forecasts that were previously miscalibrated. Nonparametric
techniques for recalibration of ensemble forecasts include rank histogram
correction [6], Bayesian model averaging [7], linear pooling [8], and
probability anomaly correction [9]. Brocklehurst et al. [10] provide a
nonparametric approach using the empirical CDF, which can recalibrate any
forecast of a scalar target. Parametric approaches include logistic regression
[11], extended linear regression [12] and beta-transform linear pooling [8].
Wilks and Hamill [13] compare the performance of different recalibration
techniques for different meteorological targets with different amounts of
training data.

Much of the work in recalibration has been applied to weather forecasting, and
thus many of the techniques are not applicable in other forecasting domains.
The most popular weather forecasting models create a distribution from a
series of point predictions, with each point being the result of a simulation
under varying initial conditions. Many of the existing recalibration methods
are defined only for this type of ensemble forecaster. For example, Bayesian
model averaging assumes that an ensemble forecast is comprised of the same $N$
forecasts in each observation. This method cannot be extended trivially to a
domain where the forecaster itself outputs a distribution. Additionally,
weather forecasts usually have a plethora of training data on which to train
recalibration methods. For example, recalibration has been applied to a set of
weather forecasts generated daily from 1979 to at least 2006, almost 10,000
days [14]. In settings like these, techniques need not be robust to small
amounts of recalibration training data.

To be clear on nomenclature, throughout this paper, we use the term forecast
to refer to the predicted probability distribution of a variable and the term
forecaster to refer to an algorithm that produces a forecast for a variable
given a context. Common examples of forecasters are an algorithm that
forecasts the amount of precipitation two days in advance given current
meteorological information, one that forecasts the price of a certain stock
given the stock’s historical trend, or one that forecasts the statewide
influenza incidence given historical incidence data. We also distinguish
between calibration and recalibration; calibration refers to the property of a
forecaster, and recalibration refers to a method whose goal is to make a
forecaster more calibrated. Specifically, recalibration takes as input a set
of a forecaster’s forecasts and corresponding observations (“training data”),
and outputs a forecaster which should be more calibrated on a different set of
forecasts and observations (“test data”).

In what follows, we present a generalized approach to forecast recalibration
and show its performance when applied to forecasters in the FluSight Network.
We demonstrate that across the diverse set of FluSight forecasters,
recalibration consistently improves not just calibration but accuracy as well.

##  2 Methods

Consider the following setup. At each $i=1,2,3\ldots$, a forecaster $M$
outputs a density forecast $f_{i}$ given observations $x_{i}$ for a
continuously distributed scalar random variable $y_{i}$ whose true
distribution is $h_{i}$. We assume that the corresponding cumulative
distribution functions (CDFs) $F_{i}$ and $H_{i}$ are continuous and strictly
increasing (for simplicity). The forecaster $M$ is evaluated according to a
proper scoring rule, such as the Brier score [15] or the logarithmic score
[16].

The goal of a forecaster is to produce ideal forecasts, i.e., to forecast
$f_{i}=h_{i}$, the true distribution of $y_{i}$, for each $i$, though this is
usually unattainable. We can inspect how close a forecaster is to being ideal
with the distribution of the probability integral transform (PIT) values [17].
For each forecast $f_{i}$ and observed value $y_{i}$, the PIT is defined as

$\mathrm{PIT}(f_{i},y_{i})=F_{i}(y_{i}),$

where $F_{i}$ is the CDF of $f_{i}$. A necessary (but not sufficient)
condition for a forecaster to be ideal is probabilistic calibration [3]:

$\frac{1}{N}\sum_{i=1}^{N}H_{i}\circ F_{i}^{-1}(p)\to p\;\;\text{as
$N\to\infty$},\;\;\text{for all $p\in(0,1)$}.$

(Here and throughout we interpret convergence in the almost sure sense.) An
example of a probabilistically calibrated forecaster that is not ideal is the
so-called climatological forecaster, which for each $i$ outputs the marginal
distribution of $y_{i}$ over $i=1,2,3,\ldots$. (To make this concrete, suppose
that each $y_{i}$ is distributed as $\mathcal{N}(\mu_{i},1)$, a normal
distribution with mean $\mu_{i}$ and variance 1, and each $\mu_{i}$ itself
follows $\mathcal{N}(0,1)$, then the climatological forecaster simply outputs
$\mathcal{N}(0,2)$ for each $i$.)

Note that the PIT distribution of a probabilistically calibrated forecaster is
close to uniform in large samples. The expected CDF of the PIT distribution is

$G(p)=\mathbb{E}[\mathbb{P}[F_{i}(y_{i})\leq
p]]=\mathbb{E}[\mathbb{P}[y_{i}\leq F^{-1}_{i}(p)]]=\mathbb{E}[H_{i}\circ
F_{i}^{-1}(p)],$

where here $\mathbb{E}$ denotes the sample average operator over
$i=1,\ldots,N$. This expression converges to $p$ as $N\to\infty$ when the
forecaster is probabilistically calibrated. Thus an examination of the
distribution of PIT values—looking for potential deviations from
uniformity—serves as a good diagnostic tool to assess probabilistic
calibration. Many use a PIT histogram to examine the PIT distribution because
it is easy to read and understand [3]. For example, if the PIT distribution is
bell-shaped, then the forecaster does not put enough weight in the middle of
its distribution and is underconfident. In general, we can compare the PIT
density to the horizontal line at 1, which corresponds to the uniform density.
The greater the deviation from this line (which can be quantified via
Kullback-Leibler divergence from the uniform distribution to the PIT
distribution, or equivalently, negrative netropy of the PIT distribution), the
greater the miscalibration; see Fig 2 for examples.

Our recalibration method uses $G$ as a CDF-CDF transform. The recalibrated
forecaster, denoted $M^{*}$, is defined by a recalibrated forecast CDF of
$F_{i}^{*}(y)=G(F_{i}(y))$, for each $i$. By the chain rule, the recalibrated
forecast density is $f_{i}^{*}(y)=g(F_{i}(y))\cdot f_{i}(y)$, for each $i$.
Thus the recalibrated forecast $f_{i}^{*}$ is the original forecast $f_{i}$
weighted by the PIT density $g$. An illustration of this method is provided in
Fig 2. In practice, of course, we do not have access to the true distributions
$H_{i}$, so we need to estimate $G$. The ultimate estimate of $G$ that we
propose in this paper will be an ensemble (weighted linear combination) of
three estimates: a nonparametric method, a parametric method, and a null
method. First, we will motivate calibration as a tool to increase forecast
accuracy, and then, we explain the individual estimation methods.

Fig 1: Densities of PIT distributions for five sample forecasters, when the
true distribution is a standard normal.

Fig 2: An illustration of recalibration. The original, underconfident forecast
density is $f(y)=\mathcal{N}(0,2)$ while the true density is
$h(y)=\mathcal{N}(0,1)$. By calculating the PIT density $g$ and producing a
recalibrated forecast as the product $g(F(y))\cdot f(y)$, we recover the true
$h(y)$.

###  2.1 Calibration and log score

In order to quantify how well a forecaster is calibrated, we calculate the
entropy of the distribution of PIT values. As above, $G$ is the CDF of the PIT
distribution of $M$. The entropy of the PIT density $g$ is defined as

$H(g)=-\int_{p=0}^{1}g(p)\log g(p)\,dp.$

If $M$ is probabilistically calibrated, then (asymptotically, as $N\to\infty$)
the PIT values are uniform and the entropy is zero because $g(p)$ is $1$
everywhere. When the PIT values are not uniform, the entropy is negative.

Entropy is also useful because it provides an understanding of how
miscalibration penalizes the expected log score, as shown below. First observe
that

$g(p)=\frac{d}{dp}G(p)=\frac{d}{dp}\mathbb{E}[H_{i}\circ
F^{-1}_{i}(p)]=\mathbb{E}\bigg{[}\frac{h_{i}(F^{-1}_{i}(p))}{f_{i}(F^{-1}_{i}(p))}\bigg{]},$

where the last step assumes the smoothness and integrability conditions on
$h_{i},f_{i}$ needed to exchange expectation and differentiation (the Leibniz
rule). Next observe that

$\displaystyle\mathbb{E}[\log f^{*}_{i}(y_{i})]-\mathbb{E}[\log f_{i}(y_{i})]$
$\displaystyle=\mathbb{E}[\log g(F_{i}(y_{i}))]$
$\displaystyle=\mathbb{E}\Bigg{[}\int_{-\infty}^{\infty}\log
g(F_{i}(y_{i}))h_{i}(y_{i})\,dy_{i}\Bigg{]}$
$\displaystyle=\mathbb{E}\Bigg{[}\int_{0}^{1}\log
g(F_{i}(F_{i}^{-1}(p)))\frac{h_{i}(F_{i}^{-1}(p))}{f_{i}(F_{i}^{-1}(p))}\,dp\Bigg{]}$
$\displaystyle=\int_{0}^{1}\mathbb{E}\Bigg{[}\log
g(p)\frac{h_{i}(F_{i}^{-1}(p))}{f_{i}(F_{i}^{-1}(p))}\Bigg{]}\,dp$
$\displaystyle=\int_{p=0}^{1}g(p)\log g(p)=-H(g),$ (1)

where the third line is obtained by a variable substitution, and fourth by
applying the Leibniz rule again assuming the needed regularity conditions.

For any forecaster, if the PIT distribution is the same for the training data
and the test data, then the improvement of the recalibrated forecast’s log
score can be estimated by estimating the negative entropy of $g$ (note that
the entropy of any distribution on $[0,1]$ is nonpositive). We can explain
this intuitively as well: the more negative $H(g)$ is, the more it indicates
that there is information lying in the structure of $g$ that can be extracted
to improve forecasts.

###  2.2 Nonparametric correction

Given an observed training set of PIT values for a forecaster, $F_{i}(y_{i})$,
$i=1,\ldots,N$, the empirical PIT CDF is

$\hat{G}(x)=\frac{1}{N}\sum_{i=1}^{N}\mathbb{I}[F_{i}(y_{i})\leq x].$

As $\hat{G}$ is discrete, it does not admit a well-defined density, and hence
to use this for recalibration we can first smooth $\hat{G}$ using a monotone
cubic spline interpolant, and then it will have a bonafide density $\hat{g}$,
which is itself smooth (twice continuously differentiable, to be precise).
Using this for recalibration produces $f_{i}^{*}(y)=\hat{g}_{i}(F_{i}(y))\cdot
f_{i}(y)$.

In practice, with a large amount of training data, recalibration using the
empirical CDF as described above can be effective. However, with little
training data, or a lot of diversity within the training data among the
distributions of $y_{i}$, it can be ineffective for assuring calibration on
the test set. This is in line with the practical difficulties of using
nonparametric, distribution-free methods in general.

###  2.3 Parametric correction

Gneiting and Ranjan [8] present a recalibration method originally motivated by
redistributing weights on the components of an ensemble forecast, but their
method can applied generally to recalibrate any black box forecaster. Given an
observed training set of PIT values, $F_{i}(y_{i})$, $i=1,\ldots,N$, we fit a
beta density $\hat{g}$ via maximum likelihood estimation. This in fact
corresponds to the beta transform that maximizes the log score of the
recalibrated forecaster on the training data [8].

This parametric model is more resilient to minimal training data, and a beta
distribution is usually an effective estimate of the PIT distribution: because
a beta density can be either convex or concave, it is flexible enough to fit
the PIT distribution of overconfident and underconfident forecasters; and
because the mean can be in the interval $(0,1)$, it can fit biased forecasters
as well. However, problematic behaviors arise at the tails. Except in
exceptional cases (one or both of its two shape parameters is exactly 1), the
beta density is 0 or $\infty$ at the endpoints of its support, which can cause
problmes for recalibration (there can be a big gap between the true PIT
density and $\hat{g}$ in the tails).

###  2.4 Null correction

The final component of the recalibration ensemble is a null correction, in
which there is no recalibration at all, i.e., we simply set
$f_{i}^{*}(y)=f_{i}(y)$. This prevents overfitting and decreases variance of
the overall ensemble correction, to be described next.

###  2.5 Recalibration ensemble

The final recalibration system uses the three components described previously
and weights them in an ensemble. The ensemble weights are calculated to
maximize the overall log score. Letting $f^{*}_{ij}$ denote the forecast
density for sample $i$ and component $j$, the weights ensemble $w$ are defined
by solving the optimization problem:

$\mathop{\mathrm{minmize}}_{w}\;\;\frac{1}{N}\sum_{i=1}^{N}\log\Bigg{(}\sum_{j=1}^{p}w_{j}f^{*}_{ij}(y_{i})\Bigg{)}\;\;\mathop{\mathrm{subject\,\,to}}\;\;w\geq
0,\;\sum_{j=1}^{p}w_{j}=1,$ (2)

where $p$ is the number of ensemble components (for us, $p=3$) and the
constraint $w\geq 0$ is to be interpreted componentwise.

###  2.6 Recalibration under seasonality

Epidemic forecasting presents a new challenge for recalibration. The
methodology discussed above assumes that the previous behavior of a forecaster
is indicative of future behavior, or more concretely, that the PIT
distribution on the training set will be similar to that on the test set.
However, this is not necessarily the case in epidemic forecasting, due to the
fact that a forecaster’s behavior generally changes across the different
phases of an epidemic. For example, some forecasters are too conservative and
do not predict enough of a change in disease incidence from one week to the
next. For such a forecaster, the PIT values are usually too high between a
season’s onset and peak, because incidence increases more quickly than
forecasted. Conversely, after the season peaks, the PIT values are too low,
because incidence decreases more quickly than forecasted.

In order to account for such nonstationarity in the PIT distribution, and at
the same time take advantage of the seasonal nature of epidemic forecasting,
we can recalibrate a forecast made in a particular week $i$ of a given
epidemic season by forming and using a special training set based on forecasts
made at nearby weeks in different seasons. For example, a forecast made in
week 6 can be recalibrated based on forecasts in other seasons made in weeks
in between 3 and 9\. This is what we do in our experiments in this paper, with
more details given in the next section.

##  3 Results

We apply this ensemble recalibration method to data from influenza forecasting
in the US. In an effort to better prepare for seasonal influenza, the US CDC
has organized a seasonal influenza forecasting challenge every year since
2013, called the FluSight Challenge [18]. In 2017, a group of forecasters
formed the FluSight Network [19] and began submitting an ensemble forecast of
27 component forecasters. As part of this collaboration, each of these
forecasters produced and stored retrospective forecasts spanning 9 seasons,
from 2010-11 to 2018-19. These forecasters include mechanistic and non-
mechanistic forecasters, as well as baseline forecasters. They are diverse in
behavior, accuracy, and calibration, and therefore provide an interesting
challenge for our recalibration method, which treats the forecaster as a black
box.

First, we summarize the retrospective forecasts in the FluSight data set. Each
week, a forecast is produced for seven forecasting targets, all of which are
based on weighted ILI (wILI), a population-weighted average of the percentage
of outpatient visits with influenza-like illness derived from reports to the
CDC from a network of healthcare providers called ILINet [20]. The forecasting
targets are:

  * •

season onset (the first week where wILI is above a predefined baseline for
three consecutive weeks);

  * •

season peak week (week of maximum wILI);

  * •

season peak percentage (maximum wILI value);

  * •

the wILI value at 1, 2, 3, and 4 weeks ahead of the current week.

The first three targets are referred to as seasonal targets and the last four
targets are referred to as short-term targets. Each forecast is submitted as a
binned probability distribution, in intervals of weeks for season onset and
season peak week, and intervals of 0.1% wILI for season peak percentage and
all short-term targets. Forecasts are produced for each of the 10 HHS Regions
as well as the US as a whole, for a total of 9 seasons, from 2010-11 to
2018-19. Thus to be clear, the forecasts in this FluSight data set are indexed
by forecaster, target, season, forecast week, and location.

Next, we describe the training setup we use for recalibrating the forecasts in
this data set, which is a kind of nested leave-one-season-out cross-
validation. This is laid out in the steps below for a given forecaster and
forecasting target, and a particular season $s$.

  1. 1.

Create recalibrated forecasts for all seasons $r\neq s$, using each of the
three methods: nonparametric, parametric, and null. For a forecast in season
$r$ at week $i$ and at location $\ell$, we build a training set using PIT
values from all seasons other than $r$ and $s$, all forecast weeks in
$[i-3,i+3]$ (within three weeks of $i$), and all locations.

  2. 2.

Optimize the ensemble weights $w$ by solving (2) using the recalibrated
forecasts from Step 1.

  3. 3.

Create recalibrated forecasts for season $i$, again using each of the three
methods: nonparametric, parametric, and null. This is just as in Step 1,
except we use one more season in the training set. Explicitly, for a forecast
in season $s$ at week $i$ and at location $\ell$, we build a training set
using PIT values from all seasons other than $s$, all forecast weeks in
$[i-3,i+3]$ (within three weeks of $i$), and all locations.

  4. 4.

Create ensemble recalibrated forecasts from season $i$, using the
recalibration components from Step 3 and the weights from Step 2.

In what follows, we present and discuss the results. The code and data used to
produce all of these results is publicly available online [21].

###  3.1 Effect of varying window size

The training procedure just presented assumes a window of $k=3$ weeks on
either side of a given week $i$ in order to build the set of PIT values used
for recalibration (using forecast data from other seasons). However, we could
consider varying $k$, which would navigate something like a bias-variance
tradeoff. We would expect the optimal window $k$ to be larger for the
nonparametric recalibration method versus the parametric one. It turns out
that $k=3$ is typically a reasonable choice for both, as displayed in Fig 3,
an example using data from the forecaster scoring highest on the short-term
targets.

Fig 3: Mean log score for the most accurate forecaster of the short-term
targets, for the different recalibration methods. A window size of $k$
corresponds to training recalibration on forecasts within $k$ weeks of the
given forecast week, inclusive. Log score is averaged over 9 seasons, 11
locations, and 29 weeks (higher log score is better). The window size hardly
affects performance of the parametric recalibration model. However, the
smallest window sizes hurt the nonparametric model.

###  3.2 Forecast accuracy and calibration

For the short-term targets, the ensemble recalibration method improves the
mean log score for almost all forecasters. Both the nonparametric and
parametric recalibration methods significantly improve the mean log score, and
the ensemble improves it even further. For the seasonal targets, some
component recalibration methods do not improve accuracy, although the ensemble
method does improve accuracy, averaged over all forecasters. However, the
ensemble improves accuracy for seasonal targets in only about three-quarters
of forecasters. See Fig 5 and Fig 5.

Fig 4: Improvement in mean log score, for the different recalibration methods.
Log score is averaged over all 27 forecasters in the FluSight, 9 seasons, 11
locations, and 29 weeks (higher log score is better). The ensemble
recalibration method improves accuracy for every target.

Fig 5: Proportion of forecasters for which recalibration improves mean log
score (left) and entropy of the PIT values (right). The ensemble method
improves accuracy for the short-term targets for all forecasters, and most
forecasters for the seasonal targets. It also improves calibration (as
measured by entropy) for most forecasters and most targets. The ensemble
method outperforms both the nonparametric and parametric methods.

Furthermore, in Fig 5 we see that the recalibrated forecaster’s improvement in
calibration (as measured by entropy) is quite close to the improvement in mean
log score. This confirms our expectations from (1). Fig 7 gives a more direct
comparison of improvements in accuracy versus calibration, i.e., in mean log
score versus entropy, for the short-term forecasts. (Note that we estimate the
entropy of the distribution of PIT values using a simple histogram estimator
with 100 equal bins along the interval $[0,1]$.) We see a clear linear trend,
with slope approximately 1, again confirming (1).

Finally, in Fig 7, we show that our ensemble recalibration method increases
the entropy of the PIT distribution to nearly zero for nearly every
forecaster. The two exceptions, the line segments towards the bottom of Fig 7,
correspond to particularly poor forecasters (so poor that are outperformed by
a baseline forecaster that outputs a uniform distribution).

Fig 6: Improvement in mean log score versus improvement in entropy for each of
the 27 FluSight forecasters and short-term targets. There is a clear linear
trend (with slope approximately 1) between the improvement in calibration and
the improvement in accuracy.

Fig 7: Entropy and mean log score before and after recalibration, for each of
the 27 FluSight forecasters and short-term targets. The tail of arrow
represents a quantity before recalibration, and the head after recalibration.
The dotted lines show the central 90% interval of the entropy of a comparably-
sized sample of standard uniform random variables for comparison. For all but
two forecasters (the eight line bottom-most line segment), the ensemble
recalibration method achieves almost perfect calibration as evidenced by a
near-zero PIT entropy, and this is accompanied by significant improvements in
accuracy.

###  3.3 Recalibrating the FluSight ensemble

As we just saw, recalibration improves the performance of the individual
forecasters in the FluSight Network. A natural follow up is therefore to
investigate whether it can improve the performance of the FluSight ensemble, a
forecaster that combines 27 component forecasters (the individual FluSight
forecasters), whose construction is described in [19].

As both recalibration and ensembling are post-processing methods (i.e., that
can be applied in post-processing of forecast data), we are left with two
options to explore. We can recalibrate the component forecasters and then
ensemble (C-E), or ensemble the components and then recalibrate (E-C). In the
C-E model, we train ensemble weights in a leave-one-season-out format, on the
recalibrated component forecasts. In the E-C model, we train ensemble weights
in a leave-one-season-out format on the original component forecasts, and then
recalibrate the ensemble forecasts.

Fig 8 reveals that E-C model performs better than the C-E model. This is in
line with established forecasting theory, which states that linear ensembles
(which take a linear combination of component forecasters, such as the
FluSight ensemble approach) themselves are generally miscalibrated, even when
the individual component forecasters are themselves calibrated [5, 22, 8].

Fig 8: Mean log score for the two different approaches to recalibrating the
FluSight ensemble forecaster, with C-E and E-C reflecting the order of
recalibration and ensembling. Both the C-E and E-C models outperform the
original ensemble (with no recalibration), but ensembling followed by
recalibration performs best. By viewing forecast performance as a function of
time, recalibration increases performance as much as roughly two days’ time
would.

##  4 Discussion

Even in a domain as complex as epidemic forecasting, relatively simple
recalibration methods such as those described in this paper can significantly
improve both calibration and accuracy. A forecaster’s performance for any
proper score can be decomposed into three components: the inherent uncertainty
of the target itself, the resolution of the forecaster (concentration of the
forecasts), and the reliability of the forecaster to the target (calibration)
[1]. In epidemic forecasting, without seasonality-aware recalibration training
(such as that proposed and implemented in this paper), recalibration will not
affect the resolution term, which is left to the individual forecasters, but
it will improve the reliability term. However, using seasonality-aware
recalibration, it can also improve the resolution term.

Over 9 seasons of forecast data from 27 forecasters in the FluSight Challenge,
we found that recalibration was especially helpful for the short-term targets
(1-4 week ahead forecasts). With the exception of two very similar forecasters
that have poor performance, the ensemble recalibration method was able to
reduce the entropy of the PIT distribution to nearly zero (not or barely
statistically significantly different than a uniform distribution). The
recalibrated forecasts are therefore more accurate and more reliable. This is
true across a diverse set of forecasters, including mechanistic, statistical,
baseline, and ensemble models; indeed, as our recalibration method treats the
forecaster as a black box, it can be applied to any forecaster, given access
to suitable training data (retrospective historical forecasts).

The performance of recalibration with respect to the seasonal targets (onset,
peak week, and peak percentage) was less conclusive. Although the mean log
score averaged over all of the forecasters was improved, recalibration only
improved the performance of about three-quarters of the forecasters. Seasonal
targets are inherently more difficult to recalibrate because at the end of the
season, the true value has almost certainly been observed, and the forecasts
are highly confident. For these forecasts, the correct bin has a mass of
almost 1, and the observed PIT value then is approximately 0.5. At the end of
the season, the PIT distribution is very concentrated at 0.5, which indicates
underconfidence and poor calibration. If these PIT values of 0.5 are used to
train forecasts for recalibration earlier in the season, before the target is
observed, then recalibration incorrectly makes the forecast more confident.
Because one is unsure whether the season peak has occurred or not for several
weeks after the peak occurs, recalibration training is a nontrivial task. In
general, more work is required to reliably improve accuracy and calibration
for seasonal targets, which is a topic for future work.

## Acknowledgements

AR was supported by a fellowship from the Center for Machine Learning at
Carnegie Mellon University and a gift from the McCune Foundation. RT and RR
were support by Centers for Disease Control and Prevention grant U01IP001121.

## References

  * 1. Bröcker J.  Reliability, sufficiency, and the decomposition of proper scores.  Quarterly Journal of the Royal Meteorological Society. 2009;135(643):1512–1519.
  * 2. Dawid AP.  Calibration-based empirical probability.  Annals of Statistics. 1985;13(4):1251–1274.
  * 3. Gneiting T, Balabdaoui F, Raftery AE.  Probabilistic forecasts, calibration and sharpness.  Journal of the Royal Statistical Society, Series B. 2007;69:243–268.
  * 4. Gneiting T, Raftery AE.  Strictly proper scoring rules, prediction, and estimation.  Journal of the American Statistical Association. 2007;102(477):359–378.
  * 5. Hora SC.  Probability judgements for continuous quantities.  Management Science. 2004;50(5):597–604.
  * 6. Hamill TM, Colucci SJ.  Verification of Eta-RSM short-range ensemble forecasts.  Monthly Weather Review. 1997;125:1312–1327.
  * 7. Raftery AE, Gneiting T, Balabdaoui F, Polakowski M.  Using Bayesian model averaging to calibrate forecast ensembles.  Monthly Weather Review. 2005;133:1155–1174.
  * 8. Gneiting T, Ranjan R.  Combining predictive distributions.  Electronic Journal of Statistics. 2013;7:1747–1782.
  * 9. van den Dool H, Becker E, Chen LC, Zhang Q.  The probability anomaly correlation and calibration of probabilistic forecasts.  Weather and Forecasting. 2017;32:199–206.
  * 10. Brocklehurst S, Chan PY, Littlewood B, Snell J.  Recalibrating software reliability models.  IEEE Transactions on Software Engineering. 1990;16(4):458–470.
  * 11. Hamill TM, Whitaker JS, Wei X.  Ensemble reforecasting: Improving medium-range forecast skill using retrospective forecasts.  Monthly Weather Review. 2004;132:1434–1447.
  * 12. Gneiting T, Raftery AE, Westveld III AH, Goldman T.  Calibrated probabilistic forecasting using ensemble model output statistics and minimum CRPS estimation.  Monthly Weather Review. 2005;133:1098–1118.
  * 13. Wilks DS, Hamill TM.  comparison of ensemble-MOS methods using GFS reforecasts.  Monthly Weather Review. 2007;135:2379–2390.
  * 14. Hamill TM, Whitaker JS, Mullen SL.  Reforecasts: An important dataset for improving weather predictions.  Bulletin of the American Meteorological Society. 2006;87:33–46.
  * 15. Brier GW.  Verification of forecasts expressed in terms of probability.  Monthly Weather Review. 1950;78(1):1–3.
  * 16. Good IJ.  Rational decisions.  Journal of the Royal Statistical Society. 1952;14:107–114.
  * 17. Dawid AP.  Statistical theory: The prequential approach (with discussion).  Journal of the Royal Statistical Society Series A. 1984;147:278–292.
  * 18. Centers for Disease Control and Prevention. FluSight: Flu forecasting;.  https://www.cdc.gov/flu/weekly/flusight/index.html.
  * 19. Reich NG, McGowan CJ, Yamana TK, Tushar A, Ray EL, Osthus D, et al.  Accuracy of real-time multi-model ensemble forecasts for seasonal influenza in the U.S.  PLOS Computational Biology. 2019;15(11).
  * 20. Centers for Disease Control and Prevention. National, regional, and state level outpatient illness and viral surveillance;.  https://gis.cdc.gov/grasp/fluview/fluportaldashboard.html.
  * 21. Rumack A, Brooks LC, Hyun S, Tibshirani RJ, Rosenfeld R. epiforecast-R/recalibration repository;.  https://github.com/cmu-delphi/epiforecast-R/tree/recalibration.
  * 22. Ranjan R, Gneiting T.  Combining probability forecasts.  Journal of the Royal Statistical Society, Series B. 2010;72(1):71–91.


# Calculation of Heavy Meson Potential Coefficients by Solving Lippman-
Schwinger Equation

P Sadeghi Alavijeh , M Monemzadeh and N Tazimi

Department of Physics, University of Kashan, Kashan, Iran
ph.parva@yahoo.commonem<EMAIL_ADDRESS>
###### Abstract

In the present work, we study meson systems consisiting of quark-antiquark. We
solve Lippman-Schwinger equation numerically for heavy meson systems . We
attempt to find a non-relativistic potential model through which we can solve
the quark-antiquark bound state problem. The co-efficients so obtained are in
agreement with Martin potential coefficients. Via this method we also
determine the strong coupling constant for the mesons $t\bar{t}$, $t\bar{c}$
and $t\bar{b}$ which is a coefficient of Cornel potential.

Keywords: binding energy, lippman-Schwinger, Martin potential, strong coupling
constant

PACS Numbers: 12.39.Jh, 12.39.Pn, 14.40.Pq, 14.40.Rt, 21.10.Dr

##  1 Introduction

The heavy quark potential is an important quantity associated with
confinement. This quantity can be almost accurately calculated via lattice
simulations [1]. This potential is speciﬁed through a non-relativistic
effective theory, yet it also deserves serious investigation in perturbative
QCD. Lattice QCD simulations are used for analyzing at long distances, but for
short distances, perturbation theory works better.

There exists a variety of experimental data on light and heavy mesons, and it
is still increasing [2]. When the heavy quark mass ($m_{Q}$) remarkably
exceeds $\Lambda_{QCD}$ (the QCD scale), the running coupling
$\alpha_{s}(m_{Q})$ is small. This means that at this scale of the order of
the Compton wavelengh $\lambda\sim 1/m_{Q}$, perturbative QCD can serve to
describe hadrons. In the non-relativistic limit, it has been shown that the
interaction between the two $\bar{Q}$ and $Q$ states can be described by a
local potential $V(r)$, where $r$ is the relative coordinate between $Q$ and
$\bar{Q}$ (spin is ignored for the moment) [1].

We exploit Lippman-Schwinger equation to analyze the quark- antiquark bound
state. This equation has recently been used in some studies of non-
relativistic bound state [3, 4, 5]. We have investigated some potential models
for heavy mesons. We have also identified the stability intervals of some
suitable potentials for these systems and obtained the mass spectrum of these
mesons. In our recent work [6], we have studied the tetraquark and solved the
two-body problem for tetraquark systems of diquark-antidiquark.

Having studied some types of appropriate potentials such as Martin [7] and
Cornel [8, 9], we found Martin potential as the best one for heavy mesons.
This is because of the potential’s larger stability interval [3]. In the
present work, we calculate the coefficients of Cornell and Martin potentials
(which are suitable for describing heavy mesons) through numerical solution of
Lippman-Schwinger equation and obtain the physical eigenvalue. It must be
notified that we have used input parameters diﬀerent from those of ref [3]
(quark mass, binding energy, and potential). In ref [3], we utilized a revised
potential, but we have made use of Cornell potential in the present paper.

It is worth mentioning that in our previous work [3], we detected appropriate
potentials through this method and obtained $\lambda=1$ within an extensive
interval of $r$; therefore, we have been certain about the appropriateness of
Cornell and Martin potentials for heavy mesons in our method. That is, we took
the binding energy as already identified and solved Schrodinger equation to
find two potential coefficients for each meson.

The procedure of the study is given in sec 1. We explain how to find Martin
potential coefficients in sec 2. In sec 3, we have also used Coulomb potential
confinement to predict the strong coupling constant for $t\bar{t}$, $t\bar{c}$
and $t\bar{b}$. The results and discussion appear in sec. 4 and finally, sec.
5 gives a concise conclusion of the study.

##  2 Method

In this section, we obtain the coefficients of a suitable potential for heavy
meson systems by using heavy mesons’ binding energies and by solving Lippman-
Schwinger equation numerically.

The two-body bound state of Schrodinger equation with potential $V$ is:

$\mid\psi_{b}>=G_{0}V\mid\psi_{b}>$ (1)

where $G_{0}$ is the propagator of a free particle. In configuration space, it
is shown as:

$\psi_{b}(r)=-m\frac{1}{4\pi}{\int_{0}^{\infty}}{dr^{\prime}}{r^{\prime}}^{2}{\int_{-1}^{1}}dx^{\prime}{\int_{0}}^{2\pi}d\phi^{\prime}\frac{exp(-\sqrt{m|E_{b}|}|{r-r^{\prime}}|)}{|{r-r^{\prime}}|}V(r^{\prime}){\psi_{b}(r^{\prime})}$
(2)

${\psi_{b}(r)}=\int_{0}^{\infty}dr^{{}^{\prime}}\int_{-1}^{1}dx^{{}^{\prime}}M(r,r^{{}^{\prime}},x^{{}^{\prime}}){\psi_{b}(r^{{}^{\prime}})}$
(3)

where:

$M(r,r^{{}^{\prime}},x^{{}^{\prime}})=-\frac{m}{2}\frac{exp((-\sqrt{m|E_{b}|)}\sqrt{r^{2o}+\acute{r}^{2}-2r\acute{r}\acute{x}})}{\sqrt{r^{2}+\acute{r}^{2}-2r\acute{r}\acute{x}}}{}\\\
{V(\acute{r}^{2})}$ (4)

$E_{b}$ stands for the binding energy of two-body bound system. The eigenvalue
form of equation (3) is:

$K(E_{b})|\psi_{b}>=\lambda(E_{b})|\psi_{b}>$ (5)

To solve this equation, we use a Fortran code. It is used to calculate the
eigenvalues for a real non-symmetric square matrix. The magnitude of
eigenvalue $\lambda$ depends on energy. The energy that results in $\lambda=1$
is the binding energy. $\lambda=1$ is the highest positive eigenvalues. The
eigenvalue equation can be solved by direct iteration method [10]. To
discretize it, the integral Gauss-legendre method [11] is used. It is worth
mentioning that to get better results out of the calculations, we consider
Gaussian quadrature points for $r$, $r^{{}^{\prime}}$, and $x^{{}^{\prime}}$
as 100 (The more the points, the more accurate the results, although this
lowers the running speed of the program).
The required input for the program includes system mass, $E_{b}$, potential
coefficients, and r-cutoff. R-cutoff, which is of the same order of magnitude
as the meson radius, is supposed to be the point where the potential tends to
zero. $E_{b}$ is obtained via solving an eigenvalue equation by using quark-
antiquark interaction (spin-spin interaction in the potential and spin
splitting are ignored). The masses of constituent quarks (MeV) are fitted as:

$m_{c}=1800,m_{b}=5174,m_{t}=174000$

Ref [12] also makes use of the same inputs for Martin potential. Table 1 shows
heavy meson masses and heavy meson binding energies. Binding energy is deﬁned
as the energy used when breaking a meson into its components, i.e.quark and
antiquark, so it is negative 111In order to separate the quark from the
antiquark, huge energy of the order of several MeV is required, i.e. the rest
mass energy of the quarks in free state exceeds their total energy while they
are inside the meson. The energy required to separate the quarks is the
binding energy of the meson.. As we know, meson mass is calculated through
equation (6)

$M(meson)=m_{quark}+m_{antiquark}+E_{b}$ (6)

Table 1: Heavy meson mass and binding energies (MeV)  Spin | $\quad Meson$ | Meson Mass (PDG[13]) | $E_{b}$
---|---|---|---
S=1 | $\quad c\bar{c}(J/\psi)$ | $\quad 3096.9$ | $\quad-503.1$
S=0 | $\quad c\bar{c}(\eta_{c})$ | $\quad 2979.9$ | $\quad-620.1$
S=1 | $\quad b\bar{b}(\Upsilon)$ | $\quad 9460$ | $\quad-888$
S=0 | $\quad b\bar{b}(\eta_{b})$ | $\quad 9398$ | $\quad-850$

$E_{b}$ in table (1) is obtained from equation (6). Based on the information
presented in table 1, we introduce the reduced mass of the meson’s constituent
quarks and also the energies into the program [3]. First, we place the
potential in the program:

$V(r)=(ra_{n}/\hbar c)^{z}$ (7)

where $a_{n}=1000MeV$ is an input .We have introduced the coefficient $a_{n}$
to correct the dimension so that the dimensionof each term of the potential is
given in MeV. The parameter $\hbar c$ is included in the equation for
dimension coordination. $V$ is expressed in terms of $MeV$ in the program and
$r$ in terms of $fm$. In this part, we assigned $z$ an arbitrary quantity. The
purpose was to obtain $\lambda=1$, so $z$ was assigned lower or higher values.
We carried out this process of iteration until $\lambda=1$ was obtained as the
program output. There is a direct relationship between the coefficient $z$ and
the program output, i.e. the $\lambda$ spectrum. The higher the coefficient,
the higher the values obtained in the $\lambda$ spectrum. The size and
direction of each iteration depends on what $\lambda$ spectrum is obtained
following each iteration. If the values in the spectrum are much
larger/smaller than 1, the next assigned $z$ should naturally be much
larger/smaller than the previous $z$. But when we obtain a $\lambda$ spectrum
very near to 1, the next assigned $z$ will be very near to the previous one.
Of course, the rate of accuracy set in this stage of the study is 0.01, then
$0.99\leq\lambda\leq 1.01$ is acceptable. Thus, the quantities that yielded
eigenvalues either greater than 1.01 or smaller than 0.99 were ignored. The
results of running the program show that at $z$ = 0.1, the desired output is
obtained.

In the next stage, we introduce the coefficient $k$ and take the potential as
$k(r^{0.1})$. We then follow the above procedure, i.e. we obtain $\lambda=1$
by assigning different values to $k$. The rate of accuracy; however, is 0.001.
Therefore, any $0.999\leq\lambda\leq 1.001$ is acceptable. Finally, we
introduce the constant $d$ and increase the accuracy up to 0.0001 . In fact,
we don’t introduce the potential at once but rather stage by stage. Initially,
the coefficient $z$ is introduced and the predetermined acceptable error is
0.01. Later as we introduce the other coefficients, we approach the full form
of the potential and the tolerance is enhanced to 0.01 for $k$ and finally to
0.0001 for $d$.

$V(r)=k(a_{n}r)^{z}+d$ (8)

We are using tables 1 and 2 as input. Table 2 illustrates the results. The
cells marked with an asterisk $(*)$ contain results that comply with the
coefficients of Martin potential. In this potential, $r$ is defined in terms
of $GeV^{-1}$. The cells containing a dash $(-)$ indicate potential
coefficients which didn’t enjoy the desired accuracy. Through the procedure
already described, we have managed to obtain values of $z$, $k$, and $d$ which
are signiﬁcant for different potentials. Our results are spin-dependent
because the potential models we are considering are spin-dependent.

Table 2: Potential exponent and coefficients obtained from the program Meson | z | k(MeV) | d(MeV)
---|---|---|---
| 0.0931 | 6886 | -
cc̄(S=1) |  | 6890 | -
| 0.0972* | 6899* | 8071*
|  | 6910 | -
| 0.0981* | 6990 | -
cc̄(S=0) |  | 6930* | 8980*
|  |  | 8591
| 0.1* | 6849 | 6849
|  |  | 6862
|  | 6862 | -
bb̄(S=0) | 0.0983* | 6883* | 8669 *
|  |  | 9465
Martin potential coefficients [7] | 0.1 | 6898 | 8093

##  3 Calculating Heavy Meson Strong Coupling Constant

The running coupling in QCD is small at high energies, while it increases to
order 1 at low energies. This feature is asymptotic freedom of QCD. Mesons
containing top quark are unstable because of their short lifetime. However,
within this short lifetime, the strong coupling constant can be anticipated
(we are aware that this application may be purely academic because of the non-
existence of the corresponding bound states).

In this section, our anticipated strong coupling constant for $t\bar{t}$,
$t\bar{c}$ and $t\bar{b}$ are presented, thus here we use Cornell potential.
Cornell potential is composed of two terms. One is Coloumb potential used for
interaction in small distances i.e. $V\propto-\alpha_{s}/r$ at $r\rightarrow
0$. The other is a linear potential capable of describing conﬁnement. It is
used for large distances, i.e. $V\propto r$ at large $r$. This potential has
absorbed much attention in hadron physics and is definitely of strong
interaction type. This potential is of the following form [8, 9]:

$V(r)=fr-\frac{4}{3}\frac{\alpha_{s}\hbar c}{r}$ (9)

where $f=150000MeV^{2}$ [15] is the string tension . Here, $\hbar c$ is also
used for dimension coordination. $\alpha_{s}$ is the coefficient to be worked
out. Now, we introduce the energy (table (3)) and potential (eq. 9) into the
program. We assign an arbitrary initial value to $\alpha_{s}$ (we assigned
0.01) and then search through the immediate neighbourhood iteratively with a
step-size of 0.0001 until $\lambda=1\pm 0.01$ is encountered. The results are
presented in table (5). Varying $\alpha_{s}$ leads to diﬀerent $\lambda$
spectra. However, we seek to arrive at results (i.e. spectrum) which are
closest to 1. In other words, accuracy depends on the minimum allowed error of
calculation outcome. It is worth mentioning that strong coupling constant for
$c\bar{c}$ and $b\bar{b}$ mesons can also be calculated by this method.

We are using table 3 as input. Our results are available in tables 4 and 5.
According to [14], coupling constant is inversely correlated with mass, i.e.

$m_{t\overline{c}}\leq m_{t\overline{b}}\leq
m_{t\overline{t}}\Rightarrow\alpha_{t\overline{t}}\leq\alpha_{t\overline{b}}\leq\alpha_{t\overline{c}}$

Our results in table 4 also indicate this inverse correlation. We used this
method to identify the coupling constant of $c\overline{c}$ and
$b\overline{b}$ mesons. The $\alpha_{s}$ we obtained is similar to the
$\alpha_{s}$ in [16]. The results appear in table 5.

As stated in the introduction, we have earlier investigated certain potential
models appropriate for heavy mesons to calculate the binding energy of these
systems. Since this method is potential-independent, the obtained binding
energy applies to other potentials too [17]. In [3], we arrived at the binding
energy of $t\overline{t}$ and $t\overline{c}$ systems through a revised
potential (we have also calculated the binding energy of $t\overline{b}$ by
means of this potential). In this paper, we used that binding energy as the
input of the program and obtained potential coefficients diﬀerent from the
previous work. As two examples, we obtained Martin potential coefficients for
$b\overline{b}$ and $c\overline{c}$ systems and obtained the coupling constant
in Cornell potential for $t\overline{b}$, $t\overline{c}$, $t\overline{t}$,
$b\overline{b}$ and $c\overline{c}$ systems. However, we could have
alternatively calculated Martin potential coefficients for t-quark mesons.

Table 3: Introduced energies for heavy mesons 3 Meson | S | E(MeV)
---|---|---
$t\bar{t}$ | 1 | -11.296
$t\bar{t}$ | 0 | -36.948
$t\bar{b}$ | 0 | -195.0870
$t\bar{c}$ | 0 | -550
$t\bar{b}$ | 1 | -641.7242
Table 4: Strong coupling constants obtained for $t\bar{t}$, $t\bar{b}$, and $t\bar{c}$ mesons Meson | $t\bar{t}(S=0)$ | $t\bar{t}(S=1)$ | $t\bar{b}(S=0)$ | $t\bar{b}(S=1)$ | $t\bar{c}(S=0)$
---|---|---|---|---|---
$\alpha_{s}$ | 0.2340 | 0.1998 | 0.3240 | 0.3060 | 0.4260

##  4 Results and Discussion

We made use of mass, r-cutoff, and binding energy of each meson to solve
Schrodinger equation and diagonalize the kernel. We searched iteratively
through parameter space to obtain coefficients for proposed potential model
until desired eigenvalues were met within acceptable tolerance. We found
exclusive potential coefficients for each meson. The pro- cedure led to the
identiﬁcation of two appropriate potential coefficients for heavy mesons.

It must be notified that the two potentials used in this paper are spin-
independent. The coefficient diﬀerences observed across $S=0$ and $S=1$ states
result from the different binding energies that we used for these two states
(tables 1 and 3). Thus, spin indirectly affects the mass and coefficients of
the potential.

Table 5: Strong coupling constants obtained for $c\bar{c}$ and $b\bar{b}$ mesons Meson |  $\alpha_{s}$ in ref[16] | Calculated $\alpha_{s}$
---|---|---
$c\bar{c}$ | 0.51 | 0.48
$b\bar{b}$ | 0.33 | 0.31

##  5 Conclusions

In this study, we managed to calculate the coefficients of one potential
through solving Lippman-Schwinger equation numerically and obtain the physical
eigenvalue. This po- tential’s coefficients are very close to Martin potential
coefficients. Additionally, since we arrived at acceptable results in the
program, we managed to obtain $\alpha_{s}$ considered as an unidentified
coefficient in the Cornell potential for $t\bar{t}$, $t\bar{b}$, $t\bar{c}$,
$b\bar{b}$ and $c\bar{c}$ mesons. The method practiced, therefore, is a good
one to identify the potential coefficients used for heavy mesons.

##  6 Acknowledgments

We are pleased to thank the University of Kashan for Grant No. 65500.4

## References

  * [1] G S Bali, Phys. Rept. 343 : 1-136 (2001), hep-ph/0001312.
  * [2] K A Olive, et al., Particle Data Group, Chin. Phys. C38 090001 (2014)
  * [3] M Monemzadeh, M Hadizadeh, N Tazimi Int. J. Theor. Phys. Vol 50 No . (2011)
  * [4] N Tazimi, M Monemzadeh and M R Hadizadeh Int. J. Theor. Phys. 52, 2329 (2013)
  * [5] N Tazimi, M Monemzadeh and M R Hadizadeh Int. J. Theor. Phys. 51, 2871 (2012)
  * [6] M Monemzadeh, N Tazimi, P Sadeghi Physics Letters B741 (2015) 124127
  * [7] A Martin Phys. Lett. B100 (1981) 511
  * [8] G S Bali Phys. Rept. 343, 1-136 (2001); N Brambilla, A Pineda, J Soto and A Vairo Rev. Mod. Phys. 77, 1423 (2005)
  * [9] E Eichten, S Godfrey, H Mahlke and J L Rosner Rev. Mod. Phys. 80, 1161 (2008)
  * [10] H W Wyld, W A Benjamin Mathematical Methods for Physics (1976)
  * [11] Shoichiro Nakamura Applied numerical methods software, the Ohio State University, by Prentice Hall. Inc, 1998
  * [12] A Martin Phys.Lett. 100B (1981) 511
  * [13] K Nakamura et al. (Particle Data Group), J. Phys. G37, 075021 (2010)
  * [14] D Ebert, R N Faustov, V O Galkin Mod.Phys.Lett. A20 (2005) 875-880
  * [15] M. N. Sergeenko, Eur. Phys. J. C, Vol. 72 (2012) 2128
  * [16] A M Badalian, A I Veselov and B L G Bakker Phys. Rev. D70, 016007 (2004); Yu A Simonov Phys. Atom. Nucl. 58, 107 (1995)
  * [17] M. Monemzadeh, N. Tazimi, P. Sadeghi, Phys. Lett.B741 , 124 (2015).


# Graph Rewriting for Graph Neural Networks

Adam Machowczyk  University of Leicester, UK  Reiko Heckel  University of
Leicester, UK

###### Abstract

Given graphs as input, Graph Neural Networks (GNNs) support the inference of
nodes, edges, attributes, or graph properties. Graph Rewriting investigates
the rule-based manipulation of graphs to model complex graph transformations.
We propose that, therefore, (i) graph rewriting subsumes GNNs and could serve
as formal model to study and compare them, and (ii) the representation of GNNs
as graph rewrite systems can help to design and analyse GNNs, their
architectures and algorithms. Hence we propose Graph Rewriting Neural Networks
(GReNN) as both novel semantic foundation and engineering discipline for GNNs.
We develop a case study reminiscent of a Message Passing Neural Network
realised as a Groove graph rewriting model and explore its incremental
operation in response to dynamic updates.

##  1 Introduction

Neural Networks (NNs) are among the most successful techniques for machine
learning. Starting out from application data, e.g. in a relational database,
NNs require an encoding of the data into discrete vectors. Graph Neural
Networks (GNNs) operate on graphs instead of tabular data. According to [19] a
GNN is “an optimizable transformation on all attributes of the graph (nodes,
edges, global context) that preserves graph symmetries (permutation
invariances)”. GNNs are applied to [19] images, text, molecules, social
networks, programming code [1] and even mathematical equations [14], in areas
such as antibacterial discovery [22], physics simulation [18], fake news
detection [17], traffic prediction [15] and recommendation systems [7]. They
support supervised, semi-supervised or unsupervised learning (where data
labelling is expensive or impossible) and allow links between input vectors,
resembling more closely the structure of the application data. In a
recommender system where engagements with social media posts are nodes, they
may predict such nodes based on the history of engagement, with attributes
representing the strength of engagement.

Graph rewriting is the rule-based transformation of graphs [12]. Since GNNs
are essentially graph transformations, we ask if they can be modelled by graph
rewriting, and what the benefits of such a representation might be. Drawing
from its rich theory and tool support, we would like to explore the use of
graph rewriting as a computational model for GNNs, to specify and compare GNN
variants, study their expressiveness and complexity, and as a design and
prototyping language for GNNs, evaluating specific GNN architectures and
analysing interactions between GNNs and the software and real-world systems
with which they interact. In modern data-driven applications, data is diverse,
large-scale, evolving, and distributed, demanding incremental and distributed
solutions. This has led to a variety of GNNs [25] and a need for a unifying
theory to describe their expressiveness and cost [4], evolution and
incremental updates [8].

In this short paper, we suggest that typed attributed graph rewriting with
control can address these demands. We will illustrate this claim by means of a
case study of a recommender system for posts in a social network modelled as a
_Graph Rewriting Neural Network (GReNN)_ , a controlled graph rewriting system
based on a cycle of training and inference transformations with the
possibility of dynamic data updates. Training involves the tuning of
parameters such as weights and thresholds to support the inference of graph
features or properties. The example does not follow a particular GNN
architecture, but can be seen as an elaboration of message-passing GNNs over
directed heterogeneous graphs with dynamic changes and incremental updates, a
combination we have not encountered in any existing formalism. In particular,
typed attributed graphs allow for heterogeneity, where nodes of different
types have different attributes and are connected by differently-typed edges,
while the rule-based approach supports both local inference operations of
message-passing GNNs and updates to the graph structure representing changes
in the real-world data.

##  2 A Recommender System GReNN

We consider a social network with users, posts, and engagements, e.g., reads,
likes or retweets, modelled as a Groove graph rewriting system [9]. In the
type graph

users have an _upd_ attribute that determines if training is required on
engagements with posts by this user. Posts have _weight_ s that model how
representative they are of their author, i.e., how well engagements with these
posts predict engagements with other posts by that user. Engagements have a
_strength_ attribute to the model quality of the interaction of a user with a
post, e.g. a read may be weaker than a like or retweet. Engagements also have
_error, obs_ and _count_ attributes. The error represents the difference
between the actual strength of engagement and the strength inferred from
engagements with similar posts. This is used for training the posts’ weights
by the gradient descent rule. The Boolean _upd_ and _obs_ attributes signal,
respectively, if an update is required to the training after a change of data
affecting this node and if the strength of the node is observed from real-
world engagement as opposed to being inferred. The _count_ attribute, which is
always 1, is used to calculate the cardinality of Engagement nodes in a multi-
pattern match in Groove. To control training we add an Error node with _error_
and _delta_ attributes. All attributes except _strength_ are part of the
general mechanism while _weight_ also has application-specific meaning.

Figure 1: Rules for dynamic updates: _newUser_ , _newPost_ , and
_newEngagement04_.

The rules come in two flavours: dynamic update rules to create users, posts
and engagements, and rules for training and inference. Update rules are shown
in Fig. 1 in Groove notation. As usual, elements in thin solid black outline
represent read access, dashed blue elements are deleted, and thick green
elements are created, while nodes and edges in dotted red represent negative
conditions. The dotted red line labelled “=” means that the linked nodes are
not equal. Rules to create engagements exist in three versions to create
engagements of strengths $0.2,0.4$ and $0.8$ (but only the rule for $0.4$ is
shown). When creating a new engagement, the user who authored the post has
their _upd_ attribute set to _true_ to trigger retraining.

Figure 2: Top to bottom: Rules _error, delta_ and _infer_.

The training and inference rules are shown in Fig. 2. Nodes labelled $\forall$
or $\forall^{>}$ represent nested multi rules, $@$-labelled edges associating
nodes to their nesting levels and $in$ edges showing quantifier nesting. The
_infer_ rule at the bottom calculates the expected strengths of new
engagements as weighted sums of the strengths of all observed engagements of
the same user with posts by a second user whose _upd_ attribute is _true_ (and
then changes to _false_). The red-green dotted outline says that an Engagement
node is only created if there is not already one, combining node creation with
a negative condition. The assignment shown in green inside that Engagement
node realises the following formula.

$strength=\frac{\sum_{m\in M}(p(m).weight*e(m).strength)}{\mbox{card }M}\quad$
(1)

$M$ is the set of matches of the multi pattern represented by the non-empty
quantification node $\forall^{>}$, with $e(m)$ and $p(m)$ referencing the
Engagement and Post nodes of match $m\in M$. This is nested in another
$\forall$ node to apply the inference to all eligible Engagement nodes at
once. Recall that the weights of Post nodes show how well a post represents
the average post by that user.

The _error_ rule in the top computes the inference error as the difference
between the observed value of strength and the inferred value. This is typical
of semi-supervised learning in GNNs where some of the nodes represent training
samples while others are inferred. Structure and inference formula are similar
to the _infer_ rule. We use an Error node with a global _error_ attribute to
control the training cycle. This is calculated as the sum of absolute _error_
values of all observed engagements, and training stops when the _delta_ , the
difference of global errors between this and the previous cycle is below
$0.0001$. The _delta_ rule in the middle of Fig. 2 implements gradient descent
backward propagation, updating the weights of the posts by the product of the
average error and strengths of the user’s engagements’ with this post.

We use a control program to iterate training and inference, and to simulate
the insertion of newly observed data during the runtime of the model. The
program starts with a training cycle through _error_ and _delta_ for as long
as possible, i.e., until the delta value falls below the threshold. Inference
consists of one application of the _infer_ rule:

    function training(){
       alap{error; delta;}
    }
    function inference(){
       infer;
    }
    function update(){
       node u; newUser(out u);
       node p; newPost(u, out p);
       newEngagement02(p); newEngagement04(p);
    }
    training; inference;
    update;
    training; inference;

All rules use nested universal quantification over engagements, so a single
application per rule covers the entire graph. The sample dynamic update rules
called by the _update_ function create a new user and post and two
engagements. Rule parameters ensure that new posts, engagements and users are
linked. Then we run the training and inference again, demonstrating the
possibility of incremental updates and repeated training to reflect the
interleaving of real-world changes and consequent data updates with the
operation of the GReNN model.

Figure 3: Top to bottom: start, intermediate and final graphs.

The start graph at the top of the Fig. 3 has three users, three posts and six
engagements, all representing “observed” data that can be used to train post
weights. An intermediate graph, after the first training cycle, is shown in
the middle, with final global error and delta values, and updated post
weights. The final graph extends this by, among other things, the structure at
the bottom, where the user on the right of the fragment is the same as the one
in the top left of the intermediate graph. In the final graph, it has two
newly observed and one newly inferred engagement with a new post by the new
user on the left.

##  3 Discussion: GReNN vs GNN

The purpose of this exercise was not to demonstrate how existing GNN
approaches can be realised as graph rewrite systems, but to adopt principles
of GNNs, such as the use of graphs to represent data for training and
inference, to enable semi-supervised learning by labelling parts of the graph
with real-world observations, and to employ the graph structure to direct
training and inference.

However, our model goes beyond mainstream GNNs in two important aspects.
First, the graphs employed are heterogeneous (typed) directed graphs rather
than homogeneous (untyped) undirected ones. This allows us to operate directly
on an application-oriented data model such as a knowledge graph [13] of a
social media network, rather than an encoding. Knowledge graphs are used for
data integration, graph-based analytics and machine learning [20].

“Message-passing type GNNs, also called Message Passing Neural Networks
(MPNN), propagate node features by exchanging information between adjacent
nodes. A typical MPNN architecture has several propagation layers, where each
node is updated based on the aggregation of its neighbours’ features” [4].
MPNNs realise a supervised machine learning algorithm on graph-structured data
by utilising the message-passing algorithm [10]. Notable variants of GNNs
include Convolutional NNs [6], Gated GNNs [16] and Deep Tensor NNs [21].

We relate our model to the concepts and terminology from [10], which describes
the operation of MPNNs on undirected graphs $G$ by message passing, update and
readout steps. Given node features $x_{v}$ and edge features $e_{vu}$, message
passing updates hidden states $h^{t}_{v}$ at each node $v$ for $T$ steps
according to message functions $M_{t}$ and vertex update functions $U_{t}$
using messages

$m^{t+1}_{v}=\sum\limits_{u\in N(v)}M_{t}(h^{t}_{v},h^{t}_{u},e_{vu})$

where $N(v)$ is the set of neighbours of $v$ and updates

$h^{t+1}_{v}=U_{t}(h^{t}_{v},m^{t+1}_{v}).$

After step $T$, readout $R$ computes an output vector for the graph. Hence
messages aggregated from nodes’ neighbourhoods are used in updates to compute
new states. $M_{t}$ and $U_{t}$ can be learned by neural networks. Often this
means that $M_{t}$ computes a weighted sum of node and edge features, with
weights trained by backpropagation, while $U_{t}$ implements activation
thresholds or controls the rate of change by weighing residual states against
massages.

Our Eq. (1) can be seen as a combination of $M_{t}$ and $U_{t}$. The update
function is the identity since no activation or balancing between residual and
new features is required, and the message is a weighted average of engagement
strengths with posts by the same user. However, the analogy is incomplete
because our graphs are directed and not homogeneous, and the strength of
engagement is a feature of an “engaged with post by” relation between users
rather than a user node. Indeed “one could also learn edge features in an MPNN
by introducing hidden states for all edges in the graph …and updating them
analogously” [10]. The relational nature of engagement is the reason why we
only use a single step of forward propagation: All “engagements with posts”
between the same two users are accessible locally in one step. A complex
pattern realises the relation, so we need a sum over the multi pattern to
compute the average. Hence graph patterns allow querying heterogeneous graphs
where their homogeneous encodings are processed by navigating along edges.
Apart from supporting a direct representation of the data, graph patterns
increase the expressive power of GNNs also on homogeneous graphs [3].

Like weights and features, the graph structure is usually represented by
matrices, such as adjacency or incidence and degree matrices. If the graph is
undirected, its adjacency matrix is symmetrical. If the graph is weighted, it
contains real-valued edge weights. A GReNN model supports directed edges,
class heterogeneity and incremental updates to process dynamic application
data in its original form without encoding. As of now, most GNN approaches
only support undirected and homogenous graphs, but there is a trend towards
more general models: In [11], the authors present an example of a directed
MPNN and use it in conjunction with a graph attention network. In [10] the
authors state that “It is trivial to extend the formalism to directed
multigraphs”. Multi-class approaches are utilised in IoT applications such as
[2]. Emerging support for heterogeneous graph data includes [24, 23] which
present a new algorithm called Explicit Message-Passing Heterogeneous Graph
Neural Network (EMP). While GReNNs support both directed and heterogeneous
data at the model level, such approaches may provide a platform to implement
our models efficiently.

Incremental forward and back propagation in GNNs can substantially save time
and resources [8] by reducing computations to nodes affected by updates. This
is especially useful for big-data applications, where batch processing is not
affordable and efficient. NNs and GNNs do not easily support data updates.
Hence, the second generalisation in our approach is to model the dynamic
updates of the graph and control specifically where retraining is required.
This includes adding new nodes, an operation that would break most matrix-
based representations by changing the dimensions of the matrices involved.

GReNNs reimagine GNNs free of considerations of efficient representation of
and computations on graphs. This is appropriate from a modelling perspective,
but unlikely to scale to large graphs. To support an engineering approach to
GNN development based on graph rewriting, a study of the possible
implementations of such models in mainstream GNNs is required. First
experiments show that it is possible in principle to map a model such as ours
to an undirected homogeneous graph by interpreting the edges of that simpler
graph as paths or patterns in our model. This is a reduction that GNN
developers have to perform when they present their application data as input
to a GNN, but we believe that by studying this systematically we can find
standard mappings from a typed attributed graph model into graphs supported by
suitable GNN technology.

The suggestion that graph rewriting can serve as a semantic foundation for
GNNs requires mappings in the opposite direction, taking the various GNN
approaches and modelling them as GReNNs. Then, such embeddings could be used
to compare the different graph types, inference or training rules. The theory
of graph grammars with its language hierarchies and corresponding rule formats
may well play a role here. GNNs can be based on attributed directed or
undirected, simple or multi, homogeneous or typed, binary, hyper or
hierarchical graphs, with little support for reuse of algorithms or
implementations [19]. The categorical theory and potentially generic
implementation of graph rewriting [5] can be a model to unify such diverse
approaches.

## References

  * [1] Allamanis, M., Brockschmidt, M., Khademi, M.: Learning to represent programs with graphs. arXiv:1711.00740 [cs] (May 2018), https://arxiv.org/abs/1711.00740
  * [2] Anjum, A., Ikram, A., Hill, R., Antonopoulos, N., Liu, L., Sotiriadis, S.: Approaching the internet of things (iot): a modelling, analysis and abstraction framework. Concurrency and Computation: Practice and Experience 27(8), 1966–1984 (Oct 2013)
  * [3] Bronstein, M.: Using subgraphs for more expressive GNNs. Medium (Dec 2021), https://towardsdatascience.com/using-subgraphs-for-more-expressive-gnns-8d06418d5ab
  * [4] Bronstein, M.: Beyond message passing: a physics-inspired paradigm for graph neural networks. The Gradient (2022), https://thegradient.pub/graph-neural-networks-beyond-message-passing-and-weisfeiler-lehman
  * [5] Brown, K., Patterson, E., Hanks, T., Fairbanks, J.P.: Computational category-theoretic rewriting. In: Behr, N., Strüber, D. (eds.) Graph Transformation - 15th International Conference, ICGT 2022, Held as Part of STAF 2022, Nantes, France, July 7-8, 2022, Proceedings. Lecture Notes in Computer Science, vol. 13349, pp. 155–172. Springer (2022), https://doi.org/10.1007/978-3-031-09843-7_9
  * [6] Duvenaud, D., Maclaurin, D., Aguilera-Iparraguirre, J., Gómez-Bombarelli, R., Hirzel, T., Aspuru-Guzik, A., Adams, R.P.: Convolutional networks on graphs for learning molecular fingerprints. In: Cortes, C., Lawrence, N.D., Lee, D.D., Sugiyama, M., Garnett, R. (eds.) Advances in Neural Information Processing Systems 28: Annual Conference on Neural Information Processing Systems 2015, December 7-12, 2015, Montreal, Quebec, Canada. pp. 2224–2232 (2015), https://proceedings.neurips.cc/paper/2015/hash/f9be311e65d81a9ad8150a60844bb94c-Abstract.html
  * [7] Eksombatchai, C., Jindal, P., Liu, J.Z., Liu, Y., Sharma, R., Sugnet, C., Ulrich, M., Leskovec, J.: Pixie: A system for recommending 3+ billion items to 200+ million users in real-time. arXiv:1711.07601 [cs] (Nov 2017), https://arxiv.org/abs/1711.07601
  * [8] Galke, L., Vagliano, I., Scherp, A.: Incremental training of graph neural networks on temporal graphs under distribution shift. CoRR abs/2006.14422 (2020), https://arxiv.org/abs/2006.14422
  * [9] Ghamarian, A.H., de Mol, M., Rensink, A., Zambon, E., Zimakova, M.: Modelling and analysis using GROOVE. Int. J. Softw. Tools Technol. Transf. 14(1), 15–40 (2012), https://doi.org/10.1007/s10009-011-0186-x
  * [10] Gilmer, J., Schoenholz, S.S., Riley, P.F., Vinyals, O., Dahl, G.E.: Neural message passing for quantum chemistry. CoRR abs/1704.01212 (2017), http://arxiv.org/abs/1704.01212
  * [11] Han, X., Jia, M., Chang, Y., Li, Y., Wu, S.: Directed message passing neural network (D-MPNN) with graph edge attention (GEA) for property prediction of biofuel-relevant species. Energy and AI 10, 100201 (2022), https://www.sciencedirect.com/science/article/pii/S2666546822000477
  * [12] Heckel, R., Taentzer, G.: Graph Transformation for Software Engineers. Springer International Publishing (2020), http://graph-transformation-for-software-engineers.org/
  * [13] Hogan, A., Blomqvist, E., Cochez, M., d’Amato, C., de Melo, G., Gutiérrez, C., Kirrane, S., Labra Gayo, J.E., Navigli, R., Neumaier, S., Ngonga Ngomo, A.C., Polleres, A., Rashid, S.M., Rula, A., Schmelzeisen, L., Sequeda, J.F., Staab, S., Zimmermann, A.: Knowledge Graphs. No. 22 in Synthesis Lectures on Data, Semantics, and Knowledge, Morgan & Claypool (2021), https://kgbook.org/
  * [14] Lample, G., Charton, F.: Deep learning for symbolic mathematics. arXiv:1912.01412 [cs] (Dec 2019), https://arxiv.org/abs/1912.01412
  * [15] Lange, O., Perez, L.: Traffic prediction with advanced graph neural networks. www.deepmind.com (Sep 2020), https://www.deepmind.com/blog/traffic-prediction-with-advanced-graph-neural-networks
  * [16] Li, Y., Tarlow, D., Brockschmidt, M., Zemel, R.S.: Gated graph sequence neural networks. In: Bengio, Y., LeCun, Y. (eds.) 4th International Conference on Learning Representations, ICLR 2016, San Juan, Puerto Rico, May 2-4, 2016, Conference Track Proceedings (2016), http://arxiv.org/abs/1511.05493
  * [17] Monti, F., Frasca, F., Eynard, D., Mannion, D., Bronstein, M.M.: Fake news detection on social media using geometric deep learning. arXiv:1902.06673 [cs, stat] (Feb 2019), https://arxiv.org/abs/1902.06673
  * [18] Sanchez-Gonzalez, A., Godwin, J., Pfaff, T., Ying, R., Leskovec, J., Battaglia, P.: Learning to simulate complex physics with graph networks (Nov 2020), https://proceedings.mlr.press/v119/sanchez-gonzalez20a.html
  * [19] Sanchez-Lengeling, B., Reif, E., Pearce, A., Wiltschko, A.: A gentle introduction to graph neural networks. Distill 6(8) (Aug 2021)
  * [20] Schad, J.: Graph powered machine learning: Part 1. ML Conference Berlin (October 2021), https://mlconference.ai/ml-summit/
  * [21] Schütt, K.T., Arbabzadah, F., Chmiela, S., Müller, K.R., Tkatchenko, A.: Quantum-chemical insights from deep tensor neural networks. Nature Communications 8(1) (jan 2017), https://doi.org/10.1038%2Fncomms13890
  * [22] Stokes, J.M., Yang, K., Swanson, K., Jin, W., Cubillos-Ruiz, A., Donghia, N.M., MacNair, C.R., French, S., Carfrae, L.A., Bloom-Ackerman, Z., Tran, V.M., Chiappino-Pepe, A., Badran, A.H., Andrews, I.W., Chory, E.J., Church, G.M., Brown, E.D., Jaakkola, T.S., Barzilay, R., Collins, J.J.: A deep learning approach to antibiotic discovery. Cell 180(4), 688–702.e13 (Feb 2020), https://www.cell.com/cell/fulltext/S0092-8674(20)30102-1
  * [23] Wang, Z., Wang, C., Zhao, S., Xu, Y., Hao, S., Hsieh, C.Y., Gu, B., Duan, W.: Heterogeneous relational message passing networks for molecular dynamics simulations. CoRR abs/2109.00711 (2021), https://arxiv.org/abs/2109.00711
  * [24] Xu, L., He, Z.Y., Wang, K., Wang, C.D., Huang, S.Q.: Explicit message-passing heterogeneous graph neural network. IEEE Transactions on Knowledge and Data Engineering pp. 1–13 (2022)
  * [25] Zhou, J., Cui, G., Hu, S., Zhang, Z., Yang, C., Liu, Z., Wang, L., Li, C., Sun, M.: Graph neural networks: A review of methods and applications. AI Open 1, 57–81 (2020), https://www.sciencedirect.com/science/article/pii/S2666651021000012