Deep Set Autoencoders for Anomaly Detection in Particle Physics

Bryan Ostdiek

Department of Physics, Harvard University, Cambridge, MA 02138, USA

The NSF AI Institute for Artificial Intelligence and Fundamental Interactions

<EMAIL_ADDRESS>
## Abstract

There is an increased interest in model agnostic search strategies for physics
beyond the standard model at the Large Hadron Collider. We introduce a Deep
Set Variational Autoencoder and present results on the Dark Machines Anomaly
Score Challenge. We find that the method attains the best anomaly detection
ability when there is no decoding step for the network, and the anomaly score
is based solely on the representation within the encoded latent space. This
method was one of the top-performing models in the Dark Machines Challenge,
both for the open data sets as well as the blinded data sets.

###### Contents

  1. 1 Introduction
  2. 2 The Datasets
  3. 3 Methodology
    1. 3.1 The Network
    2. 3.2 Training
  4. 4 Results
    1. 4.1 Explorations on channels 1 and 2
    2. 4.2 Application to Channel 3
  5. 5 Discussion

##  1 Introduction

Despite extensive searches for new physics, the experiments at the Large
Hadron Collider (LHC) have not yet found physics Beyond the Standard Model
(BSM). While there are many theoretically motivated BSM models, the lack of
significant evidence for these models motivates studies of broad, model-
agnostic searches. The rise of modern machine learning tools has generated
many such data driven methods to search for new physics without a specific
model, see Ref. [1] for a recent review of anomaly detection or unsupervised
methods.

We break the types of anomaly detection techniques into two categories. The
first is methods that use density estimation or look for bumps on smooth
curves [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24]. The recent LHC Olympics [17] introduced four data sets of
simulated hadronic LHC events. Many teams submitted results for new methods
looking for hidden resonances in the data. The second category of anomaly
detection searches for individual events that look different than what is
expected from the standard model (also known as out-of-distribution events)
[25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43,
44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54]. The Dark Machines Anomaly Score
Challenge [50] studied many methods to flag individual events as anomalous
after training on a simulated set of standard model background events. This
work was a contributed method to the Dark Machines challenge.

In this paper, we introduce the Deep Set $\beta$ Variational Autoencoder for
anomaly detection at the LHC. It is built upon the principle that the
reconstructed objects at the LHC can be viewed as sets of particles, or a
point cloud. In the mathematical sense, a set is a collection of objects; in
our case, the objects are the reconstructed particles. We build upon the
particle flow deep set networks introduced in Ref. [55] for supervised
classification.111Graphs have also been used in association with the point
cloud representation of the data [56]. Similarly, deep set networks have been
studied in the context of self-supervision in Ref. [57]. The particle flow
networks were built on two functions represented as neural networks. The first
function, called $\Phi$, operates on each particle within the set in parallel,
mapping from the number of input features per particle (four momentum, charge,
particle id, etc) to a defined latent output dimension. The outputs for each
particle are then summed along the number of particles to make the latent
representation of the event permutation invariant. The second function, called
$F$, maps from the latent space to the final answer, which for Ref. [55] was a
number between 0 and 1 for quark-jet versus gluon-jet discrimination. To
implement this method for anomaly detection, we adapt an autoencoder approach,
where the original data is compressed to a small latent space and then
remapped back to reconstruct the input data; anomalous events are expected to
have larger discrepancies between the input and output than the in-
distribution data which is used to train the network. Thus, after constructing
a permutation invariant latent space, we then use $F$ to generate a new set.
To compare the input and output sets, we use the Chamfer loss, which is a
method to measure the distance between sets [58].

To force a structure on the latent space, we add a variational step, such that
the output of $\Phi$ controls the parameters of a random draw for the
representation in latent space. An extra term is added to the loss function
which drives the latent representation towards a multidimensional normal
distribution. The parameter $\beta$ will control the relative importance of
the set reconstruction loss and the latent representation loss,

$\text{loss}=\beta\times(\text{latent representation
loss})+\big{(}1-\beta\big{)}\times(\text{set reconstruction loss})$ (1)

with $\beta=0$ focusing only on the reconstruction loss and $\beta=1$ focusing
only on the latent representation. When evaluating events in the test set, we
use the loss as the metric to determine how anomalous the event is, we expect
lower loss for SM events and larger loss for BSM events. Unexpectedly, we find
that $\beta=1$ offers the best anomaly detection performance for this method.
In fact, our networks with $\beta=1$ were among the top methods submitted to
the Dark Machines challenge [50]. As $\beta=1$ does not include the
reconstruction loss, there is actually no need for the decoder part of the
network, and the method is only trying to map the set to a Gaussian latent
space. The other top methods of the challenge used density estimation
techniques combined with either a variational autoencoder with $\beta=1$ or a
Deep Support Vector Data Description (Deep SVDD) [59] model. The Deep SVDD
models work by trying to map the input data to a single fixed vector. The SM
background events are expected to get mapped very close to the vector, while
BSM events will be further away. Thus, all of the top models in the Dark
Machines challenge used methods which either map the input data to a fixed
Gaussian distribution or a fixed vector. It is unclear why these so called
“fixed target” methods are outperforming models which reconstruct the input.

The paper is organized as follows. In Sec. 2, we briefly review the Dark
Machines datasets used for this study. We introduce our network architecture
and training details in Sec. 3. The networks are applied to independent test
sets in Sec. 4. Conclusions are presented in Sec. 5.

##  2 The Datasets

The data for this work come from the Dark Machines Anomaly Score Challenge
[50]. A large dataset of over 1 billion simulated proton-proton collisions
with a center of mass energy of $\sqrt{S}=13$ TeV was generated [60]. This set
of events comes from 26 different standard model (SM) processes. From these
events, four potential signal regions are examined.

  * •

Channel 1: is focused on events with a lot of hadronic activity and missing
energy. The selection requirements are

$\displaystyle H_{T}\geq 600\text{ GeV},$ $\displaystyle E_{T}^{\rm{miss}}\geq
200\text{ GeV,}$ $\displaystyle E_{T}^{\rm{miss}}/H_{T}\geq 0.2,$

along with at least four jets (including $b$-jets) with $p_{T}\geq 50$ GeV
including at least one with $p_{T}\geq 200$ GeV. There are $214\,000$ SM
events in this channel.

  * •

Channel 2: focuses on leptonic-events. There are two sub channels.

    * –

a: contains a least 3 leptons (muons/electrons) with $p_{T}\geq 15$ GeV and at
least 50 GeV of $E_{T}^{\rm{miss}}$. There are only $20\,000$ SM events in
this channel.

    * –

b: removes the necessity of one of the leptons and replaces it with hadronic
activity. The requirements are at least two leptons with $p_{T}\geq 15$ GeV,
$E_{T}^{\rm{miss}}\geq 50$ GeV, and $H_{T}\geq 50$ GeV. There are $340\,000$
SM events in this channel.

  * •

Channel 3: is more inclusive and contains events which have $H_{T}\geq 600$
GeV and $E_{T}^{\rm{miss}}\geq 100$ GeV. This channel contains $8\,500\,000$
background SM events.

In addition to the SM events, events were generated for 11 different BSM
signals for a total of 18 mass spectra. The BSM models contain a dark matter
particle and come from general classes of $Z^{\prime}$, R parity conserving
supersymmetry, and R parity violating supersymmetry (RPV). While the BSM
events for each signal may appear in multiple channels, for the most part we
focus on anomaly detection within each channel separately.

The events for each channel are available on
Zenodo.222https://zenodo.org/record/3961917 [61] The data consists of the
reconstructed physics objects: non-$b$-tagged jets, $b$-tagged jets, $e^{+}$,
$e^{-}$, $\mu^{+}$, $\mu^{-}$, and $\gamma$. For each object in the event, the
four-momentum and object type are available. In addition, the missing
transverse momentum and its azimuthal direction are recorded for each event.
There were a maximum of 20 objects in each event.

The task of the challenge was to train an anomaly detection method on the SM
background events (leaving the final 10% for testing). For any given event, we
produce an “anomaly score” which denotes how close or far away the event is
from the SM expectation.

##  3 Methodology

###  3.1 The Network

The main idea behind our method is that at the outgoing particles in collider
experiments can be thought of as a collection of four-vectors. There is no
intrinsic ordering, although they are often sorted by the magnitude of the
transverse momentum. Using a network architecture which respects the
permutation symmetry is more natural and could lead to improved results. In
Ref. [62], it was shown that “Deep Sets” with permutation-invariant functions
of variable-length inputs can be parameterized in a fully general way. This
idea was introduced to the High Energy Physics community in Ref. [55], where
the operations were generalized to include infrared and collinear (IRC)
safety. Their methods outperformed other state-of-the art classifiers, and
they found a slight improvement if the operations were not restricted to be
IRC safe.

As an unsupervised learning task, we modify the deep sets paradigm to include
an auto encoding structure, following the example of Ref. [58]. As in [55], we
map each particle to the latent space using a function $\Phi$. The functional
form of $\Phi$ is a four layer neural network which has inputs of
$(\log_{10}E,\log_{10}p_{T},\eta,\phi,\text{ID})$. There are then three hidden
layers with 264 nodes each, using ReLU activation functions. The final layer
of $\Phi$ contains 16 nodes, containing the mean ($\mu_{i}$) and variance
($\sigma_{i}$) for each of the eight ($i$) latent dimensions. This is depicted
in left panel of Fig. 1. The right panel illustrates how the same function is
applied to each particle within the event, leading to many latent
representations.

Figure 1: The left panel displays the network $\Phi$ which takes as input the
four-vector and particle identification of a particle and maps it into the
latent space. The right panel shows that how this is done for each particle in
the event. We use a featurewise sort pooling layer to combine the per-particle
latent space into an event-level representation as shown in Eq. (2).

One possible way of combining the per-particle latent spaces into an event-
level latent representation is to sum the individual components [58]. However,
Ref. [63] showed that the performance can be improved by sorting each feature
(for instance all of the $\mu_{1}$) along all of the particles followed by a
learned mapping from the sorted features to the latent space. As an example
after the first feature ($\mu_{1}$) is sorted per-particle, it is mapped to
the event-level representation as

$\includegraphics{figures/FSPOOL},$ (2)

where $\mu_{1}^{\text{max}(\text{min})}$ is the maximum (minimum) over the
particles. This layer/operation is known as featurewise sort pooling (FSPool)
[63].

After the FSPooling layer, we reparameterize as a variational autoencoder with
a Gaussian prior. The event-level standard deviation is given as

$\sigma_{1}\rightarrow e^{\frac{\sigma_{1}}{2}}.$ (3)

Then a random representation is drawn using

$z_{1}=\mathcal{G}(\mu_{1},\sigma_{1}),$ (4)

where $\mathcal{G}$ is a Gaussian.

A decoding network is then used to transform the latent data back to a set of
four-vectors and particle IDs. We do so using a dense neural network as shown
in Fig. 2. Two layers with 256 nodes with ReLU activations are used in the
network. The network then splits into two parts. In the first part, there are
80 nodes, representing 20 possible four-vectors. The second part is a series
of 180 nodes representing the probability that each of the 20 particles
belongs to one of the eight particle IDs or the probability that the four
vector should not be included as member of the final set.

Figure 2: The decoder takes latent space representations of events and outputs
a set of particles. The blue box denotes shared wights between the two final
output components. The orange box is a set of four-vectors for each of the
predicted particles. The pink box represents the classification prediction
which assigns a probability to each particle to belong to a certain class or
to be masked out.

Figure 3 shows an illustrative example representation of an input set (left)
and encoded-decoded set (right). Each panel is a two dimensional slice across
the four-vector representation $(\log_{10}E$, $\log_{10}p_{T}$, $\eta$,
$\phi)$ and the points denote the ojects. The color of the points denotes the
type of the object. We include the missing transverse momentum as an “object”
with zero energy and no $\eta$ information. In the output set, we color the
objects based on the highest classification score. This leads to sets that do
not make sense from a physics perspective (there is more than one missing
transverse momentum object). While this is just a single example, we find that
most of the reconstructed sets have similar problem, with too many output
particles and missing energy objects in the set. If we were worried about
event generation, rather than anomaly detection, procedures could be
implemented, such as having a specific output for the missing transverse
momentum rather than including it as generic object, or only including objects
if the classification score is above some threshold. While the output sets may
not look physical, we find that we are still able to achieve decent anomaly
detection.

Figure 3: Illustrative example of the input and output structure of the
network. The left panels show the original event in different two dimensional
slices. The color of each point denotes the type of object. This event is from
Channel 2b and contains two jets, a muon pair, and missing transverse
momentum. The right panels contain the encoded-decoded set, and there are more
objects than in the original set.

###  3.2 Training

To train the deep set variational autoencoder, the output four-vectors and
class predictions are compared with the input set of four-vectors and particle
IDs. This is challenging for sets on full events for two reasons. The first is
that by construction, the input and output are permutation invariant so we
cannot just compare the first element of the input to the first element of the
output. One solution to this challenge would be to use a metric such as the
Energy Mover’s distance [64] which is based on optimal transport and looks at
how much energy is needed to reshuffle one event to look like a different
event. However, this brings about the second challenge, what is the cost to
change particle types?

With no clear answer to this, we instead implement a modified version of the
Chamfer loss. For each particle in the input set, we find the Euclidian
distance (in the space of $(\log_{10}E,\log_{10}p_{T},\eta,\phi)$) to the
closest particle in the output set. To this distance, we then add the log
likelihood that the particle in the output is classified as the same particle
type. For this, we use the log of the softmax of the output of the classifier.
At this stage, if there are any particles in the output set that are not close
to any of the input particles, they may not be included in the loss because
only the closest particles are used. To account for this, the Chamfer loss
then repeats the procedure for every particle in the output set, finding the
closest particle in the input set. We again supplement this distance by adding
by the probability that they are the same class. It is unclear how to weight
the relative importance of the four-vectors being close compared the particle
ID being assigned correctly, therefore, we include an extra hyper-parameter,
$w$ which could be optimized over. Our modified Chamfer loss is then given by,

$\displaystyle L_{C}=$ $\displaystyle\sum_{x\in
S_{\rm{input}}}\Big{(}\underset{y\in
S_{\rm{output}}}{\text{min}}\big{|}\vec{x}-\vec{y}\big{|}\Big{)}+w\,\log
p(x_{\rm{ID}}=y_{\rm{ID}}^{\rm{min}})$ (5) $\displaystyle+\sum_{y\in
S_{\rm{output}}}\Big{(}\underset{x\in
S_{\rm{input}}}{\text{min}}\big{|}\vec{x}-\vec{y}\big{|}\Big{)}+w\,\log
p(x_{\rm{ID}}^{\rm{min}}=y_{\rm{ID}})~{}.$

The last component of the loss for the variational autoencoder enforces
structure in the latent space. This is done using the reparameterization trick
and computing the KL divergence between the latent space and a Gaussian prior
[65, 66]. This term helps to regularize the network and put similar events
nearby in latent space. However, it is unclear how much regularization should
be included. Therefore, we include the hyper-parameter $\beta$ to control the
relative importance of the KL term and the chamfer loss as,

$L_{\rm{total}}=\beta\times KL+\big{(}1-\beta)\times L_{C}~{}.$ (6)

The network is implemented in PyTorch [67] and trained using the Adam
optimizer [68]. We use the default parameters for Adam, but lower the learning
rate by a factor of 10 if the loss computed on the validation set has not
improved for 2 epochs. Training ends when the validation loss has not improved
for 4 epochs, and typically takes around 30 epochs. We use a batch size of 500
events. We note that the loops required in computing the Chamfer loss makes
training quite slow. Because of this, we were unable to explore a huge range
of hyper-parameter options. Despite this, our deep set autoencoders were among
the best models in the Dark Machines Challenge [50].

##  4 Results

In this work, we experiment with different values of the regularization,
$\beta$, and the importance of the classification step of the decoder, $w$. We
scan over

$\displaystyle\beta\in$
$\displaystyle\left\\{10^{-6},0.001,0.1,0.5,0.8,0.999,1.0\right\\}\text{ and
}$ (7) $\displaystyle w\in$ $\displaystyle\left\\{1.0,10.0,100.0\right\\}~{}.$

The number of events in channel 3 is so large that training the network was
prohibitively slow, and could not be performed for the entire scan (see Sec.
2). Instead, we treat channels 1, 2a, and 2b as open data sets to optimize the
values of $\beta$ and $w$ and use channel 3 as a final check to test for
generalization.

Figure 4:  Anomaly score for the networks using $\beta=10^{-3}$ and $w=10.0$.
The value of $\beta$ determines the ratio of contributions to the loss coming
from the set reconstruction and the latent space representation and $w$
controls the relative importance of the four-vector reconstruction versus the
particle identification. The networks are trained on SM background events and
then applied to a new set of background and BSM events. The background
histograms are stacked, while the BSM signals are independently overlaid. Even
though the number of dominant background changes for each channel, the same
network architecture can be applied to each to separate some signals from the
background.

###  4.1 Explorations on channels 1 and 2

We use the full loss term, including the KL divergence, as our event-by-event
anomaly score. This is found to yield good separation between the SM
background and many of the signals. For example, Fig. 4 shows the total loss
of events from the test set for $\beta=0.001$ and $w=10.0$. The colored bars
are stacked histograms for the SM background as indicated by the legends. Each
channel has different dominant backgrounds. Overlayed on the plots are
histograms for three different signals, denoted by the solid, dashed, and
dotted black lines. In channel 1, these are RPV stop production, R parity
conserving stop production, and a mono-top $Z^{\prime}$ model, respectively.
For channel 2 (both and A and B) we display signal events for R parity
conserving gluino pair production, a lepton number violating $Z^{\prime}$, and
R parity conserving chargino-neutralino production. Even though the sets are
not reconstructed well (for example see Fig. 3), we see that some of the new
physics events have larger loss values than the SM events. In each channel,
there is at least one BSM model which is easy to separate from the background,
the RPV stops for channel 1 and gluino production for channel 2. However, we
also note that there are also BSM models that the network has a hard time
distinguishing from the SM, for instance the mono-top $Z^{\prime}$ for channel
1 and the chargino-neutralino production for channel 2.

In Fig. 5, we show the similar loss histograms for the network using
$\beta=1.0$ and $w=1.0$. This value for $\beta$ essentially removes the
decoder from the network, and the network is only trying to reduce the set to
an eight dimensional Gaussian. The same BSM models are chosen as in Fig. 4 to
aid in comparison. In Channel 1, we see that both the
$\tilde{t}\tilde{t}^{*}\rightarrow t\bar{t}\tilde{\chi}^{0}\tilde{\chi}^{0}$
signal and the $Z^{\prime}+t$ signal are much more separated from the
background when the network does not have the decoder. A similar pattern is
observed in the channels 2, with more of the hard to discover signals getting
separated from the background.

Figure 5: Anomaly score for the networks using $\beta=1.0$ and $w=1.0$. The
network does not try to reconstruct the set, but only uses the KL divergence
in the latent space. This leads to better signal separation than the networks
which use a decoder.

With many unique BSM signal-Channel combinations, determining the best values
of $\beta$ and $w$ is not trivial. Sometimes networks find certain signals
better than others, but changing the parameters can trade which signals can be
separated from the background. We incorporate a holistic view of the problem
as done in the challenge [50, 69]. For each channel, we find the total loss
which allows 1%, 0.1%, or 0.01% of the background events through. We then
determine the amount of each signal that passes each of these loss thresholds
to compute the significance improvement, given by

$\text{Significance Improvement}=\epsilon_{S}/\sqrt{\epsilon_{B}},$ (8)

where $\epsilon_{S}$ is the signal efficiency and $\epsilon_{B}$ is the
background efficiency. For each signal in each channel, we keep the threshold
which maximizes the significance improvement. In some cases, the looser cut
allows through enough signal to be better than a smaller cut, while other
times the tighter cut is necessary to remove a larger portion of the
background. Some of the BSM physics models have signals which would appear in
multiple channels. When this is the case, we denote the total significance
improvement (TI) as the largest significance improvement across the channels
for that signal.

Figure 6 shows the distribution of the TI scores across the 11 BSM models (18
total mass spectra) for the different networks. The color of the marker
denotes the value of $w$ and value of $\beta$ is denoted by the region along
the $y$-axis. The boxes cover the 25th-75th quantile, with the line denoting
the median of the data. The line extend out to either the maximum or minimum,
or to 1.5 times the box width from the end of the box, in which case the
remaining data are shown as (outlier) points. The minimum TI score is less
than 1 for all models, which implies that there is some BSM model which all of
the networks cannot distinguish from the SM background.

There appears to be a clear distinction between the anomaly detectors that use
the decoder ($\beta\neq 1$) and the method without the decoder ($\beta=1$).
The distribution of TI scores are much larger for the model without the
decoder. For instance, the maximum scores are around a factor of 3-4 better,
and the median scores are nearly an order of magnitude better. The cause of
this is unknown and a subject of future research. Going forward, we will only
study the model with $\beta=1$ in this work.

###  4.2 Application to Channel 3

Figure 6: Distribution of total improvement scores from the signals of
channels 1, 2a, and 2b. Changing the value of $w$ (denoted by the colors) has
a minor impact. In contrast, setting $\beta=0$ removes the decoder and yields
significantly better anomaly detection performance.  Figure 7: Significance
improvements of the BSM signals for the models using $\beta=1.0$. The maximum
significance improvement over three anomaly score cuts are displayed for each
channel. The circle and squares mark the networks with $w=1.0$ and $w=10.0$,
respectively. However, $w$ does not enter into the loss when $\beta=1$, so the
spread between the results can be viewed as a training uncertainty.

After determining that the networks which only compress to a small latent
space have larger median and maximum total improvements on channels 1, 2a, and
2b, we train these models on channel 3. The resulting significance
improvements are shown for each BSM physics signal in Fig. 7. The circle and
square markers denote the network trained with $w=1$ and $w=10$, respectively.
However, the parameter $w$ should not affect the network when $\beta=1$, so
the differences indicate the level of training uncertainty. The blue, orange,
green, and red markers show the results for channel 1, 2a, 2b, and 3,
respectively. While channel 3 has the most training data, this did not lead to
overall better anomaly detection. In only one BSM signal (monojet +
$Z^{\prime}\rightarrow$ dark matter) is the significance improvement greater
for channel 3 than the other channels.

Overall, we find that the networks are able to detect most of the BSM models.
We note that the performance is worst for channel 2a, which was expected
because the training set is quite small. Given this, the networks are not
sensitive to the two supersymetric chargino-neutralino production signals
which only appear in channel 2a. All of the other signals which appear in
channel 2a also appear in 2b, and the significance improvement is between
5-10.

We find that the best significance improvement tends to come from channel 2b,
which is the two lepton channel. At this stage, it is unclear if the
performance is better because of the size of the training data, or if the BSM
signals that appear in this channel just happen to be easier to detect.
Investigating this is left for future work.

##  5 Discussion

We have developed an anomaly detection method for use at the LHC as a model
agnostic BSM search. The method uses a variational autoencoder setup, where
the input data is compressed to a small latent space and then decompressed and
compared with the original input. We use a deep set input data architecture,
where each reconstructed particle is viewed as an independent point in the
set, and the operations are permutation invariant.

In comparing the input and output sets, we include a term to account for the
particle type, and scan over how much to weight this term in comparison the
the four momentum reconstruction. In general, we find that placing more
emphasis on the particle identification leads to better anomaly detection
performance. This is especially important if one wants to use the deep set VAE
as a generative model.

Surprisingly, we find that the network trained without a decoder provides the
best anomaly detection. This is unintuitive because the network is only trying
to compress the data, and it does not need to do so in a way that allows for
good decompression. With no decoder or decompression, we can say that the
network is mapping the input data to a fixed representation (a
multidimensional Gaussian). In the Dark Machines challenge [50], these deep
set VAE models were among the top models, achieving median total improvements
greater than 2 in both the hackathon datasets (results from this work) and in
the secret dataset. Interestingly, the other set of models which had high
median scores for both the hackathon and secret data sets also used as part of
their anomaly score a mapping to a fixed representation [51]. The fixed
representations were either also a multidimensional Gaussian distribution or a
fixed vector. These type of networks were called fixed target methods in the
Dark Machines Challenge.

In this work, and the Dark Machines challenge as a whole, the anomaly
detection models were trained on simulation and applied to simulation using
the same underlying parameters. The models which do not use a reconstruction
loss, but solely focus on the compression of data to the latent space may be
harder to validate when applying to real LHC data. To our knowledge, these
fixed target methods have not been applied to single object studies, such as
anomalous jet tagging, it would be interesting to see if the fixed target
aspect also works well there, or if the methods are better for full event
analyses. These questions are left for future study.

## Acknowledgements

We thank Baptiste Ravina and Nathan Simpson for early collaboration on this
work. We thank Katie Fraser, Samuel Homiller, Rashmish Mishra, and Sascha
Caron for insightful comments on a previous version of this paper. The work of
B.O. is supported by the U.S. Department of Energy under contracts DE-
SC0013607 and DE-SC0020223. This work is supported by the National Science
Foundation under Cooperative Agreement PHY-2019786 (The NSF AI Institute for
Artificial Intelligence and Fundamental Interactions, http://iaifi.org/).
Computations in this paper were run on the FASRC Cannon cluster supported by
the FAS Division of Science Research Computing Group at Harvard University.

## References

  * [1] B. Nachman,  _Anomaly Detection for Physics Analysis and Less than Supervised Learning_ (2020),  2010.14554.
  * [2] J. H. Collins, K. Howe and B. Nachman,  _Anomaly Detection for Resonant New Physics with Machine Learning_ ,  Phys. Rev. Lett. 121(24), 241803 (2018),  10.1103/PhysRevLett.121.241803,  1805.02664.
  * [3] R. T. D’Agnolo and A. Wulzer,  _Learning New Physics from a Machine_ ,  Phys. Rev. D 99(1), 015014 (2019),  10.1103/PhysRevD.99.015014,  1806.02350.
  * [4] A. De Simone and T. Jacques,  _Guiding New Physics Searches with Unsupervised Learning_ ,  Eur. Phys. J. C 79(4), 289 (2019),  10.1140/epjc/s10052-019-6787-3,  1807.06038.
  * [5] A. Casa and G. Menardi,  _Nonparametric semisupervised classification for signal detection in high energy physics_ ,  arXiv e-prints arXiv:1809.02977 (2018),  1809.02977.
  * [6] J. H. Collins, K. Howe and B. Nachman,  _Extending the search for new resonances with machine learning_ ,  Phys. Rev. D 99(1), 014038 (2019),  10.1103/PhysRevD.99.014038,  1902.02634.
  * [7] B. M. Dillon, D. A. Faroughy and J. F. Kamenik,  _Uncovering latent jet substructure_ ,  Phys. Rev. D 100(5), 056002 (2019),  10.1103/PhysRevD.100.056002,  1904.04200.
  * [8] A. Mullin, S. Nicholls, H. Pacey, M. Parker, M. White and S. Williams,  _Does SUSY have friends? A new approach for LHC event analysis_ ,  JHEP 02, 160 (2021),  10.1007/JHEP02(2021)160,  1912.10625.
  * [9] R. T. D’Agnolo, G. Grosso, M. Pierini, A. Wulzer and M. Zanetti,  _Learning multivariate new physics_ ,  Eur. Phys. J. C 81(1), 89 (2021),  10.1140/epjc/s10052-021-08853-y,  1912.12155.
  * [10] B. Nachman and D. Shih,  _Anomaly Detection with Density Estimation_ ,  Phys. Rev. D 101, 075042 (2020),  10.1103/PhysRevD.101.075042,  2001.04990.
  * [11] A. Andreassen, B. Nachman and D. Shih,  _Simulation Assisted Likelihood-free Anomaly Detection_ ,  Phys. Rev. D 101(9), 095004 (2020),  10.1103/PhysRevD.101.095004,  2001.05001.
  * [12] G. Aad _et al._ ,  _Dijet resonance search with weak supervision using $\sqrt{s}=13$ TeV $pp$ collisions in the ATLAS detector_,  Phys. Rev. Lett. 125(13), 131801 (2020),  10.1103/PhysRevLett.125.131801,  2005.02983.
  * [13] B. M. Dillon, D. A. Faroughy, J. F. Kamenik and M. Szewc,  _Learning the latent structure of collider events_ ,  JHEP 10, 206 (2020),  10.1007/JHEP10(2020)206,  2005.12319.
  * [14] K. Benkendorfer, L. L. Pottier and B. Nachman,  _Simulation-assisted decorrelation for resonant anomaly detection_ ,  Phys. Rev. D 104(3), 035003 (2021),  10.1103/PhysRevD.104.035003,  2009.02205.
  * [15] V. Mikuni and F. Canelli,  _Unsupervised clustering for collider physics_ ,  Phys. Rev. D 103(9), 092007 (2021),  10.1103/PhysRevD.103.092007,  2010.07106.
  * [16] G. Stein, U. Seljak and B. Dai,  _Unsupervised in-distribution anomaly detection of new physics through conditional density estimation_ ,  In _34th Conference on Neural Information Processing Systems_ (2020), 2012.11638.
  * [17] G. Kasieczka _et al._ ,  _The LHC Olympics 2020: A Community Challenge for Anomaly Detection in High Energy Physics_ (2021),  2101.08320.
  * [18] J. Batson, C. G. Haaf, Y. Kahn and D. A. Roberts,  _Topological Obstructions to Autoencoding_ ,  JHEP 04, 280 (2021),  10.1007/JHEP04(2021)280,  2102.08380.
  * [19] A. Blance and M. Spannowsky,  _Unsupervised Event Classification with Graphs on Classical and Photonic Quantum Computers_ (2021),  2103.03897.
  * [20] B. Bortolato, B. M. Dillon, J. F. Kamenik and A. Smolkovič,  _Bump Hunting in Latent Space_ (2021),  2103.06595.
  * [21] J. H. Collins, P. Martín-Ramiro, B. Nachman and D. Shih,  _Comparing weak- and unsupervised methods for resonant anomaly detection_ ,  Eur. Phys. J. C 81(7), 617 (2021),  10.1140/epjc/s10052-021-09389-x,  2104.02092.
  * [22] T. Dorigo, M. Fumanelli, C. Maccani, M. Mojsovska, G. C. Strong and B. Scarpa,  _RanBox: Anomaly Detection in the Copula Space_ (2021),  2106.05747.
  * [23] S. Volkovich, F. De Vito Halevy and S. Bressler,  _The Data-Directed Paradigm for BSM searches_ (2021),  2107.11573.
  * [24] A. Hallin, J. Isaacson, G. Kasieczka, C. Krause, B. Nachman, T. Quadfasel, M. Schlaffer, D. Shih and M. Sommerhalder,  _Classifying Anomalies THrough Outer Density Estimation (CATHODE)_ (2021),  2109.00546.
  * [25] J. A. Aguilar-Saavedra, J. H. Collins and R. K. Mishra,  _A generic anti-QCD jet tagger_ ,  JHEP 11, 163 (2017),  10.1007/JHEP11(2017)163,  1709.01087.
  * [26] J. Hajer, Y.-Y. Li, T. Liu and H. Wang,  _Novelty Detection Meets Collider Physics_ ,  Phys. Rev. D 101(7), 076015 (2020),  10.1103/PhysRevD.101.076015,  1807.10261.
  * [27] T. Heimel, G. Kasieczka, T. Plehn and J. M. Thompson,  _QCD or What?_ ,  SciPost Phys. 6(3), 030 (2019),  10.21468/SciPostPhys.6.3.030,  1808.08979.
  * [28] M. Farina, Y. Nakai and D. Shih,  _Searching for New Physics with Deep Autoencoders_ ,  Phys. Rev. D 101(7), 075021 (2020),  10.1103/PhysRevD.101.075021,  1808.08992.
  * [29] O. Cerri, T. Q. Nguyen, M. Pierini, M. Spiropulu and J.-R. Vlimant,  _Variational Autoencoders for New Physics Mining at the Large Hadron Collider_ ,  JHEP 05, 036 (2019),  10.1007/JHEP05(2019)036,  1811.10276.
  * [30] T. S. Roy and A. H. Vijay,  _A robust anomaly finder based on autoencoders_ (2019),  1903.02032.
  * [31] A. Blance, M. Spannowsky and P. Waite,  _Adversarially-trained autoencoders for robust unsupervised new physics searches_ ,  JHEP 10, 047 (2019),  10.1007/JHEP10(2019)047,  1905.10384.
  * [32] M. Romão Crispim, N. F. Castro, R. Pedro and T. Vale,  _Transferability of Deep Learning Models in Searches for New Physics at Colliders_ ,  Phys. Rev. D 101(3), 035042 (2020),  10.1103/PhysRevD.101.035042,  1912.04220.
  * [33] O. Amram and C. M. Suarez,  _Tag N’ Train: a technique to train improved classifiers on unlabeled data_ ,  JHEP 01, 153 (2021),  10.1007/JHEP01(2021)153,  2002.12376.
  * [34] M. Crispim Romão, N. F. Castro, J. G. Milhano, R. Pedro and T. Vale,  _Use of a generalized energy Mover’s distance in the search for rare phenomena at colliders_ ,  Eur. Phys. J. C 81(2), 192 (2021),  10.1140/epjc/s10052-021-08891-6,  2004.09360.
  * [35] O. Knapp, O. Cerri, G. Dissertori, T. Q. Nguyen, M. Pierini and J.-R. Vlimant,  _Adversarially Learned Anomaly Detection on CMS Open Data: re-discovering the top quark_ ,  Eur. Phys. J. Plus 136(2), 236 (2021),  10.1140/epjp/s13360-021-01109-4,  2005.01598.
  * [36] M. Crispim Romão, N. F. Castro and R. Pedro,  _Finding New Physics without learning about it: Anomaly Detection as a tool for Searches at Colliders_ ,  Eur. Phys. J. C 81(1), 27 (2021),  10.1140/epjc/s10052-020-08807-w,  2006.05432.
  * [37] T. Cheng, J.-F. Arguin, J. Leissner-Martin, J. Pilette and T. Golling,  _Variational Autoencoders for Anomalous Jet Tagging_ (2020),  2007.01850.
  * [38] C. K. Khosa and V. Sanz,  _Anomaly Awareness_ (2020),  2007.14462.
  * [39] P. Thaprasop, K. Zhou, J. Steinheimer and C. Herold,  _Unsupervised Outlier Detection in Heavy-Ion Collisions_ ,  Phys. Scripta 96(6), 064003 (2021),  10.1088/1402-4896/abf214,  2007.15830.
  * [40] J. A. Aguilar-Saavedra, F. R. Joaquim and J. F. Seabra,  _Mass Unspecific Supervised Tagging (MUST) for boosted jets_ ,  JHEP 03, 012 (2021),  10.1007/JHEP03(2021)012,  [Erratum: JHEP 04, 133 (2021)],  2008.12792.
  * [41] A. A. Pol, V. Berger, G. Cerminara, C. Germain and M. Pierini,  _Anomaly Detection With Conditional Variational Autoencoders_ ,  In _Eighteenth International Conference on Machine Learning and Applications_ (2020), 2010.05531.
  * [42] M. van Beekveld, S. Caron, L. Hendriks, P. Jackson, A. Leinweber, S. Otten, R. Patrick, R. Ruiz De Austri, M. Santoni and M. White,  _Combining outlier analysis algorithms to identify new physics at the LHC_ (2020),  2010.07940.
  * [43] S. E. Park, D. Rankin, S.-M. Udrescu, M. Yunus and P. Harris,  _Quasi Anomalous Knowledge: Searching for new physics with embedded knowledge_ ,  JHEP 21, 030 (2020),  10.1007/JHEP06(2021)030,  2011.03550.
  * [44] P. Chakravarti, M. Kuusela, J. Lei and L. Wasserman,  _Model-Independent Detection of New Physics Signals Using Interpretable Semi-Supervised Classifier Tests_ (2021),  2102.07679.
  * [45] D. A. Faroughy,  _Uncovering hidden new physics patterns in collider events using Bayesian probabilistic models_ ,  PoS ICHEP2020, 238 (2021),  10.22323/1.390.0238,  2012.08579.
  * [46] T. Finke, M. Krämer, A. Morandini, A. Mück and I. Oleksiyuk,  _Autoencoders for unsupervised anomaly detection in high energy physics_ ,  JHEP 06, 161 (2021),  10.1007/JHEP06(2021)161,  2104.09051.
  * [47] O. Atkinson, A. Bhardwaj, C. Englert, V. S. Ngairangbam and M. Spannowsky,  _Anomaly detection with Convolutional Graph Neural Networks_ (2021),  2105.07988.
  * [48] B. M. Dillon, T. Plehn, C. Sauer and P. Sorrenson,  _Better Latent Spaces for Better Autoencoders_ (2021),  2104.08291.
  * [49] A. Kahn, J. Gonski, I. Ochoa, D. Williams and G. Brooijmans,  _Anomalous jet identification via sequence modeling_ ,  JINST 16(08), P08012 (2021),  10.1088/1748-0221/16/08/P08012,  2105.09274.
  * [50] T. Aarrestad _et al._ ,  _The Dark Machines Anomaly Score Challenge: Benchmark Data and Model Independent Event Classification for the Large Hadron Collider_ (2021),  2105.14027.
  * [51] S. Caron, L. Hendriks and R. Verheyen,  _Rare and Different: Anomaly Scores from a combination of likelihood and out-of-distribution models to detect new physics at the LHC_ (2021),  2106.10164.
  * [52] E. Govorkova, E. Puljak, T. Aarrestad, M. Pierini, K. A. Woźniak and J. Ngadiuba,  _LHC physics dataset for unsupervised New Physics detection at 40 MHz_ (2021),  2107.02157.
  * [53] J. Gonski, J. Lai, B. Nachman and I. Ochoa,  _High-dimensional Anomaly Detection with Radiative Return in $e^{+}e^{-}$ Collisions_ (2021),  2108.13451.
  * [54] E. Govorkova _et al._ ,  _Autoencoders on FPGAs for real-time, unsupervised new physics detection at 40 MHz at the Large Hadron Collider_ (2021),  2108.03986.
  * [55] P. T. Komiske, E. M. Metodiev and J. Thaler,  _Energy flow networks: deep sets for particle jets_ ,  Journal of High Energy Physics 2019(1), 121 (2019),  10.1007/JHEP01(2019)121,  1810.05165.
  * [56] H. Qu and L. Gouskos,  _ParticleNet: Jet Tagging via Particle Clouds_ ,  Phys. Rev. D 101(5), 056019 (2020),  10.1103/PhysRevD.101.056019,  1902.08570.
  * [57] B. M. Dillon, G. Kasieczka, H. Olischlager, T. Plehn, P. Sorrenson and L. Vogel,  _Symmetries, Safety, and Self-Supervision_ (2021),  2108.04253.
  * [58] Y. Zhang, J. Hare and A. Prügel-Bennett,  _Deep Set Prediction Networks_ ,  arXiv e-prints arXiv:1906.06565 (2019),  1906.06565.
  * [59] L. Ruff, R. Vandermeulen, N. Goernitz, L. Deecke, S. A. Siddiqui, A. Binder, E. Müller and M. Kloft,  _Deep one-class classification_ ,  In J. Dy and A. Krause, eds., _Proceedings of the 35th International Conference on Machine Learning_ , vol. 80 of _Proceedings of Machine Learning Research_ , pp. 4393–4402. PMLR (2018).
  * [60] D. community,  _Lhc simulation project_ ,  10.5281/zenodo.3685861 (2020).
  * [61] D. community,  _Unsupervised-hackathon_ ,  10.5281/zenodo.3961917 (2020).
  * [62] M. Zaheer, S. Kottur, S. Ravanbakhsh, B. Poczos, R. Salakhutdinov and A. Smola,  _Deep Sets_ ,  arXiv e-prints arXiv:1703.06114 (2017),  1703.06114.
  * [63] Y. Zhang, J. Hare and A. Prügel-Bennett,  _FSPool: Learning Set Representations with Featurewise Sort Pooling_ ,  arXiv e-prints arXiv:1906.02795 (2019),  1906.02795.
  * [64] P. T. Komiske, E. M. Metodiev and J. Thaler,  _Metric Space of Collider Events_ ,  Phys. Rev. Lett. 123(4), 041801 (2019),  10.1103/PhysRevLett.123.041801,  1902.02346.
  * [65] D. P. Kingma and M. Welling,  _Auto-Encoding Variational Bayes_ ,  arXiv e-prints arXiv:1312.6114 (2013),  1312.6114.
  * [66] D. Jimenez Rezende and S. Mohamed,  _Variational Inference with Normalizing Flows_ ,  arXiv e-prints arXiv:1505.05770 (2015),  1505.05770.
  * [67] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf _et al._ ,  _Pytorch: An imperative style, high-performance deep learning library_ ,  In H. Wallach, H. Larochelle, A. Beygelzimer, F. dAlché Buc, E. Fox and R. Garnett, eds., _Advances in Neural Information Processing Systems 32_ , pp. 8024–8035. Curran Associates, Inc. (2019).
  * [68] D. P. Kingma and J. Ba,  _Adam: A Method for Stochastic Optimization_ ,  arXiv e-prints arXiv:1412.6980 (2014),  1412.6980.
  * [69] B. Ostdiek,  _bostdiek/DarkMachines-UnsupervisedChallenge: arXiv_v1_ ,  10.5281/zenodo.4897467 (2021).


# Brief Notes on Hard Takeoff, Value Alignment, and Coherent Extrapolated
Volition

Gopal P. Sarma1
1\. _Emory University, Atlanta, GA USA_
Email<EMAIL_ADDRESS>
## I. On Hard Takeoff

The distinction between _hard takeoff_ and _soft takeoff_ has been used to
describe different possible scenarios following the arrival of human-level
artificial intelligence. The basic premise underlying these concepts is that
software-based agents would have the ability to improve their own intelligence
by analyzing and rewriting their source code, whereas biological organisms are
significantly more restricted in their capacity for self-improvement [1, 2, 3,
4, 5].

There is no precise boundary between the two scenarios, but in broad strokes,
a hard takeoff refers to a transition from human level intelligence to
superintelligence in a matter of minutes, hours or days. A soft takeoff refers
to a scenario where this transition is much more gradual, perhaps taking many
months or years. The practical importance of this qualitative distinction is
that in a soft takeoff, there may opportunities for human intervention in the
event that the initial AI systems have problematic design flaws.

The purpose of this brief note is simply to point out the following: takeoff
speed refers to the rate of change of the agents’ _level of intelligence_ ,
and not our perceived changes in the world around us. Because the notion of an
“intelligence explosion” has been constructed in analogy to a physical
explosion [6], it gives rise to an inaccurate mental picture in people’s
minds. If self-improving AI systems111As this analysis neither requires nor
implies that the driving force of change is a unitary agent, I have chosen to
use the plural terms “software-based agents,” “AI systems,” or
“superintelligent AI systems.” It may very well be a collection of agents /
systems possessing powerful AI capabilities in aggregate. are thought to be
the intellectual analogue of nuclear chain reactions, then the natural image
of an intelligence explosion that this metaphor creates is a scenario in which
massive, disruptive changes take place in the world that are difficult for
individuals and for society to handle.

However, the premise of intelligent agents with capacities in substantial
excess of any human being, which are able to process the sum total of human
knowledge in the form of books, video, and ongoing contemporary events implies
that greater levels of intelligence in the AI systems will be accompanied by
actions taken with a corresponding level of information, insight, and
operational skill. Therefore, if the initial systems are designed correctly
with respect to value alignment and goal structure stability, it is in fact a
hard takeoff scenario which would be _less disruptive_ than a soft takeoff,
not the other way around.

I reiterate this claim for emphasis: _The takeoff speed of an intelligence
explosion refers to the rate of change of intelligence in the AI systems, and
not our perceived changes in the world around us. Therefore, under the
assumption of correctly designed systems, a hard takeoff is preferable to a
soft takeoff because the resultant changes that take place in the world will
be executed with greater precision, thoughtfulness, and insight._

## II. On Value Alignment and Coherent Extrapolated Volition

The preceding argument relied on a key assumption, namely that the AI systems
capable of self-improvement were designed correctly with respect to value
alignment and goal structure stability. Value alignment refers to the
construction of systems that take actions consistent with human values.
Russell states 3 design principles which encapsulate the notion of value
alignment [7]:

  1. 1.

The machine’s purpose must be to maximize the realization of human values. In
particular, it has no purpose of its own and no innate desire to protect
itself.

  2. 2.

The machine must be initially uncertain about what those human values are. The
machine may learn more about human values as it goes along, but it may never
achieve complete certainty.

  3. 3.

The machine must be able to learn about human values by observing the choices
that we humans make.

A related notion is Yudkowsky’s “coherent extrapolated volition” [8, 9]. The
basic premise of this proposal is that sophisticated AI systems will be
capable of extrapolating and resolving disagreements between the value systems
of individuals and groups, ultimately arriving at a goal structure that
represents the collective desires of humanity. This process of iterated
reflection is analogous to Rawls’ “reflective equilibrium” [10].

Like the metaphor of hard takeoff, the notions of value alignment and coherent
extrapolated volition can also give rise to an inaccurate mental picture,
namely, that the aligned goal structure would either require or result in all
humans arriving at complete agreement on all issues. However, with adequate
resources, it may very well be that value-aligned AI systems shape a world in
which groups of individuals co-exist who disagree about object-level issues.
Certainly we can point to many examples in contemporary human society where
individuals maintain divergent preferences without conflict.

The purpose of this brief note is simply to point out the following: _Implicit
in any practical analysis of value alignment are the physical resources
available to the AI systems. In particular, the construction of a human
compatible goal structure does not mean that all human disagreements have been
resolved. Rather, it means that a mutually satisfactory set of outcomes has
been achieved, subject to resource constraints_.

It may be impossible to arrive at a consensus goal structure without adequate
resources. As a trivial example, if we have two individuals each of whom
desires at least one apple, there is no disagreement if we have two apples. On
the other hand, if there is only one apple, conflict may very well be
inevitable in the absence of other factors with which to resolve the imbalance
between specific individual desires and total available resources. In the
context of superintelligent AI systems capable of exerting substantial
influence on the world and shaping society on a global scale, adequate
analysis of value alignment requires an understanding of the sum total of
physical resources that the AI systems have at their disposal.

### Acknowledgements

I would like to thank Eric Drexler and Anders Sandberg for valuable
discussions and feedback on the manuscript.

### ORCID

Gopal P. Sarma  0000-0002-9413-6202

## References

  * [1] N. Bostrom, Superintelligence: Paths, Dangers, Strategies.  OUP Oxford, 2014.
  * [2] I. J. Good, “Speculations Concerning the First Ultraintelligent Machine,” Advances In Computers, vol. 6, no. 99, pp. 31–83, 1965.
  * [3] D. Chalmers, “The Singularity: A Philosophical Analysis,” Journal of Consciousness Studies, vol. 17, no. 9-10, pp. 7–65, 2010.
  * [4] M. Shanahan, The Technological Singularity.  MIT Press, 2015.
  * [5] C. Shulman and A. Sandberg, “Implications of a software-limited singularity,” in Proceedings of the European Conference of Computing and Philosophy, 2010\.
  * [6] E. Yudkowsky, “Artificial intelligence as a positive and negative factor in global risks,” in Global Catastrophic Risks (N. Bostrom and M. Cirkovic, eds.), ch. 15, pp. 308–345, Oxford: Oxford University Press, 2011\.
  * [7] S. Russell, “Should We Fear Supersmart Robots?,” Scientific American, vol. 314, no. 6, pp. 58–59, 2016.
  * [8] E. Yudkowsky, “Coherent Extrapolated Volition,” Singularity Institute for Artificial Intelligence, 2004.
  * [9] N. Tarleton, “Coherent extrapolated volition: A meta-level approach to machine ethics,” Machine Intelligence Research Institute, 2010.
  * [10] J. Rawls, A Theory of Justice.  Belknapp, 1971.


# STR-GODEs: Spatial–Temporal-Ridership Graph ODEs for Metro Ridership
Prediction

Chuyu Huang
Peking University
<EMAIL_ADDRESS>

###### Abstract

The metro ridership prediction has always received extensive attention from
governments and researchers. Recent works focus on designing complicated graph
convolutional recurrent network architectures to capture spatial and temporal
patterns. These works extract the information of spatial dimension well, but
the limitation of temporal dimension still exists. We extended Neural ODE
algorithms to the graph network and proposed the STR-GODEs network, which can
effectively learn spatial, temporal, and ridership correlations without the
limitation of dividing data into equal-sized intervals on the timeline. While
learning the spatial relations and the temporal correlations, we modify the
GODE-RNN cell to obtain the ridership feature and hidden states. Ridership
information and its hidden states are added to the GODESolve to reduce the
error accumulation caused by long time series in prediction. Extensive
experiments on two large-scale datasets demonstrate the efficacy and
robustness of our model. 111Our code and benchmarks are available at
https://github.com/universe-hcy/STR-GODEs.

##  1 Introduction

With the increasing population of cities, metro lines have become more and
more crowded. For improving the service efficiency of the metro system and
facilitating passenger travel, metro ridership prediction has become hotspot
research in the community of Intelligent Transportation Systems (ITSs). Recent
works focus on designing complicated graph convolutional recurrent network
architectures to capture spatial and temporal patterns. These works depend on
dividing the ridership data into equal-sized intervals on the timeline and
counting the ridership as the time period. In metro ridership data, there are
distinct characteristics of peak and valley periods. The dense ridership
during the peak period needs to be divided more finely, while the ridership in
the low valley period is less, and there is not enough data when the time
interval is small. This preprocessing destroys information about the ridership
over time, which may be information about potential variables. At the same
time, the problem is limited to time period and time period prediction, which
is not conducive to real-time data collection and feedback to passenger users.

A more practical approach is to construct a continuous-time model that defines
a potential state at all times. Chen et al.[5] proposed a new family of
continuous time-series models (NODE) which directly models the dynamics of
network hidden states with an ODE solve. Rubanova et al.[29] generalize state
transitions in RNNs to continuous-time dynamics specified by Neural ODEs.and
proposes latent ODE network defined in continuous time, which makes it more
effective in irregular sampling data. Zhou et al.[45, 46] Applied the NODE to
the high-resolution reconstruction and prediction of urban flow grid map.
These work provide a new direction for metro passenger flow prediction.

Previous work on NODE mainly focused on Euclidean data such as sequences and
matrices. We modified it to make it suitable for non-Euclidean data, such as
graph networks, called Graph ODEs(GODE).To facilitate performance comparison
between models, we used the three complementary graphs proposed by PVCGN[22]
as spatial information to learn the relationship between subway stations.
Inspired by latent ODEs, we use a GODE-RNN structure as encoder to learn the
initial state $Z_{0}$, and hidden states containing the characteristics of
ridership. Subsequently, we modify the GODEs network in prediction so that the
ridership information and its hidden vectors can be incorporated into the
prediction process, which improves the cumulative error caused by the long
time series. To verify the effectiveness of our STR-GODEs, we conduct
experiments on two large-scale benchmarks(i.e., Shanghai Metro and Hangzhou
Metro) with widely used baselines and state-of-the-art models and the
evaluation results show that our approach outperforms existing state-of-the-
art methods under various comparison circumstances.

In summary, our contribution is threefold:

  * •

We further extend the neural ODEs network to the time series prediction of
graph networks, and to the best of our knowledge, firstly apply it to the
metro ridership prediction.

  * •

By combining the ridership data and its hidden state, we improve the
performance of the prediction model of neural grouph ODEs network so that it
can reduce the cumulative error over a long time series.

  * •

our STR-GODEs has better effect and robustness, and is no longer limited by
the same time slice and equal distance prediction in previous work. It can be
better applied in real life. Since we have only done some early work, this
direction is very prospective.

##  2 Background and Related Works

Traffic States Prediction Traffic forecasting depends on the combination of
spatial-temporal feature and has been studied for decades. Traditional
methods[24, 27, 34, 31] can only consider the relation of a single station in
temporal dimension, lack spatial information, and can only learn simple time
series models.In recent years, deep neural network has become the mainstream
method in this field. The early work [40, 39, 23] mainly divided the studied
cities into regular grid maps, and transformed the raw traffic data into
tensors. CNN is used to capture the spatial correlation among near regions.
This preprocessing method will bring structural noise information into the
traffic system with irregular topology, which will affect the effect of the
model. Graph convolution network[8, 19] greatly improves the disadvantages of
previous work on spatial information. Recent works[44, 11, 2] focus on
designing complicated graph convolutional recurrent network architectures to
capture the spatial and temporal patterns. DCRNN [21] re-formulates the
spatial dependency of traffic as a diffusion process and extends the previous
GCN to a directed graph. Graph-WaveNet[36] captures the hidden spatial
dependency in the data by constructing a adaptive dependency matrix and
learning it through node embedding. Liu et al.[22] construct three
complementary graphs to utilize the metro physical topology information. A
Graph Convolution Gated Recurrent Unit (GC-GRU) is applied to learn the
spatial-temporal representation and a Fully-Connected Gated Recurrent Unit
(FC-GRU) is applied to capture the global evolution tendency. However, due to
the limitation of RNN structure, there are still some deficiencies in the
temporal dimension.

Graph Neural Network Unlike traditional Euclidean data, graph data is
difficult to be embedded into Euclidean space losslessly. Graph Convolution
Networks (GCN)[8, 19] have been proposed to automatically learn feature
representation on graphs,which is widely used in node classification[19], link
prediction[43], and graph classification[42].[35] There are two mainstreams of
graph convolution networks, the spectralbased approaches [4, 8, 19]and the
spatial-based approaches[1, 13, 14]. In these approaches, the adjacency matrix
is considered as prior knowledge and is fixed throughout training. Recently
researchers[20, 36, 2] pay more attention on operating on both spatial and
temporal dimensions without the pre-defined graph structure.

Neural Ordinary Differential Equations Built upon previous works[33, 25], Chen
et al.[5] proposed a new family of continuous-time models named Neural
ODEs,which can be interpreted as a continuous equivalent of Residual
networks[15]. A basic formulation of Neural ODEs is shown as:

$\displaystyle\frac{dh(t)}{dt}=f_{\theta}(h(t),t),\quad h(t_{0})=h_{0}$ (1)

where $f_{\theta}$ is parametrized as a neural network.The hidden state h(t)
is defined at all times, and can be evaluated at any desired times using a
numerical black-box ODE solver:

$\displaystyle
h_{0},\ldots,h_{N}=ODESolve(f_{\theta},h_{0},(t_{0},\ldots,t_{N}))$ (2)

By using adjoint method[5],the gradients w.r.t.$\theta$ can be computed
memory-efficiently during back propagation. Recent works[10, 26, 7, 38, 12]
have tried to analyze this framework theoretically, overcome the instability
issue and improve memory-efficiency. Researchers also pay attention to apply
NODE to other fields such as medical image[28], reinforcement learning[9],
video generation[17, 41] and graph data[37].

##  3 Method

###  3.1 Problem Definition

Previous work did not support predictions over continuous or unequal density
time intervals. To facilitate performance comparison between models, we follow
the traditional prediction framework. Assuming N is the number of Metro
stations, the time period T data can be expressed as
$X_{t}=(X_{t}^{1},X_{t}^{2},...,X_{t}^{N})$, where
$X_{t}^{i}\in\mathbb{R}^{2}$ represents the the passenger counts of
inflow/outflow of the station i at time interval t. our target is to predict
the future ridership sequence based on the observed historical values:

$\displaystyle\hat{X}_{t+1},\hat{X}_{t+2},...,\hat{X}_{t+m}=STR-
GODEs(X_{t-n+1},X_{t-n+2},...,X_{t})$ (3)

###  3.2 STR-GODEs

Suppose $Z(t)$ is the latent state including the ridership feature at time
interval t, and $f(t,Z(t),\theta)$ is the function of ridership changing with
spatial,temporal and ridership information. By Euler method, there is
$Z(t+\Delta t)=Z(t)+\Delta t*f(t,Z(t),\theta)$.

The primal problem can be abstracted as we learn the ridership trend function
$f(t,Z(t),\theta)$ according to the training data and pre-defined graph
(Physical-Virtual Graph[22]), obtain the initial state $Z_{0}$ through the
observation sequence, and predict the subsequent sequence.

Chen et al.[5] present a continuous-time, generative approach to modeling time
series. Each trajectory is determined from a local initial state, $Z_{t_{0}}$
, and a global set of latent dynamics shared across all time series:

$\displaystyle Z_{t_{0}}\sim{p(Z_{t_{0}})}$ (4) $\displaystyle
Z_{t_{0}},Z_{t_{1}},\ldots,Z_{t_{n}}=ODESolve(Z_{t_{0}},f,\theta_{f},t_{0},\ldots,t_{N})$
(5) $\displaystyle each\quad x_{t_{i}}\sim{p(x|Z_{t_{i}},\theta_{x})}$ (6)

We follow the framework of Chen et al.[5] for both training and prediction. We
modify the original ODE Block to make it suitable for graph network data,
called GODE Block. As shown in figure 1, inspired by Rubanova et al.[29], we
use the GODE-RNN structure to obtain the local initial state $Z_{t_{0}}$ by
learning the spatial,temporal and ridership feature of
$X_{t},X_{t+1},...,X_{t+n}$. In particular, due to the need to gather
information from adjacent metro stations in GODE Block, we use GODE-RNN to
learn the approximate posterior of local initial state $Z_{t_{0}}$ directly
instead of learning the mean and standard deviation.

$\displaystyle q(z_{t_{0}}|\\{x_{i},t_{i}\\}^{N}_{i=0})=Z_{t_{0}}\quad
whereZ_{t_{0}},hidden_{t_{0}}=g(ODE-RNN_{\phi}(x_{i},t_{i})^{N}_{i=0})$ (7)

Furthermore, by combining the ridership information and its hidden state, we
further improve the performance of the prediction model so that it can reduce
the cumulative error over a long time series.

Figure 1: The illustration of encoder and decoder used in STR-GODEs.In
GODESolve, the green line represents ground truth, and the pink line
represents the predicted result obtained by GODE. At some special time points,
we correct the predicted result by introducing ridership information and its
hidden state through GRU cell.

###  3.3 Physical-Virtual Graphs

In this section, we follows the pre-defined graph proposed by PVCGN[22] to
learn spatial information. In our work, the physical graph, similarity graph
and correlation graph are denoted as
$\mathcal{G}_{p}=(\mathcal{V},\mathcal{E}_{p},W_{p})$,
$\mathcal{G}_{s}=(\mathcal{V},\mathcal{E}_{s},W_{s})$ and
$\mathcal{G}_{c}=(\mathcal{V},\mathcal{E}_{c},W_{c})$, where $\mathcal{V}$ is
a set of nodes represent real-world metro stations($|\mathcal{V}|=N$),
$\mathcal{E}$ is a set of edges and $W\in\mathcal{R}^{N\times N}$ denotes the
weights of all edges.

Physical Graph:$\mathcal{G}_{p}$is directly built based on the physical
topology of the studied metro system. We first construct a physical connection
matrix $P\in\mathbb{R}^{N\times N}$, where $P(i,j)=1$ if there exists an edge
between node i and j, or else $P(i,j)=0$. Note that each diagonal value P(i,
i) is directly set to 0.Finally, the edge weight $W_{p}$ is obtained by
performing a linear normalization on each row:

$\displaystyle W_{p}(i,j)=\frac{P(i,j)}{\sum_{k=1}^{N}P(i,k)}$ (8)

Similarity Graph:The similarities of metro stations are used to guide the
construction of $\mathcal{G}_{s}$. We construct a similarity score matrix
$S\in\mathbb{R}^{N\times N}$, where the score S(i, j) between station i and j
is computed with Dynamic Time Warping[3]: $S(i,j)=exp(-DWT(X^{i},X^{j}))$.
Based on the matrix S, we select some station pairs to build edges
$\mathcal{E}_{s}$ with a predefined similarity threshold 0.1 for HZMetro and
by choosing the top ten stations with high similarity scores for SHMetro.
Finally, we calculate the edge weights $W_{s}$ by conducting row normalization
on S:

$\displaystyle W_{s}(i,j)=\frac{S(i,j)}{\sum_{k=1}^{N}S(i,k)\cdot
L(E_{s},i,k)}$ (9)

where $L(\mathcal{E}_{s},i,k)=1$ if $\mathcal{E}_{s}$ contains an edge
connecting node i and k, or else $L(\mathcal{E}_{s},i,k)=0$.

Correlation Graph: We utilize the origin-destination distribution of ridership
to build the virtual graph $\mathcal{G}_{c}$. First, we construct a
correlation ratio matrix $C\in\mathbb{R}^{N\times N}$:
$C(i,j)=\frac{D(i,j)}{\sum_{k=1}^{N}D(i,k)}$ where D(i, j) is the total number
of passengers that traveled from station j to station i in the whole training
set. Based on the matrix C, we select some station pairs to build edges
$\mathcal{E}_{c}$ with a predefined similarity threshold 0.02 for HZMetro and
by choosing the top ten stations with high similarity scores for SHMetro.
Then, the edge weights $W_{c}$ is calculated by:

$\displaystyle W_{c}(i,j)=\frac{C(i,j)}{\sum_{k=1}^{N}C(i,k)\cdot
L(E_{c},i,k)}$ (10)

###  3.4 GODE Block

In general, the graph convolutional networks is a method to aggregate the
information of neighbor nodes around the target node and can be expressed
as[8, 19, 14, 32]:

$\displaystyle
x_{i}^{k}=\gamma(x_{i}^{k-1},\square_{j\in\aleph(i)}\phi^{k}(x_{i}^{k-1},x_{j}^{k-1},e_{i,j}))$
(11)

where $x_{i}^{k}$ refers to the hidden state of node i in the k-th layer,
$e_{i,j}$ refers to the eigenvector of the edge from node i to node j,
$\square$ is the aggregation method of neighbor information, $\phi$ is the
transformation method of neighbor information and $\gamma$ is a transformation
that fuses one’s own feature with its aggregated neighbors’information.

These and related methods can be understood as special cases of a simple
differentiable message-passing framework (Gilmer et al. 2017[13]):

$\displaystyle
h_{i}^{l+1}=\sigma\left(\sum\limits_{m\in\mathcal{M_{i}}}g_{m}(h_{i}^{l},h_{j}^{l})\right)$
(12)

where $h_{i}^{l+1}\in\mathcal{R}^{d}$ is hidden state of node $v_{i}$ in the
l-th layer of the neural network. Incoming messages of the form
$g_{m}(\cdot,\cdot)$ are accumulated and passed through an element-wise
activation function$\sigma(\cdot)$. $\mathcal{M}_{i}$denotes the set of
incoming messages for node $v_{i}$ and is often chosen to be identical to the
set of incoming edges. It can also be expressed in the following form:

$\displaystyle h_{t+1}=h_{t}+f(h_{t},\theta_{t},\mathcal{G})$ (13)

According to Neural ODEs[5], the hidden state h(t) can be defined as the
solution to an ODE initial-value problem:

$\displaystyle\frac{dh(t)}{dt}=f_{\theta}(h(t),t,\mathcal{G}),where\quad
h(t_{0})=h_{0}$ (14)

Inspired by R-GCN[30], we define our GODE-function by using learnable weights
$\mathbf{\Theta}_{0}$ to learn the changes of stations latent states over
time, and learnable weights $\mathbf{\Theta}_{r}$ to learn the association
relationship between stations. For every station,there are:

$\displaystyle\mathbf{z}^{t+\Delta
t}_{i}=relu(\mathbf{\Theta}_{0}\cdot\mathbf{z}^{t}_{i}+\sum_{r\in\mathcal{R}}\sum_{j\in\mathcal{N}_{r}(i)}\frac{1}{|\mathcal{N}_{r}(i)|}\mathbf{\Theta}_{r}\cdot\mathbf{z}^{t}_{j})$
(15)

For any time $t_{i}$, we can calculate the latent representation $Z_{t_{i}}$
by:

$\displaystyle
Z_{t_{i}}=Z_{t_{i-1}}+\int_{t_{i-1}}^{t_{i}}{f_{\theta}(h(t),t,\mathcal{G})}\,{\rm
d}t$ (16)

Input: Data points and their timestamps$\\{(x_{i},t_{i})\\}_{i=1...N}$

$h_{N},hidden=0,0$

for i in N,N-1,…,1 do

$time\\_points=utils.linspace\\_vector(t_{i-1},t_{i},n\\_intermediate\\_tp)$

$h^{\prime}={\color[rgb]{0,0,1}GODESolve}(f_{\theta},h_{i-1},time\\_points)$

$h_{i},{\color[rgb]{0,0,1}hidden}=RNNCell(h^{\prime},{\color[rgb]{0,0,1}hidden},x_{i})$

end for

Return $h_{0},hidden$

Algorithm 1 The GODE-RNN. The difference is highlighted in blue

###  3.5 STR-GODEs encoder

Rubanova et al. proposed using an ODE-RNN as the encoder for a latent ODEs
model to learn the mean $\mu_{z_{0}}$ and standard deviation $\sigma_{z_{0}}$
of the approximate posterior $q(z_{t_{0}}|\\{x_{i},t_{i}\\}^{N}_{i=0})$ and
the local initial state $Z_{t_{0}}=\mathcal{N}(\mu_{z_{0}},\sigma_{z_{0}})$.

As shown in algorithm 1, we made some modifications on ODE-RNN. We replaced
ODESolve with GODESolve, which can be applied to graph data. We use GODE-RNN
to learn the approximate posterior of local initial state $Z_{t_{0}}$ directly
instead of learning the mean and standard deviation. The hidden state of
ridership information is learned and stored in hidden, which is combined with
$Z_{i}$ in RNNCell. In our work, we use the classical GRU[6] structure as the
RNNCell. Specifically, the reset gate
$r_{t}=\\{r_{t}^{1},r_{t}^{2},\ldots,r_{t}^{N}\\}$, update gate
$z_{t}=\\{z_{t}^{1},z_{t}^{2},\ldots,z_{t}^{N}\\}$, new information
$N_{t}=\\{N_{t}^{1},N_{t}^{2},\ldots,N_{t}^{N}\\}$ are computed by:

$\displaystyle r_{t}=\sigma(W_{r}\cdot[Z_{t-1},h_{t-1},x_{t}])$ (17)
$\displaystyle z_{t}=\sigma(W_{z}\cdot[Z_{t-1},h_{t-1},x_{t}])$ (18)
$\displaystyle N_{t}=tanh(W_{N}\cdot[r_{i}*h_{t-1},x_{t}])$ (19)
$\displaystyle h_{t}=(1-z_{t})*N_{t-1}+z_{t}*h_{t-1}$ (20) $\displaystyle
Z_{t}=(1-z_{t})*N_{t-1}+z_{t}*Z_{t-1}$ (21)

Input: prediction timestamps$\\{t_{i}\\}_{i=1...M},Z_{0}$,hidden,Data
points$\\{x_{i}\\}_{i=1...N}$

for i in 1,2,…,M do

$Z^{\prime}=odeint(f_{\theta},Z_{i-1},(t_{i-1},t_{i}))$

$y_{i}=output\\_layer(Z^{\prime})$

$x_{i}=x_{i}$ if $x_{i}$ is available Data points else $y_{i}$

$Z_{i},hidden=RNNCell(Z^{\prime},hidden,x_{i})$

end for

Return $\\{y_{i}\\}_{i=1...M}$

Algorithm 2 The GODESolver_prediction.

###  3.6 STR-GODEs decoder

Before prediction, we convert the latent state $h_{0},hidden$ obtained from
GODE-RNN to get $Z_{0},hidden$ with a transform block which is composed of a
fully-connected layer, a tanh activation function, and a fully-connected
layer.

Previous work uses Neural ODEs as decoder for prediction and when the sequence
is long, the error will accumulate gradually. Similar to the GODE-RNN, We
introduce ridership information and its hidden state into GODESolver through
RNNCell structure, as shown in algorithm 2. The experimental results show that
it can reduce the cumulative error over a long time series. Specifically, we
use the classical GRU[6] structure described in session 3.5 as the RNNCell.

##  4 Experiment

###  4.1 Datasets

To evaluate the performance of our work, we conduct experiments on two real-
world traffic: SHMetro and
HZMetro222https://tianchi.aliyun.com/competition/entrance/231708/information.

SHMetro: The SHMetro dataset refers to the ridership data on the metro system
of Shanghai, China. There are 288 metro stations connected by 958 physical
edges. A total of 811.8 million transaction records were collected from Jul.
1st 2016 to Sept. 30th 2016, with 8.82 million ridership per day. Following
the setting of PVCGN[22], for each station, we measured its inflow and outflow
every 15 minutes by counting the number of passengers entering or exiting the
station. The ridership data of the first two months and that of the last three
weeks are used for training and testing, while the ridership data of the
remaining days are used for validation.

HZMetro: The HZMetro dataset refers to the ridership data on the metro system
of Hangzhou, China. There are 80 metro stations connected by 248 physical
edges with 2.35 million ridership per day. The time interval of this dataset
is also set to 15 minutes. Similar to SHMetro, this dataset is divided into
three parts, including a training set (Jan. 1st - Jan. 18th), a validation set
(Jan. 19th - Jan. 20th), and a testing set (Jan. 21th - Jan. 25th).

###  4.2 Experimental Settings

We have trained the official code of PVCGN[22], which has a small gap with the
results of corresponding papers. To reduce the effect difference caused by
different parameters and facilitate the model performance comparison, we keep
the parameters in STR-GODEs the same as those in the official code of PVCGN.

We implement our STR-GODEs in Python with Pytorch 1.6.0 and executed on a
server with one NVIDIA Titan X GPU card. The lengths of input and output
sequences are set to 4 simultaneously. The input data and the ground-truth of
output are normalized with Z-score
Normalization333https://en.wikipedia.org/wiki/Standard_score before being fed
into the network. The batch size is set to 8 and the feature dimensionality d
is set to 256. The initial learning rate is 0.001 and its decay ratio is 0.1.
We optimize the models by Adam[18] optimizer for a maximum of 200 epochs by
minimizing the mean absolute error between the predicted results and the
corresponding ground-truths. The best parameters for all deep learning models
are chosen through a carefully parameter-tuning process on the validation set.

###  4.3 Compared Methods

To evaluate the overall performance of our work, we compare STR-GODEs with
widely used baselines and state-of-the-art models:

DCRNN: Huang et al.[16] proposed a deep learning framework specially designed
for traffic prediction, which uses bidirectional random walks on graphs to
capture spatial dependencies and an encoder-decoder architecture to learn
temporal correlation. We implement this method with its official code.

Graph-WaveNet: Wu et al.[36] developed an adaptive dependency matrix to
capture the hidden spatial dependencies, and used a stacked dilated 1D
convolution component to handle very long sequences. We implement this method
with its official code.

ST—ODEs: Zhou et al. applied the NODE to the prediction of urban flow grid
maps. They concats the data over different time intervals, feeds it into CNN
layer to get the initial state $Z_{0}$, and then uses the neural ODEs network
for prediction.In this experiment, we reproduced the code and extended the ODE
network to GODE network which can be used for graph data.

latent ODEs: Rubanova et al. refine the Latent ODE model of Chen et al. [5] by
using the ODE-RNN as a recognition network, where ODE-RNN is generalized from
RNNs with continuous-time hidden dynamics defined by ordinary differential. In
the experiment, we modify its official code and extended the ODE network to
GODE network which can be used for graph data.

PVCGN: Liu et al.[22] developed a Seq2Seq model with GC-GRU and FC-GRU to
forecast the future metro ridership sequentially where Graph Convolution Gated
Recurrent Unit (GC-GRU) is applied to learn the spatial-temporal
representation and Fully-Connected Gated Recurrent Unit (FC-GRU) is applied to
capture the global evolution tendency. We implement this method with its
official code.

Time | Metric | DCRNN | Graph-wavenet | ST—ODEs | latent ODEs | PVCGN | STR-GODEs
---|---|---|---|---|---|---|---
| MAE | 23.81 | 24.01 | 24.54 | 22.79 | 22.70 | 22.87
15min | RMSE | 40.52 | 40.68 | 42.43 | 38.61 | 37.78 | 37.71
| MAPE | 14.13% | 14.31% | 15.01% |  13.54% |  13.54% | 13.78 %
| MAE | 25.20 | 25.45 | 25.34 | 23.53 | 23.71 | 23.24
30min | RMSE | 42.51 | 42.02 | 43.49 | 39.63 | 39.96 | 38.69
| MAPE | 14.89% | 15.18% | 15.43% | 13.85% | 14.06% |  14.03 %
| MAE | 26.86 | 27.21 | 26.56 | 24.44 | 24.53 | 23.79
45min | RMSE | 45.93 | 46.03 | 46.00 | 41.22 | 41.31 | 39.68
| MAPE | 16.02% | 17.39% | 16.50% | 14.67% |  14.66% |  14.66 %
| MAE | 28.41 | 29.21 | 27.45 | 24.99 | 24.88 | 24.35
60min | RMSE | 49.23 | 49.77 | 47.61 | 42.20 | 41.61 | 40.83
| MAPE | 18.06% | 19.41% | 19.01% | 16.05% | 15.68% |  15.56 %
Table 1: results of conventional prediction experiment on HZMetro Time | Metric | DCRNN | Graph-wavenet | ST—ODEs | latent ODEs | PVCGN | STR-GODEs
---|---|---|---|---|---|---|---
| MAE | 24.09 | 24.86 | 24.79 | 23.51 | 23.41 | 23.21
15min | RMSE | 46.11 | 46.82 | 46.96 | 44.64 | 45.31 | 44.58
| MAPE | 17.87% | 19.98% | 18.56% | 17.64% | 17.04% |  16.99 %
| MAE | 25.23 | 26.73 | 25.77 | 25.06 | 24.36 | 23.63
30min | RMSE | 49.92 | 52.01 | 50.08 | 48.91 | 48.12 | 46.28
| MAPE | 18.32% | 20.42% | 19.29% | 18.37% | 17.35% |  17.12 %
| MAE | 26.84 | 28.71 | 26.80 | 26.18 | 25.29 | 24.65
45min | RMSE | 55.07 | 58.32 | 53.59 | 52.60 | 51.42 | 49.93
| MAPE | 19.37% | 21.93% | 20.59% | 19.22% | 17.84% |  17.58 %
| MAE | 28.22 | 30.61 | 27.63 | 26.83 | 26.21 | 25.56
60min | RMSE | 60.18 | 64.02 | 56.04 | 55.78 | 55.17 | 53.39
| MAPE | 20.57% | 24.22% | 22.38% | 20.28% | 18.83% |  18.25 %
Table 2: results of conventional prediction experiment on SHMetro

###  4.4 Performance Comparison and Analysis

We measure the performance of predictive models with three widely used metrics
- Mean Absolute Error (MAE), Root Mean Square Error (RMSE), and Mean Absolute
Percentage Error (MAPE).

Conventional prediction experiment: Table 1 and Table 2 present the prediction
performances in conventional prediction experiments of our STR-GODEs and
representative comparison methods in HZMetro and SHMetro datasets. Compared
with latent ODE, We can observe that the GODE network combined with ridership
data and hidden states can indeed reduce the cumulative error. Compared with
all other methods, STR-GODEs has achieved some improvement in most metrics.

Further more, we extract the prediction of ridership in peak periods(7:30-9:30
and 17:30-19:30) in conventional experiments, as shown in table 6 and table 7.
We can observe that our STR-GODEs outperforms all comparative methods
consistently on most metrics on both datasets. On HZMetro, we achieve a
relative improvement of 4.13% in MAE metric compared with other optimal
algorithms on average, 3.59% in RMSE and 4.28% MAPE. On SHMetro, we achieve a
relative improvement of 4.32% in MAE metric compared with other optimal
algorithms on average, 4.57% in RMSE and 1.84% MAPE.

Time | Metric | DCRNN | Graph-wavenet | ST—ODEs | latent ODEs | PVCGN | STR-GODEs
---|---|---|---|---|---|---|---
| MAE | 34.22 | 36.42 | 34.89 | 33.17 | 32.68 | 31.68
15min | RMSE | 53.09 | 56.13 | 54.71 | 52.20 | 50.23 | 48.78
| MAPE | 10.02% | 10.41% | 10.02% | 9.65% | 9.67% |  9.31 %
| MAE | 36.79 | 38.47 | 36.37 | 34.64 | 33.83 | 32.37
30min | RMSE | 57.13 | 59.05 | 56.84 | 53.82 | 52.59 | 50.33
| MAPE | 10.27% | 10.74% | 10.13% | 9.89% | 9.58% |  9.18 %
| MAE | 37.53 | 38.71 | 36.15 | 34.22 | 33.85 | 32.06
45min | RMSE | 59.89 | 59.92 | 57.11 | 53.89 | 53.28 | 50.69
| MAPE | 10.71% | 11.63% | 10.78% | 10.44% | 10.04% |  9.62 %
| MAE | 37.16 | 37.89 | 34.85 | 32.71 | 32.05 | 30.82
60min | RMSE | 58.93 | 60.13 | 54.98 | 52.36 | 51.64 | 50.45
| MAPE | 11.13% | 12.04% | 11.86% | 11.20% | 10.74% |  10.20 %
Table 3: results of conventional prediction experiment in peak periods on HZMetro. The peak periods refer to 7:30-9:30 and 17:30-19:30 Time | Metric | DCRNN | Graph-wavenet | ST—ODEs | latent ODEs | PVCGN | STR-GODEs
---|---|---|---|---|---|---|---
| MAE | 38.12 | 38.84 | 38.32 | 35.31 | 36.57 | 35.17
15min | RMSE | 67.74 | 68.01 | 66.89 | 61.87 | 64.87 | 61.93
| MAPE | 14.07% | 14.23% | 13.99% | 13.27% | 13.32% |  13.03 %
| MAE | 39.72 | 42.73 | 40.46 | 38.12 | 37.93 | 35.81
30min | RMSE | 73.14 | 77.48 | 72.20 | 69.02 | 68.91 | 64.34
| MAPE | 14.19% | 14.92% | 14.57% | 13.94% | 13.49% |  13.20 %
| MAE | 41.52 | 45.33 | 41.71 | 38.44 | 38.51 | 36.21
45min | RMSE | 79.03 | 85.14 | 77.12 | 71.12 | 73.55 | 67.16
| MAPE | 15.13% | 15.79% | 15.53% | 14.67% | 14.11% |  13.78 %
| MAE | 41.13 | 45.74 | 41.10 | 37.51 | 37.87 | 35.44
60min | RMSE | 78.81 | 89.21 | 78.73 | 71.43 | 73.13 | 66.99
| MAPE | 16.12% | 17.21% | 16.92% | 15.89% | 15.03% |  14.87 %
Table 4: results of conventional prediction experiment in peak periods on
SHMetro. The peak and valley periods refer to 7:30-9:30 and 17:30-19:30

Irregular prediction experiment: In order to compare the performance of
different methods in continuous time or unequal interval time, we use a
compromise experiment that the traditional method can also work.

In this experiment, assuming that $X_{t_{1}},X_{t_{2}},\ldots,X_{t_{N}}$ are
the observed ridership sequence, where $t_{1},t_{2},\ldots,t_{N}$ is the
random subsequence of $t-n+1,t-n+2,\ldots,t+1,\ldots,t+m$, our goal is to
predict a the ridership sequence:

$\displaystyle\hat{X}_{t_{n+1}},\hat{X}_{t_{n+2}},...,\hat{X}_{t_{n+m}}=STR-
GODEs(X_{t_{1}},X_{t_{2}},...,X_{t_{n}})$ (22)

Time | Metric | ST—ODEs | latent ODEs | PVCGN | STR-GODEs
---|---|---|---|---|---
| MAE | 31.53 | 19.95 | 22.19 | 17.70
$t_{n+1}$ | RMSE | 63.46 | 40.12 | 38.91 | 37.31
| MAPE | 16.34% | 11.42% | 13.12% |  10.14%
| MAE | 29.70 | 20.83 | 23.57 | 18.25
$t_{n+2}$ | RMSE | 56.69 | 40.21 | 41.22 | 36.04
| MAPE | 15.90% | 12.34% | 13.98% |  10.75%
| MAE | 30.04 | 21.10 | 22.22 | 18.19
$t_{n+3}$ | RMSE | 57.51 | 42.67 | 40.32 | 37.52
| MAPE | 15.82% | 11.92% | 13.21% |  10.08%
| MAE | 33.24 | 22.02 | 20.63 | 18.45
$t_{n+4}$ | RMSE | 68.08 | 44.71 | 38.18 | 37.41
| MAPE | 18.54% | 13.13% | 13.08% |  10.74%
Table 5: results of irregular prediction experiment on HZMetro Time | Metric | ST—ODEs | latent ODEs | PVCGN | STR-GODEs
---|---|---|---|---|---
| MAE | 27.32 | 16.61 | 25.18 | 16.11
$t_{n+1}$ | RMSE | 59.35 | 38.46 | 51.60 | 38.66
| MAPE | 19.32% | 12.55% | 17.61% |  11.94%
| MAE | 27.55 | 16.36 | 25.36 | 15.69
$t_{n+2}$ | RMSE | 59.45 | 38.39 | 52.73 | 37.36
| MAPE | 19.99% | 12.66% | 18.00% |  11.96%
| MAE | 27.94 | 15.63 | 25.01 | 14.96
$t_{n+3}$ | RMSE | 61.24 | 37.24 | 53.52 | 36.33
| MAPE | 21.13% | 12.21% | 18.21% |  11.19%
| MAE | 28.93 | 16.12 | 22.29 | 15.66
$t_{n+4}$ | RMSE | 61.88 | 39.88 | 50.83 | 39.72
| MAPE | 24.08% | 12.24% | 16.81% |  11.02%
Table 6: results of irregular prediction experiment on SHMetro

Table 3 and Table 4 present the prediction performances in irregular
prediction experiments of our STR-GODEs and representative comparison methods
in HZMetro and SHMetro datasets. Through the experimental results, we can see
that Neural GODEs network is superior to the traditional method in irregular
prediction and compared with other work, the superiority in effectiveness of
the proposed STR-GODEs is verified.

##  5 Conclusion

In this paper, we extended Neural ODE algorithms to the graph network and
proposed the STR-GODEs network, which can effectively learn spatial, temporal,
and ridership correlations without dividing data into equal-sized intervals on
the timeline. Specifically, we use a GODE-RNN structure as encoder to learn
the initial state $Z_{0}$, and hidden states containing the characteristics of
ridership. Subsequently, GODE network which can incorporate the ridership
information and its hidden vectors is used to predict the future sequence.
Experimental results on two public large-scale datasets demonstrate the
superior performance of our algorithm.

In future works, we will make further use of the continuous time series-based
characteristics of the model to meet the requirements of real-time prediction
and apply the model to the passenger. At the same time, we will further
optimize the time and space efficiency of the model to enable it to run on
mobile devices such as mobile phones.

## References

  * [1] James Atwood and Don Towsley.  Diffusion-convolutional neural networks.  In Advances in neural information processing systems, pages 1993–2001, 2016.
  * [2] Lei Bai, Lina Yao, Can Li, Xianzhi Wang, and Can Wang.  Adaptive graph convolutional recurrent network for traffic forecasting.  arXiv preprint arXiv:2007.02842, 2020.
  * [3] Donald J Berndt and James Clifford.  Using dynamic time warping to find patterns in time series.  In KDD workshop, volume 10, pages 359–370. Seattle, WA, USA:, 1994\.
  * [4] Joan Bruna, Wojciech Zaremba, Arthur Szlam, and Yann LeCun.  Spectral networks and locally connected networks on graphs.  arXiv preprint arXiv:1312.6203.
  * [5] Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud.  Neural ordinary differential equations.  arXiv preprint arXiv:1806.07366, 2018.
  * [6] Kyunghyun Cho, Bart Van Merriënboer, Caglar Gulcehre, Dzmitry Bahdanau, Fethi Bougares, Holger Schwenk, and Yoshua Bengio.  Learning phrase representations using rnn encoder-decoder for statistical machine translation.  arXiv preprint arXiv:1406.1078, 2014.
  * [7] Jared Quincy Davis, Krzysztof Choromanski, Jake Varley, Honglak Lee, Jean-Jacques Slotine, Valerii Likhosterov, Adrian Weller, Ameesh Makadia, and Vikas Sindhwani.  Time dependence in non-autonomous neural odes.  arXiv preprint arXiv:2005.01906, 2020.
  * [8] Michaël Defferrard, Xavier Bresson, and Pierre Vandergheynst.  Convolutional neural networks on graphs with fast localized spectral filtering.  arXiv preprint arXiv:1606.09375, 2016.
  * [9] Jianzhun Du, Joseph Futoma, and Finale Doshi-Velez.  Model-based reinforcement learning for semi-markov decision processes with neural odes.  arXiv preprint arXiv:2006.16210, 2020.
  * [10] Emilien Dupont, Arnaud Doucet, and Yee Whye Teh.  Augmented neural odes.  arXiv preprint arXiv:1904.01681, 2019.
  * [11] Xu Geng, Yaguang Li, Leye Wang, Lingyu Zhang, Qiang Yang, Jieping Ye, and Yan Liu.  Spatiotemporal multi-graph convolution network for ride-hailing demand forecasting.  In Proceedings of the AAAI conference on artificial intelligence, volume 33, pages 3656–3663, 2019.
  * [12] Arnab Ghosh, Harkirat Singh Behl, Emilien Dupont, Philip HS Torr, and Vinay Namboodiri.  Steer: Simple temporal regularization for neural odes.  arXiv preprint arXiv:2006.10711, 2020.
  * [13] Justin Gilmer, Samuel S Schoenholz, Patrick F Riley, Oriol Vinyals, and George E Dahl.  Neural message passing for quantum chemistry.  In International Conference on Machine Learning, pages 1263–1272. PMLR, 2017.
  * [14] William L Hamilton, Rex Ying, and Jure Leskovec.  Inductive representation learning on large graphs.  arXiv preprint arXiv:1706.02216, 2017.
  * [15] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770–778, 2016.
  * [16] Yujun Huang, Yunpeng Weng, Shuai Yu, and Xu Chen.  Diffusion convolutional recurrent neural network with rank influence learning for traffic forecasting.  In 2019 18th IEEE International Conference On Trust, Security And Privacy In Computing And Communications/13th IEEE International Conference On Big Data Science And Engineering (TrustCom/BigDataSE), pages 678–685. IEEE, 2019.
  * [17] David Kanaa, Vikram Voleti, Samira Kahou, Christopher Pal, and Mila CIFAR.  Simple video generation using neural odes.  In Annu. Conf. Neural Inf. Process. Syst, 2019.
  * [18] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [19] Thomas N Kipf and Max Welling.  Semi-supervised classification with graph convolutional networks.  arXiv preprint arXiv:1609.02907, 2016.
  * [20] Ruoyu Li, Sheng Wang, Feiyun Zhu, and Junzhou Huang.  Adaptive graph convolutional neural networks.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 32, 2018.
  * [21] Yaguang Li, Rose Yu, Cyrus Shahabi, and Yan Liu.  Diffusion convolutional recurrent neural network: Data-driven traffic forecasting.  arXiv preprint arXiv:1707.01926, 2017.
  * [22] Lingbo Liu, Jingwen Chen, Hefeng Wu, Jiajie Zhen, Guanbin Li, and Liang Lin.  Physical-virtual collaboration modeling for intra-and inter-station metro ridership prediction.  IEEE Transactions on Intelligent Transportation Systems, 2020.
  * [23] Lingbo Liu, Ruimao Zhang, Jiefeng Peng, Guanbin Li, Bowen Du, and Liang Lin.  Attentive crowd flow machines.  In Proceedings of the 26th ACM international conference on Multimedia, pages 1553–1561, 2018.
  * [24] Yang Liu, Zhiyuan Liu, and Ruo Jia.  Deeppf: A deep learning based architecture for metro passenger flow prediction.  Transportation Research Part C: Emerging Technologies, 101:18–34, 2019.
  * [25] Yiping Lu, Aoxiao Zhong, Quanzheng Li, and Bin Dong.  Beyond finite layer neural networks: Bridging deep architectures and numerical differential equations.  In International Conference on Machine Learning, pages 3276–3285. PMLR, 2018.
  * [26] Stefano Massaroli, Michael Poli, Michelangelo Bin, Jinkyoo Park, Atsushi Yamashita, and Hajime Asama.  Stable neural flows.  arXiv preprint arXiv:2003.08063, 2020.
  * [27] Miloš Milenković, Libor Švadlenka, Vlastimil Melichar, Nebojša Bojović, and Zoran Avramović.  Sarima modelling approach for railway passenger flow forecasting.  Transport, 33(5):1113–1120, 2018.
  * [28] Hans Pinckaers and Geert Litjens.  Neural ordinary differential equations for semantic segmentation of individual colon glands.  arXiv preprint arXiv:1910.10470, 2019.
  * [29] Yulia Rubanova, Ricky TQ Chen, and David Duvenaud.  Latent odes for irregularly-sampled time series.  arXiv preprint arXiv:1907.03907, 2019.
  * [30] Michael Schlichtkrull, Thomas N Kipf, Peter Bloem, Rianne Van Den Berg, Ivan Titov, and Max Welling.  Modeling relational data with graph convolutional networks.  In European semantic web conference, pages 593–607. Springer, 2018\.
  * [31] Man-Chun Tan, Sze Chun Wong, Jian-Min Xu, Zhan-Rong Guan, and Peng Zhang.  An aggregation approach to short-term traffic flow prediction.  IEEE Transactions on Intelligent Transportation Systems, 10(1):60–69, 2009.
  * [32] Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Lio, and Yoshua Bengio.  Graph attention networks.  arXiv preprint arXiv:1710.10903, 2017.
  * [33] E Weinan.  A proposal on machine learning via dynamical systems.  Communications in Mathematics and Statistics, 5(1):1–11, 2017.
  * [34] Billy M Williams and Lester A Hoel.  Modeling and forecasting vehicular traffic flow as a seasonal arima process: Theoretical basis and empirical results.  Journal of transportation engineering, 129(6):664–672, 2003.
  * [35] Zonghan Wu, Shirui Pan, Fengwen Chen, Guodong Long, Chengqi Zhang, and S Yu Philip.  A comprehensive survey on graph neural networks.  IEEE transactions on neural networks and learning systems, 2020\.
  * [36] Zonghan Wu, Shirui Pan, Guodong Long, Jing Jiang, and Chengqi Zhang.  Graph wavenet for deep spatial-temporal graph modeling.  arXiv preprint arXiv:1906.00121, 2019.
  * [37] Louis-Pascal Xhonneux, Meng Qu, and Jian Tang.  Continuous graph neural networks.  In International Conference on Machine Learning, pages 10432–10441. PMLR, 2020.
  * [38] Hanshu Yan, Jiawei Du, Vincent YF Tan, and Jiashi Feng.  On robustness of neural ordinary differential equations.  arXiv preprint arXiv:1910.05513, 2019.
  * [39] Huaxiu Yao, Xianfeng Tang, Hua Wei, Guanjie Zheng, and Zhenhui Li.  Revisiting spatial-temporal similarity: A deep learning framework for traffic prediction.  In Proceedings of the AAAI conference on artificial intelligence, volume 33, pages 5668–5675, 2019.
  * [40] Huaxiu Yao, Fei Wu, Jintao Ke, Xianfeng Tang, Yitian Jia, Siyu Lu, Pinghua Gong, Jieping Ye, and Zhenhui Li.  Deep multi-view spatial-temporal network for taxi demand prediction.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 32, 2018.
  * [41] Cagatay Yildiz, Markus Heinonen, and Harri Lähdesmäki.  Ode2vae: Deep generative second order odes with bayesian neural networks.  2019\.
  * [42] Rex Ying, Jiaxuan You, Christopher Morris, Xiang Ren, William L Hamilton, and Jure Leskovec.  Hierarchical graph representation learning with differentiable pooling.  arXiv preprint arXiv:1806.08804, 2018.
  * [43] Muhan Zhang and Yixin Chen.  Link prediction based on graph neural networks.  arXiv preprint arXiv:1802.09691, 2018.
  * [44] Ling Zhao, Yujiao Song, Chao Zhang, Yu Liu, Pu Wang, Tao Lin, Min Deng, and Haifeng Li.  T-gcn: A temporal graph convolutional network for traffic prediction.  IEEE Transactions on Intelligent Transportation Systems, 21(9):3848–3858, 2019.
  * [45] Fan Zhou, Liang Li, Kunpeng Zhang, and Goce Trajcevski.  Urban flow prediction with spatial–temporal neural odes.  Transportation Research Part C: Emerging Technologies, 124:102912, 2021.
  * [46] Fan Zhou, Liang Li, Ting Zhong, Goce Trajcevski, Kunpeng Zhang, and Jiahao Wang.  Enhancing urban flow maps via neural odes.