# Pedestrian Trajectory Prediction Using Dynamics-based Deep Learning

Honghui Wang1, Weiming Zhi2, Gustavo Batista3, Rohitash Chandra1
Correspondence to: H. Wang<EMAIL_ADDRESS>Transitional Artificial
Intelligence Research Group, School of Mathematics and Statistics, UNSW
Sydney, Australia2 Robotics Institute, Carnegie Mellon University, Pittsburgh,
USA3 School of Computer Science and Engineering, UNSW Sydney, Australia

###### Abstract

Pedestrian trajectory prediction plays an important role in autonomous driving
systems and robotics. Recent work utilising prominent deep learning models for
pedestrian motion prediction makes limited a priori assumptions about human
movements, resulting in a lack of explainability and explicit constraints
enforced on predicted trajectories. This paper presents a dynamics-based deep
learning framework where a novel asymptotically stable dynamical system is
integrated into a deep learning model. Our novel asymptotically stable
dynamical system is used to model human goal-targeted motion by enforcing the
human walking trajectory converges to a predicted goal position and provides a
deep learning model with prior knowledge and explainability. Our deep learning
model utilises recent innovations from transformer networks and is used to
learn some features of human motion, such as collision avoidance, for our
proposed dynamical system. The experimental results show that our framework
outperforms recent prominent models in pedestrian trajectory prediction on
five benchmark human motion datasets.

##  I INTRODUCTION

Human motion analysis is required for safe human-human and human-robot
interaction systems. Human trajectory prediction plays an important part in
human motion analysis [1] and therefore has been widely used in various fields
such as autonomous driving systems [2][3] and robot navigation [4].

The research on predicting human movements starts from physics-based methods,
such as the social force model [5] and the constant velocity/acceleration
model [6]. In recent years, deep learning methods have become appealing for
pedestrian (human) trajectory prediction with the public availability of
large-scale data [7]. In these approaches, pedestrian motion patterns
extracted from trajectories of previously observed humans are used to predict
human future movements. Deep learning methods applied to pedestrian trajectory
prediction have evolved from Recurrent Neural Networks (RNN) and their
variants [8] to Transformer models [9]. Recent efforts also have found notable
performance improvements by encoding goals (position or desired location) in
the deep neural network together with historically observed trajectories[10].

However, since coefficients cannot be physically understood and interpreted in
neural networks, the above-mentioned deep learning methods lack
explainability, making it unclear why the predicted trajectories have the
given shapes [11]. In addition, there is a lack of constraints which are
explicitly enforced when predicting desired trajectories [12]. The lack of
explainability and the explicit enforcement of known rules of desired
predicted trajectories can be addressed by introducing prior knowledge into
the deep learning methods [12][13].

Therefore, in this paper, we propose a dynamics-based deep learning method by
introducing a novel asymptotically stable dynamical system into a learning-
based transformer model. Our proposed asymptotically stable dynamical system
is used to model human motion and can enforce the human trajectory to converge
to the equilibrium point [14]; i.e. the destination of the pedestrian in the
domain of pedestrian trajectory prediction. This property is not only in line
with the goal-targeted feature of human movements, as depicted in Fig. 1, but
also can provide our used transformer network with prior knowledge and
explainability. In addition, our proposed method learns regularity properties
of human motion such as continuity, smoothness and boundness [15], and the
temporal and spatial reaction of human movements via the transformer-based
model. Overall, our method has the potential to forecast pedestrian
trajectories more precisely than conventional methods by accurately capturing
the properties of human motion.

Figure 1: Pedestrian trajectories are generally _goal-driven_ , that is
pedestrians seek to arrive around a determined goal.

To summarize, the main contributions of this paper are as follows:

(1) We propose to integrate asymptotically stable dynamical systems within a
neural network to model the human goal-directed motion by, and thereby
enforcing the human walking trajectory converges to a predicted goal position,
and providing a neural network with prior knowledge and explainability.

(2) We provide extensive evaluation on our proposed method against established
motion prediction benchmark datasets.

##  II Related Work

Human motion prediction is a key part of autonomous entities, such as robots
and driverless vehicles. The development of deep learning results in learning-
based methods where the human observed walking trajectory, seen as a temporal
sequence, is used to learn motion patterns for predicting future positions
[12]. In particular, long short-term memory (LSTM) [16] model, which is a type
of recurrent neural network (RNN), has been widely used to represent the
temporal information and spatial information of human movements in different
scenes accordingly [17] [18] [19] [20]. Transformer networks, as predominant
methods in Natural Language Processing domains (NLP) [21][22], also have been
applied to pedestrian trajectory prediction. The transformer-based methods
show better prediction accuracy of human future movements than LSTM-based
approaches via the self-attention mechanism for capturing temporal dependency
directly and spatial information of human movements indirectly [23] [24][25].

Encoding goal information in neural network models helps predict human motion
more precisely by reducing the uncertainty in human motion since humans are
always walking toward their destination. There are two predominant ways to use
goal information in deep learning methods. The first approach is to model the
relationship between semantic goals via computer vision methods, such as
right-turn, to achieve trajectory prediction [26] [27]. The other research
branch is to use the goal position derived by the expert system or sampling in
neural network models [28] [29] [10].

Deep learning approaches have shown state-of-the-art performance in
forecasting human trajectories[30]. However, limited a priori assumptions
about human movements in these methods make them lack explicit constraints on
predicted trajectories and be hard to explain. Explainability and explicit
constraints are essential for autonomous entities since these allow autonomous
entities to analyze their behaviour, attribute responsibility, and engender
trust [31]. We introduce our proposed asymptotically stable dynamical system
into a transformer-based method – Spatio-Temporal grAph tRansformer (STAR)[23]
for imposing explicit constraints and providing prior knowledge and
explainability. In our proposed algorithm, the goal estimator from [10] is
used to estimate the endpoint in advance for our proposed dynamical system,
and the STAR method is directly used to learn an asymptotically stable
dynamical system that is behaviourally suitable to human motion in the real
world so as to predict the future human walking trajectory more precisely.

##  III Dynamics-based deep learning

###  III-A Problem formulation

Given the previously observed trajectories of $M$ pedestrians in 2D
environments during time steps $1$ to $t_{obs}$, the task is to forecast the
correct future human motion trajectories during time steps $t_{obs+1}$ to
$t_{end}$. In this paper, unless otherwise specified, we denote $i^{th}$
pedestrian trajectory $\xi_{i}(t^{\prime\prime},p_{i}(t^{\prime}))$ starting
from $t^{\prime}$ to $t^{\prime\prime}$ by a set of
$\\{p_{i}(t)\\}|_{t=t^{\prime}}^{t=t^{\prime\prime}}$, where
$p_{i}(t)=(x_{i}(t),y_{i}(t))\in\mathbb{R}^{2}$ is the position of $i^{th}$
pedestrian at time $t\in\\{t^{\prime},t^{\prime}+1,..,t^{\prime\prime}\\}$,
and $i\in\\{1,2....,M\\}$.

###  III-B Framework

Figure 2: Our framework uses the goal-shift algorithm to replace the original
input to the transformer network. The figure shows the one-step prediction
only.

Our framework consists of three parts as shown in Fig. 2: First, we use a
separate goal estimator to infer the end position for our proposed algorithm.
Second, we propose a new asymptotically stable dynamical system where
endpoints are estimated by the goal estimator to construct human motion
representation. Finally, we use a deep learning model and train our proposed
dynamical system to obtain a well-developed representation for human walking
prediction.

###  III-C Human motion representation

The $i^{th}$ human motion, represented by a mapping between the position
$p_{i}(t)\in\mathbb{R}^{2}$ and its velocity $v_{i}(t)\in\mathbb{R}^{2}$, can
be modelled as a time-invariant first-order dynamical system [32]:

$\displaystyle v_{i}(t)=f(p_{i}(t))$ (1)

where $f:\mathbb{R}^{2}\rightarrow\mathbb{R}^{2}$ is a non-linear mapping,
representing a dynamical system, and $p_{i}(1)$ is the initial position. The
trajectory of human motion can be represented as:

$\displaystyle\xi_{i}(t,p_{i}(1))=p_{i}(1)+\int_{1}^{t}v_{i}(s)ds$ (2)

where $\xi_{i}(t,p_{i}(1))$ represents the trajectory of $i^{th}$ human motion
starting from the initial time.

A system (human motion) $f$ is locally asymptotically stable in a region
$S\subset R^{2}$ if trajectories starting in $S$ converge to equilibrium
points $p_{i}^{\ast}=p_{i}(t_{end})\in R^{2}$, a point where the velocity is
zero ($f(p_{i}^{\ast})=0$), that is

$\displaystyle\lim_{t\to\infty}\xi_{i}(t,p_{i}(1))=p_{i}^{\ast},\forall
p_{i}(1)\in S$ (3)

###  III-D Goal estimator

As mentioned earlier, we use the goal estimator from [10] to predict endpoints
for pedestrians in advance for determining equilibrium points. The process is
composed of three parts: $\gamma$-Soft-Dynamic Time Warping (DTW) [33] for
matching similar trajectories, construction of the expert repository by
selecting the $N$ most similar trajectories, and prediction of endpoints by
K-means clustering on the expert repository.

####  III-D1 $\gamma$-soft-DTW algorithm

Unlike the DTW algorithm [34], the $\gamma$-soft-DTW algorithm is
differentiable and, therefore, widely used in many applications related to
time series. The formula of $\gamma$-soft-DTW is shown as follows:

$\begin{split}D(X,Y)&=DTW_{\gamma}(X,Y)\\\
&=min_{\gamma}\\{<A,\Delta(X,Y)>\\}\end{split}$ (4)

where $\Delta(\cdot)$ is the cost function to measure the Euclidean distance
between the temporal sequence $X\in R^{p\times n}$ and the temporal sequence
$Y\in R^{p\times m}$. $A\in R^{n\times m}$ is the alignment matrix.
$<A,\Delta(X,Y)>$ is the inner product of the cost matrix $\Delta(X,Y)$ with
the alignment matrix, reflecting the similarity between two temporal
sequences. $min_{\gamma}\\{\cdot\\}$ is $\gamma$-soft-DTW operator, defined as
follows:

$\displaystyle
min_{\gamma}(a_{1},a_{2},...a_{n})=\left\\{\begin{array}[]{rcl}min_{i<=n}a_{i}&&{\gamma=0}\\\
-\gamma log\sum_{i=1}^{n}e^{-a_{i}/\gamma}&&{\gamma>0}\end{array}\right.$ (7)

In addition to the relative positions of the neighbours, humans naturally tend
to focus on the neighbours’ velocities when navigating in crowded environments
[35]. Compared with absolute coordinates, the velocity is more informative.
Therefore, unlike the previous work [10], in our proposed algorithm, we use
pedestrian velocity as the input of $\gamma$-soft-DTW algorithm for a similar
search between two temporal sequences, and $\gamma=1$.

####  III-D2 Construction of expert repository

According to the result calculated by the $\gamma$-soft-DTW algorithm, the $N$
most similar trajectories are used to construct the expert repository. In our
proposed algorithm, we normalize all trajectories in the expert repository for
more accurate prediction, shown as follows:

$\begin{split}\tilde{\xi}_{i}(t,\tilde{p}_{i}(1))=\xi_{i}(t,p_{i}(1))-p_{i}(1)\end{split}$
(8)

where $\tilde{p}_{i}(t)$ represents the normalized position of the $i^{th}$
pedestrian from the expert expository at time $t\in\\{1,2,.....t_{end}\\}$. In
our proposed algorithm, $N=100$.

####  III-D3 Goal estimation by k-means clustering

K-means clustering[36] is used to cluster all ground truth endpoints of the
normalized trajectories from the expert repository, and $K$ cluster centres
are regarded as all candidate predicted endpoints. Finally, the cluster canter
with a minimum distance from the ground truth endpoint on the tested
trajectory is selected as the estimated endpoint.

###  III-E A novel asymptotically stable dynamical system

Human are always moving towards their destinations, and human motion
trajectories are generally irregular. In order to represent human movements
more precisely, a new dynamical system with the guarantee of stability
property is proposed in this paper. Our proposed dynamical system is based on
a stable attractor and natural gradient descent [37].

####  III-E1 Stable attractor

The stable attractor in Euclidean space is used in the gradient descent
dynamical system for the guarantee of stability property. The stable attractor
is shown as follows:

$\displaystyle\Phi(p_{i}(t))=||p_{i}(t)-p_{i}^{\star}||_{2}$ (9)

where $p_{i}^{\star}$ is the equilibrium point, indicating the destination of
the $i^{th}$ pedestrian.

The gradient descent dynamical system that uses the stable attractor is shown
as follows:

$\displaystyle v(p_{i}(t)):=-\nabla_{p_{i}(t)}\Phi(p_{i}(t))$ (10)

This dynamical system can be guaranteed to be stable since the stable
attractor is a valid Lyapunov function [38] for this system. Geometrically,
the human motion trajectory generated by this stable dynamical system is
straight-line with the unit-velocity length and converges to the equilibrium
point after a certain time.

####  III-E2 Natural gradient descent

Generally, human movement trajectories are generally curve rather than
straight lines since they avoid collision with obstacles or pedestrians by
altering speed and direction at any time when humans are walking toward their
destinations. It means that the stable gradient descent dynamical system is
hard to represent human complex motion. In order to imitate various types of
pedestrian movements as much as possible, we propose a new asymptotically
stable dynamical system based on the natural gradient descent dynamics [37],
shown as follows:

$\displaystyle v(p_{i}(t)):=-P(p_{i}(t))\nabla_{p_{i}(t)}\Phi(p_{i}(t))$ (11)

Our proposed dynamical system can generate the curved trajectory because of
natural gradient descent dynamics that are derived from the stable attractor
$\Phi(p_{i}(t))$. In addition, our proposed dynamical system is able to
maintain stability when $P(p_{i}(t))$ is a positive-definite matrix.

In our proposed stable dynamical system, the walking velocity of a person,
which dominates human walking motion behaviours, can be determined by the
positive-definite matrix $P(p_{i}(t))$. In order to imitate various behaviours
of human movement more realistically, the positive-definite matrix
$P(p_{i}(t))$ needs to be endowed with some properties of human motion such as
collision avoidance.

###  III-F Dynamics-based Trajectory Predictor

The irregularity of human movements makes it challenging to determine
positive-definite matrix $P(p_{i}(t))$. Recently, data-driven methods have
been widely used to represent complex behaviours of human motion. In this
paper, a transformer [23] is used to endow the positive-definite matrix
$P(p_{i}(t))$ with the feature of human motion in a data-driven way.
Specifically, the goal-shift algorithm is firstly conducted on the pedestrian
trajectory, and the processed trajectory is input to the transformer network.
Secondly, our asymptotically stable dynamical system is trained by training
positive-definite matrix $P(p_{i}(t))$ via the transformer-based method with
random Gaussian noise injection for generating the velocity that can predict
future positions of pedestrians. Finally, the loss function is determined to
imitate the human motion trajectory in the real world as accurately as
possible.

####  III-F1 Goal shift

Sightly different from [10], in our proposed algorithm, the whole goal-shifted
human motion trajectory is used as an input. The goal-shift algorithm on the
human walking trajectory is shown as follows:

$\begin{split}\bar{{\xi}}_{i}(t,\bar{p}_{i}(1))=\xi_{i}(t,p_{i}(1))-p_{i}(t_{end})\end{split}$
(12)

where $\bar{p}_{i}(t)$ represents the goal-shifted position of the $i^{th}$
pedestrian at time $t\in\\{1,2....,t_{end}\\}$.

During the training phase, $p_{i}(t_{end})=(x_{i}(t_{end}),y_{i}(t_{end}))$ is
the ground truth endpoint, while it is predicted by goal estimator we use
during the test phase.

####  III-F2 Training the novel asymptotically stable dynamical system

In our work, the positive-definite matrix $P(p_{i}(t))$ is decomposed based on
lower–upper (LU) decomposition for using data-driven methods to learn
behaviours of human movements for our proposed dynamical system while keeping
the positive-definite property, shown as follows:

$\displaystyle P(p_{i}(t))=L(p_{i}(t))L^{T}(p_{i}(t))+\sigma I$ (13)

where $L(p_{i}(t))$ is a lower triangular matrix. $I$ is the identity matrix.
We set $\sigma$ to be $10^{-8}$.

Human walking motion trajectories can be represented as two-dimension temporal
sequences. It is clear that for any
$a(p_{i}(t)),b(p_{i}(t)),c(p_{i}(t))\in\mathbb{R}$ in the lower triangular
matrix $L(p_{i}(t))\in\mathbb{R}^{2\times 2}$,
$P(p_{i}(t))\in\mathbb{R}^{2\times 2}$ can keep positive-definite property.

Therefore, learning the positive-definite matrix $P(p_{i}(t))$ can be
transformed into learning three elements in the lower triangular matrix
decomposed from the above formula. In our proposed algorithm, the transformer
model predicts the positive-definite matrix that shapes the motion, by
learning three elements in the lower triangular matrix. Specifically,

$f_{T}(p_{i}(t))=L(p_{i}(t))=\begin{bmatrix}a(p_{i}(t))&0\\\
b(p_{i}(t))&c(p_{i}(t))\end{bmatrix}$ (14)

where $f_{T}$ is the transformer-based method.

In our proposed algorithm, the goal-shifted position is inputted to the
transformer-based method. Pedestrian future positions can be predicted by
integrating over the velocity derived by our well-trained asymptotically
stable dynamical system. The generated pedestrian trajectories can match the
real human motion closely for more precise prediction.

####  III-F3 Loss function

A mean squared error (MSE) loss between generated trajectories by our proposed
algorithm and ground truth trajectories is optimised during training. In
particular, we have,

$\displaystyle\mathcal{L}_{\text{MSE}}=\frac{1}{M*(t_{end}-1)}\sum_{i=1}^{M}\sum_{t=2}^{t_{end}}||p_{i}(t)-\hat{p}_{i}(t)||_{2}^{2}$
(15)

where $\hat{p}_{i}(t)$ represents the predicted position of the $i^{th}$
pedestrian at time $t$.

##  IV Empirical Evaluation

###  IV-A Setup

The ETH/UCY dataset is a collection of benchmark datasets commonly used for
evaluating trajectory prediction algorithms. The dataset is divided into two
parts: the ETH dataset [39] and the UCY dataset [40]. The ETH dataset includes
three video sequences, each with a different scene and a different number of
pedestrians, ranging from 28 to 60 individuals. The UCY dataset includes five
video sequences, each with a different scene and a different number of
pedestrians, ranging from 30 to 94 individuals.

Two well-known metrics, Average Displacement Error (ADE) and Final
Displacement Error (FDE), are commonly used for assessing prediction accuracy.
ADE and FDE are given as:

$ADE=\frac{\sum_{i=1}^{M}\sum_{t=t_{obs}+1}^{t_{end}}\|p_{i}(t)-\hat{p}_{i}(t)\|_{2}}{M\times
T}$ (16)

$FDE=\frac{\sum_{i=1}^{M}\|p_{i}(t_{end})-\hat{p}_{i}(t_{end})\|_{2}}{M}$ (17)

where $T$ is the number of predicted timesteps, and $\hat{p}_{i}(t)$
represents the predicted position of the $i^{th}$ pedestrian at time $t$.

In the ETH/UCY dataset, we use a leave-one-out approach similar to previous
works [41], where the model trains on four datasets and tests on the remaining
set. 8 frames (3.2s) as an observed sequence and 12 frames (4.8s) as a
predicted sequence are used for evaluation.

###  IV-B Implementation Details

We use the same data processing method as [23] for training and testing with
the same hyperparameters, while the training epoch is set to 200.

TABLE I: Quantitative results of trajectory prediction for comparison on ETH/UCY. The lower error indicates better performance. The best performance is marked in bold. Method | minimum ADE/FDE of 20 predicted samples
---|---
Dataset | ETH | HOTEL | UNIV | ZARA1 | ZARA2 | AVERAGE
Social-STGCNN [42] | 0.64/1.11 | 0.49/0.85 | 0.44/0.79 | 0.34/0.53 | 0.30/0.48 | 0.44/0.75
PECNet [29] | 0.54/0.87 | 0.18/0.24 | 0.35/0.60 | 0.22/0.39 | 0.17/0.30 | 0.29/0.48
Constant Velocity with noise [43] | 0.43/0.80 | 0.19/0.35 | 0.34/0.71 | 0.24/0.48 | 0.21/0.45 | 0.28/0.56
STAR [23] | 0.36/0.65 | 0.17/0.36 | 0.31/0.62 | 0.26/0.55 | 0.22/0.46 | 0.26/0.53
Goal-GAN [28] | 0.59/1.18 | 0.19/0.35 | 0.60/1.19 | 0.43/0.87 | 0.32/0.65 | 0.43/0.85
AVGCN [44] | 0.62/1.06 | 0.31/0.58 | 0.55/1.20 | 0.33/0.70 | 0.27/0.58 | 0.42/0.82
Social-Implicit [45] | 0.66/1.44 | 0.20/0.36 | 0.31/0.60 | 0.25/0.50 | 0.22/0.43 | 0.33/0.67
GA-STT [25] | 0.51/0.89 | 0.22/0.46 |  0.29/0.63 | 0.25/0.55 | 0.20/0.44 | 0.29/0.59
DDL (Ours) | 0.26/0.50 |  0.15/0.35 |  0.29/0.58 | 0.16/0.29 | 0.13/0.22 | 0.20/0.39

###  IV-C Quantitative Results and Analysis

We compare our proposed method, i.e. dynamics-based deep learning (DDL), with
a wide range of methods, including Social Spatio-Temporal Graph Convolutional
Neural Network (Social-STGCNN) [42], Predicted Endpoint Conditioned Network
(PECNet) [29], Constant velocity with noise [43], Spatio-Temporal grAph
tRansformer (STAR) [23], Goal-Generative Adversarial Networks (Goal-GAN) [28],
Attention Visual Field Constraints Graph Convolutional Networks (AVGCN) [44],
Social-Implicit [45], and Group-Aware Spatial-Temporal Transformer (GA-STT)
[25]. We use the $Best-of-N$ method, that is the minimum ADE and FDE from $N$
randomly sampled trajectories. The respective methods predict multiple future
trajectories by sampling 20 times. We select the best prediction from 20
randomly-sampled trajectories and present the results on ETH/UCY datasets in
Table I.

We can see from Table I that DDL performs better than all compared methods on
most datasets. Compared with STAR [23] method, DDL achieves better results for
all datasets and has an average increase of 23.1% in ADE and 26.4% in FDE. The
remarkable improvement suggests that the goal-shift method and our proposed
asymptotically stable dynamical system, as a provider of prior knowledge and
explainability to STAR, contribute to forecasting human movements more
accurately. Even though our proposed approach reports the same result of ADE
on UNIV dataset with high crowd density as GA-STT which considers group
interactions in the spatial-temporal transformer network[25], our method
outperforms 31.0% and 33.9% in ADE and FDE, respectively, on average. PECNet
[29] sets the average endpoint between predicted endpoints sampled from latent
future endpoints and ground truth endpoints in the loss function to learn
parameters for estimating future endpoints. In the case of the Hotel dataset
where the percentage of pedestrians who are almost waiting or ambling is high,
better performance on FDE can be achieved by PECNet via learning scene-
matching pedestrian movement patterns. Nevertheless, DDL shows an average
improvement of 31.0% on ADE and 18.8% on FDE.

###  IV-D Qualitative Results and Analysis

In this section, we present selected trajectories generated by DDL to gain
more insights into the actual various behaviours of human walking, indicating
that our proposed algorithm is able to imitate a wide range of human movements
and therefore predict socially acceptable trajectories.

####  IV-D1 Convergence to destination

DDL can make trajectories converge to endpoints because of the stability
property that existed in our proposed dynamical system. We can see (Fig. 3a
and 3b) that the convergence property can be obviously reflected by DDL on
single and group human movement predictions, showing that our asymptotically
stable dynamical system provides the transformer with prior knowledge and
explainability and imposes explicit constraints on predicted trajectories by
integrating stability theory into the deep learning model.

(a) Single person walking

(b) Group people walking

Figure 3: Trajectory visualization of human goal-targeted walking. Dots
construct observed trajectories while discrete lines construct predicted
trajectories. Continuous lines represent ground truth trajectories. The round
dots represent ground truth endpoints, while stars represent predicted
endpoints.

(a) Alter velocity

(b) Keep social distance

(c) Turn around

Figure 4: Trajectory visualization in the case of collision avoidance. Dots construct observed trajectories while discrete lines construct predicted trajectories. Continuous lines represent ground truth trajectories. The round dots represent ground truth endpoints, while stars represent predicted endpoints. TABLE II: We present quantitative results of different combinations of components from DDL tested on ETH/UCY datasets. Component A indicates the goal-shift module and Component B indicates our proposed asymptotically stable dynamical system in DDL. Finally, Component C is the STAR model. Note that the lower error indicates better performance, and we mark the best performance in bold. Components | minimum ADE/FDE of 20 predicted samples
---|---
A | B | C | ETH | HOTEL | UNIV | ZARA1 | ZARA2 | AVERAGE
|  | ✓ | 0.36/0.65 | 0.17/0.36 | 0.31/0.62 | 0.26/0.55 | 0.22/0.46 | 0.26/0.53
| ✓ | ✓ | 0.27/0.52 |  0.15/0.35 | 0.33/0.59 | 0.17/0.31 | 0.14/0.25 | 0.21/0.40
✓ |  | ✓ | 0.46/0.89 | 0.17/0.35 | 0.31/0.64 | 0.17/0.40 | 0.15/0.28 | 0.25/0.51
✓ | ✓ | ✓ |  0.26/0.50 |  0.15/0.35 |  0.29/0.58 |  0.16/0.29 |  0.13/0.23 |  0.20/0.39

####  IV-D2 Collision avoidance

Our results show that DDL can predict human movements that have the behaviour
of collision avoidance (Fig. 4). As shown in Fig. 4a, one person slows down to
avoid a collision when he/she sees another person who is walking forward in
front of him/her. Therefore, our proposed algorithm is able to imitate the one
behaviour of collision avoidance by changing the speed.

In Fig. 4b, when a group of people and another group of people walk in the
opposite direction, they keep a certain social distance from each other to
avoid collision, indicating that our proposed algorithm is able to generate
the human motion trajectory that can avoid collision by maintaining the social
distance.

We can see from Fig. 4c that when a group of people see another group of
people who are standing and a person who is walking in the opposite direction,
they alter their respective walking directions to avoid collision with this
group of people. It suggests that DDL is capable of capturing one property of
changing directions for avoiding collision.

###  IV-E Ablation Study

We conduct an ablation study to verify if our proposed asymptotically stable
dynamical system is the main contribution to the improvement of the
performance of pedestrian trajectory prediction in DDL. In our ablation study
experiments, we separately replace goal-shifted input with the original input
of [23] in DDL and use the goal-shifted input in [23]. The experiment result
conducted on ETH/UCY datasets is shown in Table II. We can see that compared
with [23], our proposed asymptotically stable dynamical system is the main
contribution (19.2% and 24.5% increases in ADE and FDE, respectively, on
average). In contrast, the goal-shift processing on the original trajectories
only shows a slight improvement in performance (both ADE and FDE increase by
3.8% on average).

##  V Conclusion and future research

We introduced a dynamics-based deep learning framework to generate
trajectories that meet the real human movements, for pedestrian trajectory
prediction. We proposed a novel asymptotically stable dynamical system to
model human goal-targeted movements with the help of a goal estimator, which
provides a transformer network with prior knowledge and explainability and
enforces explicit constraints on the prediction of trajectories, by
integrating stability theory. We used a transformer-based method to endow this
dynamical system with features of human real movements such as collision
avoidance in a learning-based way. Our method shows better performance than
related methods in human walking motion prediction on five benchmark
pedestrian walking datasets; i.e. by at least 23.1% and 18.8% in ADE and FDE,
respectively, on average.

Our method does have the ability to predict trajectories with limited motion
angle changes only due to the inherent properties of positive definite
matrices geometrically. In future research, diffeomorphism that can curve the
trajectory at a larger angle can be incorporated to overcome this limitation.

## Code and Data

We have released our code and dataset on GitHub. 111https://github.com/sydney-
machine-learning/pedestrianpathprediction

## References

  * [1] A. Rudenko, L. Palmieri, M. Herman, K. M. Kitani, D. M. Gavrila, and K. O. Arras, “Human motion trajectory prediction: A survey,” _The International Journal of Robotics Research_ , vol. 39, no. 8, pp. 895–935, 2020\.
  * [2] F. Leon and M. Gavrilescu, “A review of tracking and trajectory prediction methods for autonomous driving,” _Mathematics_ , vol. 9, no. 6, p. 660, 2021\.
  * [3] Y. Huang, J. Du, Z. Yang, Z. Zhou, L. Zhang, and H. Chen, “A survey on trajectory-prediction methods for autonomous driving,” _IEEE Transactions on Intelligent Vehicles_ , vol. 7, no. 3, pp. 652–674, 2022.
  * [4] W. Zhi, T. Lai, L. Ott, and F. Ramos, “Anticipatory navigation in crowds by probabilistic prediction of pedestrian future movements,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2021, pp. 8459–8464.
  * [5] D. Helbing and P. Molnar, “Social force model for pedestrian dynamics,” _Physical review E_ , vol. 51, no. 5, p. 4282, 1995.
  * [6] R. Schubert, E. Richter, and G. Wanielik, “Comparison and evaluation of advanced motion models for vehicle tracking,” in _2008 11th international conference on information fusion_. IEEE, 2008, pp. 1–6.
  * [7] B. I. Sighencea, R. I. Stanciu, and C. D. Căleanu, “A review of deep learning-based methods for pedestrian trajectory prediction,” _Sensors_ , vol. 21, no. 22, p. 7543, 2021.
  * [8] J. Schmidhuber, “Deep learning in neural networks: An overview,” _Neural networks_ , vol. 61, pp. 85–117, 2015.
  * [9] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, and I. Polosukhin, “Attention is all you need,” _Advances in neural information processing systems_ , vol. 30, 2017.
  * [10] H. Zhao and R. P. Wildes, “Where are you heading? dynamic trajectory prediction with expert goal examples,” in _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , 2021, pp. 7629–7638.
  * [11] R. Korbmacher and A. Tordeux, “Review of pedestrian trajectory prediction methods: Comparing deep learning and knowledge-based approaches,” _IEEE Transactions on Intelligent Transportation Systems_ , 2022.
  * [12] W. Zhi, L. Ott, and F. Ramos, “Probabilistic trajectory prediction with structural constraints,” in _2021 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2021, pp. 9849–9856.
  * [13] K. Beckh, S. Müller, M. Jakobs, V. Toborek, H. Tan, R. Fischer, P. Welke, S. Houben, and L. von Rueden, “Harnessing prior knowledge for explainable machine learning: An overview,” in _2023 IEEE Conference on Secure and Trustworthy Machine Learning (SaTML)_. IEEE, 2023, pp. 450–463.
  * [14] A. M. Lyapunov, “The general problem of the stability of motion,” _International journal of control_ , vol. 55, no. 3, pp. 531–534, 1992.
  * [15] T. Flash and N. Hogan, “The coordination of arm movements: an experimentally confirmed mathematical model,” _Journal of neuroscience_ , vol. 5, no. 7, pp. 1688–1703, 1985.
  * [16] S. Hochreiter and J. Schmidhuber, “Long short-term memory,” _Neural computation_ , vol. 9, no. 8, pp. 1735–1780, 1997.
  * [17] A. Alahi, K. Goel, V. Ramanathan, A. Robicquet, L. Fei-Fei, and S. Savarese, “Social lstm: Human trajectory prediction in crowded spaces,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2016, pp. 961–971.
  * [18] H. Xue, D. Q. Huynh, and M. Reynolds, “Ss-lstm: A hierarchical lstm model for pedestrian trajectory prediction,” in _2018 IEEE Winter Conference on Applications of Computer Vision (WACV)_. IEEE, 2018, pp. 1186–1194.
  * [19] A. Vemula, K. Muelling, and J. Oh, “Social attention: Modeling attention in human crowds,” in _2018 IEEE international Conference on Robotics and Automation (ICRA)_. IEEE, 2018, pp. 4601–4607.
  * [20] H. Cheng, W. Liao, M. Y. Yang, B. Rosenhahn, and M. Sester, “Amenet: Attentive maps encoder network for trajectory prediction,” _ISPRS Journal of Photogrammetry and Remote Sensing_ , vol. 172, pp. 253–266, 2021.
  * [21] Z. Lan, M. Chen, S. Goodman, K. Gimpel, P. Sharma, and R. Soricut, “Albert: A lite bert for self-supervised learning of language representations,” _arXiv preprint arXiv:1909.11942_ , 2019.
  * [22] Z. Yang, Z. Dai, Y. Yang, J. Carbonell, R. R. Salakhutdinov, and Q. V. Le, “Xlnet: Generalized autoregressive pretraining for language understanding,” _Advances in neural information processing systems_ , vol. 32, 2019.
  * [23] C. Yu, X. Ma, J. Ren, H. Zhao, and S. Yi, “Spatio-temporal graph transformer networks for pedestrian trajectory prediction,” in _Computer Vision–ECCV 2020: 16th European Conference, Glasgow, UK, August 23–28, 2020, Proceedings, Part XII 16_. Springer, 2020, pp. 507–523.
  * [24] F. Giuliari, I. Hasan, M. Cristani, and F. Galasso, “Transformer networks for trajectory forecasting,” in _2020 25th international conference on pattern recognition (ICPR)_. IEEE, 2021, pp. 10 335–10 342.
  * [25] L. Zhou, D. Yang, X. Zhai, S. Wu, Z. Hu, and J. Liu, “Ga-stt: Human trajectory prediction with group aware spatial-temporal transformer,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 3, pp. 7660–7667, 2022.
  * [26] A. Rasouli, I. Kotseruba, T. Kunic, and J. K. Tsotsos, “Pie: A large-scale dataset and models for pedestrian intention estimation and trajectory prediction,” in _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , 2019, pp. 6262–6271.
  * [27] N. Rhinehart, R. McAllister, K. Kitani, and S. Levine, “Precog: Prediction conditioned on goals in visual multi-agent settings,” in _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , 2019, pp. 2821–2830.
  * [28] P. Dendorfer, A. Osep, and L. Leal-Taixé, “Goal-gan: Multimodal trajectory prediction based on goal position estimation,” in _Proceedings of the Asian Conference on Computer Vision_ , 2020.
  * [29] K. Mangalam, H. Girase, S. Agarwal, K.-H. Lee, E. Adeli, J. Malik, and A. Gaidon, “It is not the journey but the destination: Endpoint conditioned trajectory prediction,” in _Computer Vision–ECCV 2020: 16th European Conference, Glasgow, UK, August 23–28, 2020, Proceedings, Part II 16_. Springer, 2020, pp. 759–776.
  * [30] X. Tan, R. Liu, S. Zhang, J. Li, and P. Ma, “Review of pedestrian trajectory prediction methods,” _Frontiers in Computing and Intelligent Systems_ , vol. 1, no. 3, pp. 68–77, 2022.
  * [31] F. Pasquale, “Toward a fourth law of robotics: Preserving attribution, responsibility, and explainability in an algorithmic society,” _Ohio St. LJ_ , vol. 78, p. 1243, 2017.
  * [32] W. Zhi, T. Lai, L. Ott, and F. Ramos, “Diffeomorphic transforms for generalised imitation learning.” in _L4DC_ , 2022, pp. 508–519.
  * [33] M. Cuturi and M. Blondel, “Soft-dtw: a differentiable loss function for time-series,” in _International conference on machine learning_. PMLR, 2017, pp. 894–903.
  * [34] H. Sakoe and S. Chiba, “Dynamic programming algorithm optimization for spoken word recognition,” _IEEE transactions on acoustics, speech, and signal processing_ , vol. 26, no. 1, pp. 43–49, 1978.
  * [35] P. Kothari, S. Kreiss, and A. Alahi, “Human trajectory forecasting in crowds: A deep learning perspective,” _IEEE Transactions on Intelligent Transportation Systems_ , vol. 23, no. 7, pp. 7386–7400, 2021.
  * [36] J. A. Hartigan and M. A. Wong, “Algorithm as 136: A k-means clustering algorithm,” _Journal of the royal statistical society. series c (applied statistics)_ , vol. 28, no. 1, pp. 100–108, 1979.
  * [37] M. A. Rana, A. Li, D. Fox, B. Boots, F. Ramos, and N. Ratliff, “Euclideanizing flows: Diffeomorphic reduction for learning stable dynamical systems,” in _Learning for Dynamics and Control_. PMLR, 2020, pp. 630–639.
  * [38] M. Vidyasagar, _Nonlinear systems analysis_. SIAM, 2002.
  * [39] S. Pellegrini, A. Ess, K. Schindler, and L. Van Gool, “You’ll never walk alone: Modeling social behavior for multi-target tracking,” in _2009 IEEE 12th international conference on computer vision_. IEEE, 2009, pp. 261–268.
  * [40] A. Lerner, Y. Chrysanthou, and D. Lischinski, “Crowds by example,” in _Computer graphics forum_ , vol. 26, no. 3. Wiley Online Library, 2007, pp. 655–664.
  * [41] A. Gupta, J. Johnson, L. Fei-Fei, S. Savarese, and A. Alahi, “Social gan: Socially acceptable trajectories with generative adversarial networks,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2018, pp. 2255–2264.
  * [42] A. Mohamed, K. Qian, M. Elhoseiny, and C. Claudel, “Social-stgcnn: A social spatio-temporal graph convolutional neural network for human trajectory prediction,” in _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , 2020, pp. 14 424–14 432.
  * [43] C. Schöller, V. Aravantinos, F. Lay, and A. Knoll, “What the constant velocity model can teach us about pedestrian motion prediction,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 2, pp. 1696–1703, 2020.
  * [44] C. Liu, Y. Chen, M. Liu, and B. E. Shi, “Avgcn: Trajectory prediction using graph convolutional networks guided by human attention,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2021, pp. 14 234–14 240.
  * [45] A. Mohamed, D. Zhu, W. Vu, M. Elhoseiny, and C. Claudel, “Social-implicit: Rethinking trajectory prediction evaluation and the effectiveness of implicit maximum likelihood estimation,” in _Computer Vision–ECCV 2022: 17th European Conference, Tel Aviv, Israel, October 23–27, 2022, Proceedings, Part XXII_. Springer, 2022, pp. 463–479.


# Surface Similarity Parameter: A New Machine Learning Loss Metric for
Oscillatory Spatio-Temporal Data

Mathies Wedler<EMAIL_ADDRESS>[ Merten Stender  Marco Klein  Svenja
Ehlers  Norbert Hoffmann  Hamburg University of Technology, Dynamics Group,
Schlossmühlendamm 30, 21073 Hamburg, Germany  Imperial College London,
Department of Mechanical Engineering, 58 Prince’s Gate, South Kensington,
London SW7 1AY, United Kingdom

###### Abstract

Supervised machine learning approaches require the formulation of a loss
functional to be minimized in the training phase. Sequential data are
ubiquitous across many fields of research, and are often treated with
Euclidean distance-based loss functions that were designed for tabular data.
For smooth oscillatory data, those conventional approaches lack the ability to
penalize amplitude, frequency and phase prediction errors at the same time,
and tend to be biased towards amplitude errors. We introduce the surface
similarity parameter (SSP) as a novel loss function that is especially useful
for training machine learning models on smooth oscillatory sequences. Our
extensive experiments on chaotic spatio-temporal dynamical systems indicate
that the SSP is beneficial for shaping gradients, thereby accelerating the
training process, reducing the final prediction error, and implementing a
stronger regularization effect compared to using classical loss functions. The
results indicate the potential of the novel loss metric particularly for
highly complex and chaotic data, such as data stemming from the nonlinear two-
dimensional Kuramoto–Sivashinsky equation and the linear propagation of
dispersive surface gravity waves in fluids.

###### keywords:

deep learning , loss function , error metric , similarity , nonlinear dynamics
, spatio-temporal dynamics

††journal: Neural Networks

url]www.tuhh.de/dyn

##  1 INTRODUCTION

Given the fact that almost any natural process is dynamic and non-stationary,
sequential data is ubiquitous across many scientific disciplines, such as
weather forecasting, finance, natural language processing, acoustics and many
more. Besides classical analytical and numerical techniques, machine learning
approaches have been adopted lately for various regression or classification
tasks on sequential data given a continuous quantity $y\left(t\right)$ sampled
at discrete instants $\mathbf{y}\left(t\right)=\left[y(t_{1}),\dots
y(t_{N})\right]$. Most data-driven approaches are rooted in an optimization
framework for model $f$ parameters $\mathbf{\theta}$, finally targeting at
minimizing an error functional
$J\left(\mathbf{y},\hat{\mathbf{y}}(\theta)\right)$ that describes the
discrepancy between the model prediction $\hat{\mathbf{y}}$ and the ground
truth $\mathbf{y}$. Some techniques utilize sequence featurization [Fulcher &
Jones, 2014], i.e. turning sequential data into vectors of descriptive scalar
values, while direct approaches perform the error computation on the raw
sequences themselves. There is a plethora of metrics for measuring the
similarity of two (multivariate) sequences, such as generic time domain
distance metrics, or customized metrics (e.g. error in frequency domain)
tailored towards specific aspects of the given data. The choice of a specific
error metric as loss function for the training of machine learning methods is
of crucial importance: the loss function shapes the gradient, and thus the
convergence behavior of the optimization.

Conventional Euclidean distance metrics such as the mean squared error and
mean absolute error do not take the sequential data characteristic into
account, nor smoothness or continuity of the predicted sequences. Instead,
considering sequences as values sampled from some distribution, the MSE aims
at fitting the mean, and the MAE at fitting the median of that distribution.
Often times, a spectral decomposition of complex (spatio-) temporal dynamics
is utilized to describe the data in a finite series of harmonics, each
described by amplitude, frequency and phase. From this viewpoint, Euclidean
distance metrics are biased towards amplitude errors, rendering them as sub-
optimal loss functions for some machine learning tasks as our work will
illustrate. We introduce the surface similarity parameter SSP [Perlin &
Bustamante, 2014] as a novel loss metric, that is capable of combining
amplitude, frequency and phase errors into a single scalar error metric.
Thereby, the SSP can provide several benefits for machine learning
applications on oscillatory sequential data.

Studying the highly chaotic dynamics of the a) one- and b) two-dimensional
Kuramoto–Sivashinsky equation (KS-eq.), and the c) propagation of dispersive
surface waves of a very broad spectrum, our work indicates that the SSP is a
particularly beneficial loss function the more complex the sequences are, that
is, the wider the underlying wave spectrum becomes (i.e. the larger the number
of relevant harmonics in the spectral decomposition). Our numerical
experiments indicate that the SSP enhances optimizer convergence through
effective gradient shaping, increases robustness against broad-banded spectra,
and introduces strong regularization effects. Further, the benefits of the SSP
loss function are shown to be rooted in an error saturation effect, that is
effectively bounding the degree to which two sequences can be dissimilar,
thereby forming sharper global minima in a flattened error functional in the
parameter optimization space.

This work is structured as follows: Section 2 introduces the surface
similarity parameter and re-visits major considerations on machine learning
loss functions. Section 3 illustrates several in-depth studies comparing
different loss metrics for learning a convolutional encoding-decoding flow map
for the KS-eq. and a broad-banded dispersive wave spectrum.

##  2 Methods

We will introduce the SSP first, then compare against conventional Euclidean
distance metrics, and finally highlight the saturation property of the
proposed loss function.

###  2.1 Surface similarity parameter

Perlin & Bustamante [2014] introduce the surface similarity parameter (SSP) as
a quantitative comparison metric for two surfaces
$\mathbf{y}_{1},\mathbf{y}_{2}$ with spatial and/or temporal extent in one or
more dimensions. In this work, we only consider surfaces with spatial extent,
such that the SSP formulation can be written as

$J_{\mathrm{SSP}}(\mathbf{y}_{1},\mathbf{y}_{2})=\frac{\sqrt{\int|F_{\mathbf{y}_{1}}(\mathbf{k})-F_{\mathbf{y}_{2}}(\mathbf{k})|^{2}d\mathbf{k}}}{\sqrt{\int|F_{\mathbf{y}_{1}}(\mathbf{k})|^{2}d\mathbf{k}}+\sqrt{\int|F_{\mathbf{y}_{2}}(\mathbf{k})|^{2}d\mathbf{k}}}\in[0,1]$
(1)

where $F_{\mathbf{y}}$ denotes the fast Fourier transform (FFT) of the signal
$\mathbf{y}(\mathbf{x},t)$. $\mathbf{k}$ denotes the wave number vector
$\mathbf{k}=\begin{bmatrix}k_{x}\end{bmatrix}$ for one-dimensional or
$\mathbf{k}=\begin{bmatrix}k_{x},k_{y}\end{bmatrix}^{\intercal}$ for two-
dimensional spatial inputs, respectively. The SSP is based on the Sobolev norm
of the individual signals
$||\mathbf{y}_{i}||^{2}=\int|F_{\mathbf{y}_{i}}(\mathbf{k})|^{2}d\mathbf{k}$
and their difference in space domain
$||\mathbf{y}_{1}-\mathbf{y}_{2}||^{2}=\int|F_{\mathbf{y}_{1}}(\mathbf{k})-F_{\mathbf{y}_{2}}(\mathbf{k})|^{2}d\mathbf{k}$,
where the definition of the Sobolev norm implies signed FFTs. By effectively
measuring the similarity of two signals in Fourier space, the SSP inherently
penalizes amplitude-, frequency- and phase errors at the same time and in a
single quantity. If both signals are identical, the Sobolev norm of their
difference is 0, such that $\mathrm{SSP}=0$ denotes perfect agreement among
the signals. With increasing difference between $\mathbf{y}_{1}$ and
$\mathbf{y}_{2}$, the SSP approaches 1. It is noteworthy that $\mathrm{SSP}=1$
(Perlin & Bustamante [2014] refer to $\mathrm{SSP}=1$ as perfect disagreement)
is only achieved for two phase-inverted signals disregarding their respective
amplitudes. The SSP is applied in the field of ocean engineering to quantify
numeric ocean wave prediction [Klein et al., 2020] and -reconstruction methods
[Desmars et al., 2021] as well as for ensuring consistent conditions in ship
movement experiments [van Essen, 2021]. Other applications include the
quantification of denoised GPS signals [Sharie et al., 2020] and demodulated
fringe pattern [Hernandez-Lopez et al., 2021]. Despite occasional usage in a
variety of fields, no publication reports an in-depth analysis on why the SSP
is preferred over established Euclidean distance metrics in those fields.
Furthermore, and to the knowledge of the authors, the SSP has not been used as
a loss metric for machine learning-related research yet. The following
paragraphs will highlight conceptual differences of the SSP and Euclidean
error metrics, and point out aspects that render the SSP a preferable loss
function for machine learning purposes.

###  2.2 Saturation, or: is there a limit to dissimilarity?

As the loss function shapes the gradients seen by the machine learning
optimizer, the training behavior, such as convergence speed, is directly
affected by the choice of the loss function. Steep gradients are generally
beneficial for fast convergence, but also flat gradient surfaces can allow the
optimizer to take larger steps and thereby help the training process. Focusing
at the predicted model output and the ground truth for long sequential data, a
loss function is all about a condensation of many pieces of information into a
scalar and abstract deviation quantity. In simple words, the loss function
tells how similar and how dissimilar can two sequences be? The former case is
trivial to answer by almost any Euclidean distance metric. The latter,
however, is nontrivial: what is the most dissimilar signal $y(t)$ to a given
non-zero signal $x(t)$? Euclidean metrics would answer that question by
increasing the amplitudes of $y(t)$ towards infinity. Alternatively, one could
state that $y(t)=-x(t)$ is the most dissimilar signal to $x(t)$, and what
about $y(t)=0$? Looking at Euclidean distance metrics, we can observe that
they are un-bounded, and hence over-sensitive to amplitude errors and under-
sensitive to phase or frequency errors if we consider signals composed of
harmonics. On the contrary, the SSP is a bounded quantity, which is most
conveniently shown by calculating the SSP for amplitude scaling signals, i.e.
$\mathbf{y}_{1}$ and
$\mathbf{y}_{2}=\kappa\cdot\mathbf{y}_{1},\;\kappa\in\mathbb{R}$. For this
case, Eq. (1) can be simplified to the $C_{0}$ continuous function

$J_{\mathrm{SSP}}(\mathbf{y}_{1},\kappa\cdot\mathbf{y}_{1})=\frac{|1-\kappa|}{1+|\kappa|}\begin{cases}\in[0,1)&\forall\kappa\in\mathbb{R}^{+}\\\
=1&\forall\kappa\in\mathbb{R}_{0}^{-}\end{cases},\quad\frac{\mathrm{d}^{2}J_{\mathrm{SSP}}}{\mathrm{d}\kappa^{2}}=\begin{cases}>0\;(\mathrm{convex})&\forall\kappa\in[0,1)\\\
<0\;(\mathrm{concave})&\forall\kappa>1\end{cases}$ (2)

i.e. we observe saturation behavior $\mathrm{SSP}=1$ for
$\kappa\in\mathbb{R}_{0}^{-}$, meaning that there is a bound to the
dissimilarity of two signals when measured in SSP. This behavior will flatten
the error functional in any considerable distance from a (local) minimum.

###  2.3 Comparison against Euclidean distance metrics

Figure 1: Evaluation of three different distance metrics under amplitude,
phase and frequency variation of the sample signal
$x(t)=A\cdot(0.5\sin(0.5t+\phi)\cdot(0.5\cos((2+3t)t)+0.5\cdot\cos(\omega t))$
(gray lines) compared against the baseline signal (solid black line in top
panel) defined by $A=1$, $\phi=\pi$, $\omega=15$. Error saturation sets in
only for the SSP metric

Figure 1 displays the MSE and MAE metrics in comparison to the SSP for a
synthetic test signal composed of three arbitrary harmonic components
including a frequency drift, a variable phase term $\phi$, a variable
frequency term $\omega$ as well as an amplitude scaling envelope term $A$.
Looking at the phase and frequency variation about the reference
configuration, there is no qualitative difference among the metrics,
particularly not in the gradients of the metrics in vicinity of the optimum.
In contrast, this observation does not hold for the amplitude variation. Per
definition, the MAE exhibits a constant gradient and the MSE shows a linearly
decreasing gradient, i.e. becoming very flat, towards the optimum. The SSP is
convex to the left, and concave to the right of the optimum according to Eq.
(2). While the Euclidean metrics increase monotonously for positive and
negative deviations from the reference amplitude, the SSP shows quick
saturation for zero and negative amplitudes, as well as for larger positives
amplitudes. In this case, there exists a large and perfectly flat error
functional for negative amplitudes $A$, hence forcing the optimizer to take
large steps and escape that desert of zero gradients. Once the optimizer
approaches the ground truth solution, suddenly gradients become steep towards
the minimum. Similarly, if the true solution is approached from far right,
i.e. positive $A$ values, gradients become steeper as the minimum is
approached. This is not the case for the convex MSE, which becomes flat
towards the minimum. In the upcoming section, we will find that the saturation
property of the SSP and the gradients in vicinity of the minimum will greatly
improve model convergence when using the SSP as loss function.

##  3 Results

We study the SSP as novel machine learning metric in regards of three central
aspects: i) the property to shape gradients (convergence), ii) the
regularization effect of the loss function (smoothness and final loss), and
iii) compatibility with conventional optimizers (instability) such as Adam
[Kingma & Ba, 2014] and SGD [Robbins & Monro, 1951, Kiefer & Wolfowitz, 1952].

The SSP is chosen as the baseline metric for assessing the similarity between
two signals in the validation data set. All validation and prediction errors
presented hereafter hence are computed in terms of SSP, even if the individual
models were trained towards a different loss function. We are well aware that
this choice can seem biased or doubtful to the reader. For this reason only,
we will consult the well-established MSE metric for reference alongside the
SSP as a validation set metric. Under outermost neutrality, we are convinced
that the SSP is superior to Euclidean metrics for comparing two signals that
are smooth and oscillatory following the reasoning presented before.

For the numerical experiments, we evaluate the performance of a fully
convolutional encoder-decoder model learning a flow map, i.e. taking a number
of states at preceding time steps and making a prediction about the system
state at the next time step. As typical test cases, we study the chaotic
dynamics of the (a) one- and (b) two-dimensional spatio-temporal dynamics in
the Kuramoto–Sivashinsky equation (KS-eq.) and the (c) time evolution of
dispersive surface waves, see Figure 2. Cases (a) and (b) have been studied as
benchmark cases in many recent works on data-driven forecasting methods
[Pathak et al., 2017, 2018, Raissi & Karniadakis, 2018, Vlachas et al., 2020].
Details on the neural architecture and the data generation are given in
Appendix 4.2, 4.3 and 4.4, respectively. Importantly, we note that our aim is
_not_ to propose an optimal neural architecture for future state forecasting,
but to study the training behavior of a given model architecture under
different loss functions.

Figure 2: Numerical test cases: (a) one-dimensional KS-eq., (b) two snapshots
($t=\\{10,100\\}\,$s) of the two-dimensional KS-eq. and (c) one-dimensional
broad-banded surface waves

The training behavior of the given neural architecture is studied under
different choices of the loss function, namely Mean Squared Error (MSE), Mean
Absolute Error (MAE), Root Mean Squared Error (RMSE), Huber loss and SSP. To
cross-validate our results and analyze their variability, each loss function
is used to train the model on 50 random (but equal across the loss functions)
data splits. The smaller the variance in the validation set scores, the more
robust a model is against weight initialization and training data
distribution. The validation set predictions of a single model are compared to
the ground truth data in terms of SSP and MSE metrics. We introduce the
convention of stating the model training loss function as a subscript of the
evaluation metric, e.g. $\mathrm{SSP}_{\mathrm{MSE}}$ indicates the validation
set error measured in terms of SSP for a model trained with MSE loss function.
As the loss function shapes the gradient used for training, a model trained on
the MSE loss is, naively, expected to minimize the MSE validation loss, while
a model trained on the SSP loss is expected to arrive at the smallest SSP
validation loss.

###  3.1 Chaotic dynamics in the 1D KS equation

The training behavior of machine learning models is usually evaluated using
the loss curve, which provides information about the final loss value a model
arrives at, the speed of convergence and the overall data fitting. We present
and discuss exemplary validation loss curves for all loss functions in the
following. Plotting the average loss curves over all 50 training runs would
smooth out certain characteristics, hence Figure 3 shows a representative
validation loss curve per loss function and evaluation metric. Particularly,
loss curves are shown for training runs on all loss functions on the exact
same data split. Model performance on the validation set is evaluated in terms
of the SSP (Figure 3 (a)) and MSE metric (b) after each epoch. Box plots in
Figure 3 depict summary statistics derived from all $50$ training runs per
loss function. For all loss functions and all model training runs there was no
overfitting observable.

Figure 3: Model training on data of the 1D KS-eq. ($\nu=0.25$, $\Delta t=0.1$
s) evaluated per (a) SSP metric and (b) MSE metric. 50 training runs on random
data splits are performed for each model. The loss curves show the validation
loss of one exemplary but representative training run. The vertical box-and-
whisker plots show the distribution of the minimum validation loss over the
last 25 epochs across all 50 training runs. The horizontal box-and-whisker
plots above each loss curve show the distribution of the epoch count at which
the loss falls below the threshold $\mathrm{SSP}=0.01$ ($\mathrm{MSE}=0.01$,
indicated by the horizontal dashed line)

To assess the effect of the different loss functions on the training behavior,
we first discuss the convergence speed. The horizontal box plots depict the
number of training epochs that was required to reach a certain loss threshold,
e.g. $\text{SSP}=0.01$ and $\text{MSE}=0.01$, respectively. On average, models
trained with SSP loss reached the threshold within $30$ epochs, earliest at
$17$ and latest after $60$ epochs. Models trained on RMSE loss show comparable
convergence speed, although the mean number of epochs is $41$. Huber loss and
MSE loss show rather similar statistics with median of $53$ and $64$ epochs,
respectively. The MAE loss shows a clearly different result with a very large
variation and $90$ epochs required on average to achieve $\text{SSP}=0.01$.
The variance in the model convergence behavior is minimal for the SSP and
maximal for the MAE: The MAE shows the largest interquartile range (IQR) of 48
epochs. MSE and Huber both have an IQR of 15 epochs. RMSE (10 epochs) and SSP
(9 epochs) have the lowest IQRs. As non-overlapping IQRs allow for statements
that are statistically significant, this experiment indicates the superiority
of the SSP loss function in terms of convergence speed. The qualitatively same
picture is resembled when evaluating the validation set predictions in terms
of MSE instead of SSP, see Figure 3 (b). Notably, models trained on SSP loss
will minimize the MSE-rated prediction error faster than models trained on MSE
loss. This observation will be confirmed in the following experiments on
different data sets. Therefore, the SSP loss function is shown to extend the
features taken into account by Euclidean loss functions, specifically the
phase and frequency information, which was exemplified in Section 2.3.

Next, we focus on the final validation loss, i.e. the performance of several
models trained for 250 epochs on the same data set but given a different loss
function. The vertical box-and-whisker plots in Figure 3 show the distribution
of the minimum loss value within epochs 225 and 250 across all training runs
per loss function. The model trained with SSP loss arrives at a median
validation loss of $\mathrm{SSP}_{\mathrm{SSP}}=0.0044$, and the model trained
on MSE loss arrives at $\mathrm{SSP}_{\mathrm{MSE}}=0.0052$. The Huber loss
performs roughly the same as the MSE, the RMSE arrives at a slightly lower
final loss value. Lastly, the MAE shows the weakest performance at a median of
$\mathrm{SSP}_{\mathrm{MAE}}=0.0061$. Considering the variance among all
training runs, the MAE has the largest IQR of $0.0012$. The remaining loss
functions exhibit similar IQRs at $0.0004$ (SSP and RMSE) and $0.0005$ (MSE
and Huber). As the IQR of the SSP does not overlap with any of the other loss
function’s IQRs, one can state that the superiority of the SSP in terms of
lowest final loss is statistically significant. Lastly, one can observe that
the curves fluctuate for all loss functions. However, in direct comparison,
the SSP loss function exhibits the smoothest loss characteristic. The loss
curves allow to glance at the optimizer paths towards a minimum: The less
fluctuating loss curve for the SSP thus indicates that the optimizer has less
trouble finding it’s path downwards on the loss surface compared to the other
loss functions.

To conclude, for the presented case, using the SSP loss function results in
smoother loss curves which point towards better regularization properties than
the competing loss functions. The SSP loss function not only speeds up the
training but also results in better validation scores after a pre-defined
number of training epochs. Training a model with SSP loss function results in
lower MSE validation errors as when training a model with MSE loss. This
behavior indicates the regularization property, i.e. showing that the SSP is
penalizing small prediction errors more strongly than other loss functions.

###  3.2 Gradient shaping through SSP loss function

The training behavior discussed in the previous paragraph indicates several
advantages of the SSP loss over other loss functions for the given prediction
task. In order to understand some of those observations, and to come back to
the discussion in Section 2.3, Figure 4 displays the loss hyperplane of a
single model (1D KS-eq. $\nu=0.25$, $\Delta t=1.0\,\mathrm{s}$, SSP loss) at
different stages of the training process, namely after epochs 1, 5, 10, and
500. For this analysis, two weights in the last layer, $\theta_{1}^{*}$ and
$\theta_{2}^{*}$, are varied while all other weights are kept frozen. The
prediction of that model for a single data sample is performed for a wide
variation of $\theta_{1},\theta_{2}$ about the final weights
$\theta_{1}^{*}=-1.07$, $\theta_{2}^{*}=0.86$ reached after 500 epochs of
training. For every choice of the weight values, the model prediction is
measured in terms of SSP, MSE, and MAE metric. As a result, one can study the
loss hyperplane and its gradients that the optimizer would see for that single
data sample if only two weights of the model were trainable. Even though this
situation does not reflect the actual training process, the error surfaces
depict salient differences that can link to the observations made before on
the training behavior of the SSP loss. The error isolines of the MAE are
separated by a constant value, indicating a linear gradient towards the
minimum. The MSE shows a different behavior, which resembles the shape of a
flat bowl as gradients get smaller towards the optimum. The SSP error surface
is not as circular, and exhibits large values particularly in the upper left
parameter region. Gradients are very large in the anti-diagonal direction of
the parameter space. The most important difference in the three metrics lies
in the area surrounding the minimum at $(\theta_{1}^{*},\theta_{2}^{*})$:
while the MSE would allow for a significant change in the weights without
crossing the lowest isoline, the SSP would increase very quickly under small
variations of the weights. For example, reducing $\theta_{1}^{*}$ by $10\%$
will increase the MSE by $0.3\%$, the MAE by $2.1\%$, and the SSP by $56.5\%$,
in relation to the full range of the respective hyperplane, for the model
after epoch 10. While the Euclidean metrics may accept a small variation in
that weight, the SSP would respond with a severe increase of the loss value.
Concluding, the local minima in the SSP error surface are very sharp and
narrow, therefore the training process may converge faster and more precisely
to the minimum, as the SSP is not as lenient to weight changes as the MAE and
MSE losses are.

Figure 4: Visualization of loss hyperplanes for a single data sample by
variation of two weights $\theta_{1}$, $\theta_{2}$ in the last layer of the
model about the final weights $\theta_{1}^{*}$, $\theta_{2}^{*}$ reached after
500 epochs. Loss hyperplanes are displayed for four stages of the training,
and model predictions are evaluated in terms of SSP (top row), MSE (middle
row) and MAE (bottom row) for the 1D KS-eq. ($\nu=0.25$, $\Delta
t=1.0\,\mathrm{s}$). The overall training process is depicted in the top panel

###  3.3 Effect of prediction task complexity

This section displays several experiments on the one- and two-dimensional KS-
eq. for different values of $\nu$, i.e. different levels of chaos in the
physical problem. Further, the time step $\Delta t$ is varied, which denotes
the time difference between consecutive data samples used as input to the flow
map model. Hence, for smaller $\nu$ and for larger $\Delta t$ the learning
task becomes more demanding, as more chaotic dynamics are to be predicted over
longer time instants into the future. The role of small perturbations and
inaccuracies is increased according to the exponential divergence driven by
positive Lyapunov spectra. Table 1 displays validation scores for five KS-eq.
cases of increasing complexity. For each case, models are trained on five
different loss functions, and the validation set scores are evaluated in terms
of the SSP and MSE metric. Cross-validation is performed against 50
realizations (10 for the 2D cases) of each experiment.

Table 1: Validation scores for all case studies, sorted by the estimated complexity of the prediction task from top (low complexity) to bottom (high complexity). For a given model and the same data set, training is performed for 250 epochs (100 for KS-2D cases) using different loss functions. In each case, the model prediction quality on the validation set is measured by metrics SSP and MSE. 50 training runs (10 for KS-2D) are performed for each set-up on different data splits, and mean $\pm$ standard deviation of validation scores are reported. Bold font and green backgrounds highlight the best score per metric and case case | loss | metric SSP | metric MSE
---|---|---|---
KS-1D ($\nu=1.0,\;\Delta t=0.1\,\mathrm{s}$) | SSP | $0.0034\pm 0.0008$ | $0.0006\pm 0.0004$
| MSE | $0.0023\pm 0.0003$ | $0.0002\pm 0.0001$
| MAE | $0.0045\pm 0.0012$ | $0.0016\pm 0.0030$
| RMSE | $0.0030\pm 0.0006$ | $0.0005\pm 0.0004$
| Huber | $\mathbf{0.0022\pm 0.0002}$ | $\mathbf{0.0002\pm 0.0000}$
KS-1D ($\nu=0.25,\;\Delta t=0.1\,\mathrm{s}$) | SSP | $\mathbf{0.0044\pm 0.0003}$ | $\mathbf{0.0014\pm 0.0002}$
| MSE | $0.0053\pm 0.0004$ | $0.0020\pm 0.0004$
| MAE | $0.0063\pm 0.0010$ | $0.0034\pm 0.0021$
| RMSE | $0.0049\pm 0.0004$ | $0.0018\pm 0.0006$
| Huber | $0.0052\pm 0.0004$ | $0.0018\pm 0.0003$
KS-1D ($\nu=0.25,\;\Delta t=1.0\,\mathrm{s}$) | SSP | $\mathbf{0.0755\pm 0.0020}$ | $\mathbf{0.4554\pm 0.0210}$
| MSE | $0.0800\pm 0.0023$ | $0.4882\pm 0.0262$
| MAE | $0.0812\pm 0.0021$ | $0.5039\pm 0.0253$
| RMSE | $0.0781\pm 0.0023$ | $0.4728\pm 0.0247$
| Huber | $0.0806\pm 0.0023$ | $0.4958\pm 0.0243$
KS-2D ($\nu=0.25,\;\Delta t=0.1\,\mathrm{s}$) | SSP | $\mathbf{0.0074\pm 0.0003}$ | $\mathbf{0.0043\pm 0.0004}$
| MSE | $0.0090\pm 0.0004$ | $0.0061\pm 0.0006$
| MAE | $0.0091\pm 0.0005$ | $0.0065\pm 0.0008$
| RMSE | $0.0090\pm 0.0011$ | $0.0064\pm 0.0018$
| Huber | $0.0089\pm 0.0008$ | $0.0064\pm 0.0019$
KS-2D ($\nu=0.25,\;\Delta t=1.0\,\mathrm{s}$) | SSP | $\mathbf{0.2336\pm 0.0044}$ | $\mathbf{4.2762\pm 0.0833}$
| MSE | $0.2483\pm 0.0033$ | $4.2998\pm 0.0619$
| MAE | $0.2473\pm 0.0025$ | $4.3478\pm 0.0347$
| RMSE | $0.2484\pm 0.0035$ | $4.3055\pm 0.0848$
| Huber | $0.2470\pm 0.0028$ | $4.3220\pm 0.0743$
dispersive waves ($\Delta t=1.3\,\mathrm{s}$) | SSP | $\mathbf{0.0162\pm 0.0015}$ | $\mathbf{0.0005\pm 0.0001}$
| MSE | $0.0210\pm 0.0018$ | $0.0008\pm 0.0001$
| MAE | $0.0242\pm 0.0016$ | $0.0010\pm 0.0001$
| RMSE | $0.0190\pm 0.0016$ | $0.0007\pm 0.0001$
| Huber | $0.0209\pm 0.0017$ | $0.0008\pm 0.0001$

The one-dimensional KS-eq. case using $\nu=1.0$, $\Delta t=0.1$ s represents
the seemingly simplest prediction task. Here, models trained on Huber loss
outperform models trained on any other loss function after the same amount of
training epochs. Huber loss reduces both validation metrics the most, and also
across all 50 different models trained on different splits of the data set.
The SSP loss achieves weaker results, while the MSE loss achieves only
slightly weaker results than the Huber loss. Moreover, the range of validation
scores in the cross validation study is minimal for the Huber loss models.
With increasing prediction task complexity ($\nu=0.25$, $\Delta
t=\\{0.1,1.0\\}\,\mathrm{s}$), we observe an expected degradation in overall
prediction quality for all models. However, both experiments show a similar
result: models trained on the SSP loss achieve the best validation scores,
even in terms of MSE. That is, if a model is trained using SSP loss, it
outperforms a model trained on MSE, even if the validation score is measured
in terms of MSE. Moreover, the standard deviation of validation scores across
all 50 models trained on SSP loss is the smallest for both experiments and
both validation metrics across all loss functions. For the two-dimensional KS-
eq. cases, a diverse observation is made: using a time increment of $\Delta
t=0.1$ s, no performance degradation is visible, and even an increase in
prediction accuracy can be read from the SSP metric when compared to the one-
dimensional case studies. However, as the time increment is increased, the
prediction task is rendered much more demanding. As a result, considerable
performance degradation is observed for $\Delta t=1.0\,\mathrm{s}$ throughout
all loss functions. Still, the SSP loss function manages to minimize the loss
in terms of both SSP and MSE metric best. The remaining loss functions perform
on par with each other, showing only minor deviations of
$\Delta_{\mathrm{SSP}}\leq 1.3\times 10^{-3}$ or $\Delta_{\mathrm{MSE}}\leq
4.8\times 10^{-2}$ between results.

The experiments conducted on the Kuramoto-Sivashinsky equation indicate that
the SSP is especially well-suited for more complex prediction tasks on
sequential data with an oscillatory character. While the Huber loss performs
best on the least complex problem, the models trained with SSP loss achieve
better validation scores for the remaining problems. The 2D implementation of
the SSP manages to pull ahead from the competition. Most importantly, the SSP
loss function manages to reduce the loss in terms of MSE better than the MSE
itself for all cases except the least complex one. While the presented results
are achieved using the Adam optimizer [Kingma & Ba, 2014], we ran a subset of
the experiments listed in Table 1 using the SGD optimizer [Robbins & Monro,
1951, Kiefer & Wolfowitz, 1952]. The loss curves for the SSP showed a very
smooth convergence compared to the remaining loss functions, however training
turned out to be unstable in some cases across all loss functions. Overall,
the results on SGD, while seeming less stable, confirm the observations
presented here.

###  3.4 Dispersive wave propagation

The third evaluation case deals with one-dimensional dispersive surface waves
governed by the linear wave equation Eq. (5). Compared to the one- and two-
dimensional nonlinear KS-eq., the surface waves are characterized by a very
broad wave spectrum, rendering this prediction case conceptually different to
the KS-eq. experiments. The dispersion relation lets long waves travel faster
than short waves, hence creating complex spatio-temporal patterns. Details on
the data generation are presented in Appendix 4.4. In analogy to the 1D KS-
eq., we cross validate our results against 50 training runs on independent and
random data splits. Again, each model is trained for 250 epochs. Figure 5 (a)
shows three exemplary model trainings using MAE, MSE and SSP as loss function.
Here, we can observe that the loss curves for MAE and MSE overlap for the most
part while the model with SSP loss undershoots both MAE and MSE for all
epochs. The validation scores (see Table 1) indicate that models with SSP loss
outperform models using different loss functions on the dispersive wave data.
On average, a mean validation set score of
$\mathrm{SSP}_{\mathrm{SSP}}=0.0162$ is achieved. The models trained on RMSE
performs second best ($\mathrm{SSP}_{\mathrm{RMSE}}=0.0190$), MSE and Huber
loss perform equally well at a mean $\mathrm{SSP}\approx 0.0210$ and lastly
the MAE performs the weakest ($\mathrm{SSP}_{\mathrm{MAE}}=0.0242$).

To quantify the convergence speed of each loss function, we define a threshold
at $\mathrm{SSP}=0.03$, indicated by the dashed line in Figure 5 (a). On
average over all 50 training runs, models with MAE loss reach the threshold
within 142 epochs, and models with MSE take 130 epochs. Models trained on RMSE
and Huber loss take 108 and 120 epochs, respectively. The models trained with
SSP loss converge considerably faster and reach the threshold after 78 epochs.
Training on Huber and SSP loss shows the smallest variability among the loss
functions with an IQR of 16 and 17 epochs, respectively. MSE and RMSE model
trainings both have an IQR of 24 epochs while the MAE loss exhibits the
largest variability with an IQR of 34 epochs. Variability is directly linked
with the sensitivity of the training process against weight initialization and
data distribution in the training sets. The low training variability in
combination with the highest convergence speed makes the SSP loss function
superior to the competition in this experiment. Again, the difference in
convergence speed can be stated statistically significant based on the non-
overlapping IQRs.

Figure 5: Model training on the 1D dispersive wave data evaluated in terms of
SSP metric. The loss curves (a) show the validation loss for models trained
with SSP, MSE and MAE loss functions. The gray area marks the first ten
epochs, where we obtain predictions for one validation sample after each
weight update. (b) shows the integral of the Fourier spectra of the respective
prediction after each weight update. The integral of the ground truth is
indicated by the dashed line. Panels (c) and (d) show the Fourier spectra of
the predictions using the model trained with MAE, MSE and SSP loss functions
(from top to bottom) after 2 (c) and 250 (d) training epochs, respectively.
The Fourier spectrum of each prediction is evaluated against the ground truth
spectrum shown in black

To shed more light onto the training process, we use the model to predict one
validation sample after each weight update within the first 10 epochs of
training (using a batch size of 2048, we obtain 81 weight updates per epoch).
The results discussed hereafter are representative for a wide range of sample
signals that were analyzed. In addition, both RMSE and Huber, being Euclidean
distance-based metrics, exhibit the same characteristics as portrayed by MSE
and MAE and are therefore not presented explicitly. Figures 5 (c) and (d) show
the spectra of the predictions after 2 and 250 epochs of training,
respectively. After 2 epochs, the models trained on MSE and MAE achieve
validation scores around $\mathrm{SSP}\approx 0.17$ while the model trained on
SSP achieves a slightly lower validation score of
$\mathrm{SSP}_{\mathrm{SSP}}=0.1584$. The spectra indicate that this marginal
lead is due to the fact that the model trained with SSP already attempts to
fit higher frequency components while the spectra associated with MAE and MSE
remain mostly flat for $k>0.2\,\mathrm{m}^{-1}$. After 250 epochs all models
fit the ground truth spectrum well over all frequencies with a slight
advantage of the SSP over MSE and MAE. Figure 5 (b) shows the integral of the
Fourier transform of each prediction within the first 10 epochs compared
against the ground truth integral value. As the integral represents a measure
of overall energy contained in the predicted wave signal, this value is mostly
related with amplitude errors. While the models trained with MAE and MSE loss
functions slowly approach the ground truth value, the model trained with SSP
exhibits an overshooting behavior in the very first training epoch. Hereafter,
the SSP model quickly approaches the ground truth. These observations link
well with the analytical discussion of the SSP metric and the conceptual
analysis presented in Section 2.3. First, as conventional Euclidean error
metrics are amplitude-biased, the overall energy in the predictions is very
small during the first training epochs. To reduce the loss value, the
optimizer reduces the overall amplitudes first, before energy is slowly
increased to match the ground truth. On the contrary, the SSP is not
amplitude-biased, such that also a frequency or phase error can increase the
loss value. The above mentioned characteristic of the SSP to shape higher
frequency components during early training results from this property. Second,
the convergence speed towards the ground truth value is linked to the
gradients in vicinity of a loss minimum: the MSE gradient becomes more flat as
the minimum is approached, such that convergence becomes slower over the
number of training epochs. Even for excessive training, the MSE exhibits a
non-vanishing difference to the ground truth, which is fundamentally rooted in
the bowl-type shape of the MSE. The constant gradient of the MAE results in a
steady approach of the minimum. Conceptually different, the gradients of the
SSP are very sharp in vicinity of a minimum, hence increasing convergence
speed. Predictions far from the optimum enter the saturating loss behavior,
hence flattening out the loss plane in sub-optimal parameter regions. Overall,
the SSP is the most unforgiving in cases of small errors in amplitude, phase,
or frequency, hence enhancing the convergence behavior of the neural network
training. This supremacy becomes especially evident for the broad-banded wave
spectrum, while it is already observable for the chaotic KS-eq. case studies.
Thus, we conclude that the SSP is a particularly promising alternative to
conventional loss functions when training machine learning models on
sequential data with complex spectra.

##  4 CONCLUSION

This work proposes the surface similarity parameter (SSP) as a novel machine
learning loss function that is designed for oscillatory sequences arising
across many fields of research related to physics and natural processes. Our
experiments cover the chaotic dynamics in the one- and two-dimensional
Kuramoto-Sivashinsky equation and dispersive surface gravity waves in fluids.

Our research investigates the effect of the loss function on the training
behavior of a machine learning flow map, i.e. a fully convolutional neural
network that advances historic state observations one time increment into the
future. The experiments indicate that the proposed SSP loss function is
superior to classical Euclidean loss functions, such as MSE and MAE, in terms
of a) convergence rate, b) robustness against weight initialization effects,
and c) final model prediction accuracy. Compared to conventional Euclidean
loss functions, the SSP loss function is a bounded error metric that exhibits
steep gradients in vicinity of local minima. This regularization property
allows to faster train more precise models. While Euclidean distance-based
metrics turn out to be amplitude-biased for oscillatory data, the SSP metric
is taking amplitude, phase and frequency information explicitly into account,
and is therefore much stricter even for small deviations between prediction
and ground truth sequences. The superiority of the SSP is especially evident
the more complex the spectral composition of the sequential data is, i.e. the
wider the spectrum of harmonic components and the more irregular the amplitude
distribution becomes. The novel SSP loss function is therefore highly
promising both for regression and classification tasks on oscillatory
sequential data.

## APPENDIX

###  4.1 Conventional loss functions

Table 2 lists an overview on all conventional Euclidean loss functions used in
the course of this work, accompanied by a short description of the
interpretation.

Table 2: Overview on common regression loss functions for sequence prediction tasks loss |  definition |  comment
---|---|---
mean absolute error |  $\mathrm{MAE}=\frac{1}{N}\sum_{i}^{N}\left|y_{i}-\hat{y}_{i}\right|$ |  arithmetic mean of absolute errors, also denoted $L$1 loss. Estimates the median of the error distribution function, therefore robust against outliers. But: constant and large (!) gradients also for small errors, requires adaptive learning rates for better convergence. Gradients are not continuous (not differentiable at 0.)
mean squared error |  $\mathrm{MSE}=\frac{1}{N}\sum_{i}^{N}\left(y_{i}-\hat{y}_{i}\right)^{2}$ |  arithmetic mean of squared errors, also denoted $L$2 loss. Estimates the mean of the error distribution function, put emphasis on large values. Large gradients at large errors, small gradients at small errors. Helps training neural networks with constant learning rate
root mean squared error |  $\mathrm{RMSE}=\frac{1}{N}\sqrt{\sum_{i}^{N}\left(y_{i}-\hat{y}_{i}\right)^{2}}$ |  arithmetic mean of squared errors, also denoted $L$2 loss. Same scale as MAE. Estimates the mean of the error distribution function, put emphasis on large values.
Huber loss |  $\mathcal{L}_{\delta}=\begin{cases}\frac{1}{2}(y-\hat{y})^{2}&\mbox{for }|y-\hat{y}|\leq\delta,\\\ \delta(|y-\hat{y}|-\frac{1}{2}\delta)&\mbox{else}\end{cases}$ |  Combination of MSE (small errors) and MAE (large errors) with small and large gradients, respectively, and differentiable at 0.

###  4.2 Encoder-decoder model architecture

We use a fully convolutional encoder-decoder model for all learning tasks
conducted in this work. The model structure for the one-dimensional cases is
presented in Figure 6. The input is first processed by three sequential
encoder blocks, where each block extracts features and down-samples the
feature maps along the spatial dimension. While the spatial dimension of the
internal input representation halves, the amount of convolutional filters is
doubled per stage. In the bottleneck, only the most relevant features remain.
The decoder path is set up the same but uses up-sampling layers to recreate
the initial spatial dimension of the input within three sequential decoder
blocks. We apply a final convolutional layer with one filter of size 1 to
merge the information from eight channels into one signal
$\mathrm{y}\in\mathbb{R}^{n_{b}\times 1024}$ where $n_{b}$ indicates the batch
size. The model for the two-dimensional cases shares the same architecture but
uses the two-dimensional counterparts of the layers depicted in Figure 6.

Figure 6: Fully convolutional network with encoder-decoder structure for 1D
data. The numbers above each feature map indicate the number of channels, i.e.
the amount of applied filters in the preceding convolutional layer. All layers
apply padding to maintain spatial dimensions. A final flatten layer that
solely eliminates the channel dimension is not depicted here. The model maps
an input $\mathrm{X}\in\mathbb{R}^{n_{b}\times 1024\times 8}$ to an output
$\mathrm{y}\in\mathbb{R}^{n_{b}\times 1024}$ where $n_{b}$ indicates the batch
size

###  4.3 Data generation – Kuramoto-Sivashinsky equation

The Kuramoto-Sivashinsky equation has been studied in several recent works
[Smaoui & Al-Enezi, 2004, Pathak et al., 2017, 2018, Raissi & Karniadakis,
2018, Vlachas et al., 2018, Bar-Sinai et al., 2019] related to machine
learning as illustrative example for complex spatio-temporal dynamics. For our
case studies, we use the KS-eq.

$u_{t}+\Delta u+\nu\Delta^{2}u+\frac{1}{2}|\nabla u|^{2}$ (3)

augmented with a viscosity term $\nu\in[0.25,1]$ that controls the amount of
chaos in the system [Smaoui & Al-Enezi, 2004, Vlachas et al., 2018]. We
initialize both our 1D and 2D KS systems with the superposition of 10 waves

$u(x)=\sum\limits_{i=1}^{10}A_{i}\sin(2\pi\frac{i}{L}x+\varphi_{i}),\quad
u(x,y)=\sum\limits_{i=1}^{10}A_{i}\sin(2\pi\frac{i}{L}\left(x\cos(\phi_{i})+y\sin(\phi_{i})\right)+\varphi_{i})$
(4)

with random amplitude $A_{i}\in[-1,1]$ and random phase
$\varphi_{i}\in[0,2\pi]$, $\phi_{i}\in[-\pi/2,\pi/2]$. For the 1D KS equation,
we use a domain $x$ with length $L_{x}=128\,\mathrm{m}$ discretized into 1024
grid points. The 2D KS equation uses a square domain with the same set-up as
the 1D case in $x$ and $y$ dimension. The definition of our initial conditions
Eq. (4) ensure that the longest wave fits the domain exactly. All derivatives
in Eq. (3) are calculated using the 4th order finite differences method. To
ensure numeric stability, we apply a low pass filter so that only the first 80
(2D: 60) frequency components are taken into account after each time step. For
time integration of the PDE we use the modified Euler method [Süli & Mayers,
2003, Burden & Faires, 2010]. The system is simulated using $\Delta
t_{\mathrm{sim}}=0.01\,\mathrm{s}$ of which every 10th system state is saved
to disk, resulting in a data time step size of $\Delta t=0.1\,\mathrm{s}$. The
KS-eq. (3) is simulated 100 times (2D: 10 times) for $200\,\mathrm{s}$ for
$\nu=\\{0.25,1.0\\}$. The simulation data is preprocessed using a sliding
window approach with stride 1 and shift 1 or 10 for $\Delta t=0.1\,\mathrm{s}$
or $\Delta t=1.0\,\mathrm{s}$, respectively. For $\Delta t=0.1\,\mathrm{s}$ we
obtain 1993, for $\Delta t=1.0\,\mathrm{s}$ we end up with 1921 input-output-
pairs $(\mathbf{X}\in\mathbb{R}^{1024\times
8},\;\mathbf{y}\in\mathbb{R}^{1024})$ per simulation (2D: 1993 input-output
pairs $(\mathbf{X}\in\mathbb{R}^{1024\times 1024\times
8},\;\mathbf{y}\in\mathbb{R}^{1024\times 1024})$ for $\Delta
t=0.1\,\mathrm{s}$, 1921 for $\Delta t=1.0\,\mathrm{s}$).

###  4.4 Data generation – dispersive surface waves

The dispersive surface wave data is generated using linear wave theory, i.e.
no interactions between waves are modeled. A superposition of $N=651$ single
harmonic waves forms the initial surface wave

$\zeta_{0}=\zeta(x,t=0)=\sum\limits_{i=1}^{N}A_{i}\cos(k_{i}x+\varphi_{i})$
(5)

with amplitude $A_{i}$, wave number $k_{i}$ and random phase
$\varphi_{i}\in[0,2\pi]$. The amplitudes are defined by a normalized Fourier
spectrum [Klein et al., 2021]

$F(\omega_{i})=\frac{27\left(\omega_{i}-\omega_{\mathrm{min}}\right)\cdot\left(\omega_{i}-\omega_{\mathrm{max}}\right)^{2}}{4\left(\omega_{\mathrm{max}}-\omega_{\mathrm{min}}\right)^{3}}$
(6)

with $\omega_{\mathrm{min}}=0.3\,\mathrm{rad}\,\mathrm{s}^{-1}$ and
$\omega_{\mathrm{max}}=2.0\,\mathrm{rad}\,\mathrm{s}^{-1}$ limiting the
relevant frequency range. The wave number $k$ is linked to the angular
frequency $\omega$ via the linear dispersion relation

$\omega(k)=\sqrt{g\cdot k\cdot\tanh(k\cdot d)}$ (7)

with gravitational acceleration $g$ and water depth $d=500\,\mathrm{m}$. In
simple terms, the dispersion relation Eq. (7) defines the propagation velocity
of a wave depending on its wave length. In consequence, long waves travel
faster than short waves. To simulate the time evolution of the initial wave
Eq. (5), the assumption of linear wave theory allows to simulate the time
evolution of each single harmonic wave according to their respective
dispersion relation and recreate the surface wave at any given time by
superposition of all single harmonic waves. Since $\zeta_{0}$ carries all
information about the space domain (wave numbers $\mathbf{k}\in\mathbb{R}^{N}$
and initial phases $\boldsymbol{\varphi}\in\mathbb{R}^{N}$), we can simulate
the time evolution of $\zeta_{0}$ in the complex space by application of the
Fourier transform

$F_{\zeta(t)}=F_{\zeta_{0}}\exp(-i\boldsymbol{\omega}t)$ (8)

and obtain the surface wave $\zeta(t)$ by means of the inverse Fourier
transform of $F_{\zeta(t)}$. As this FFT-based simulation approach introduces
periodic boundary conditions, the simulated wave exists the domain to the
right and re-enters it on the left. To counteract this behavior, we append
zeros to $\zeta_{0}$ thus expanding the domain to the right. The initial wave
now travels into the appended zeros instead of re-entering the domain.
Overall, we generate 100 initial sea states with random phase $\varphi_{i}$
and otherwise identical configuration. Every sea state is simulated for
$250\,\mathrm{s}$ with $\Delta t=0.1\,\mathrm{s}$. The data preprocessing for
ML involves the application of a sliding window with stride 1 and shift 13,
increasing the time step in the final data set to $\Delta t=1.3\,\mathrm{s}$.
Using this approach, we obtain 2396 input-output-pairs
$(\mathbf{X}\in\mathbb{R}^{1024\times 8},\;\mathbf{y}\in\mathbb{R}^{1024})$
from each simulation.

## ACKNOWLEDGMENT

This work was supported by the Deutsche Forschungsgesellschaft (DFG – German
Research Foundation) [grant number 277972093]; and Hamburg University of
Technology (TUHH) [$\text{I}^{3}$-project “Predicting Ship Hydrodynamics to
Enable Autonomous Shipping: Nonlinear Physics and Machine Learning”].

## References

  * Bar-Sinai et al. [2019] Bar-Sinai, Y., Hoyer, S., Hickey, J., & Brenner, M. P. (2019).  Learning data-driven discretizations for partial differential equations.  Proceedings of the National Academy of Sciences, 116, 15344–15349. URL: https://doi.org/10.1073/pnas.1814058116. doi:10.1073/pnas.1814058116.
  * Burden & Faires [2010] Burden, R. L., & Faires, J. D. (2010).  Numerical Analysis.  Clifton Park, NY: Cengage Learning.
  * Desmars et al. [2021] Desmars, N., Hartmann, M., Behrendt, J., Klein, M., & Hoffmann, N. (2021).  Reconstruction of Ocean Surfaces From Randomly Distributed Measurements Using a Grid-Based Method.  In Volume 6: Ocean Engineering (p. V006T06A059).  Virtual, Online: American Society of Mechanical Engineers.  URL: https://asmedigitalcollection.asme.org/OMAE/proceedings/OMAE2021/85161/V006T06A059/1121470. doi:10.1115/OMAE2021-62409.
  * van Essen [2021] van Essen, S. (2021).  Influence of Wave Variability on Ship Response During Deterministically Repeated Seakeeping Tests at Forward Speed.  In T. Okada, K. Suzuki, & Y. Kawamura (Eds.), Practical Design of Ships and Other Floating Structures Lecture Notes in Civil Engineering (pp. 899–925).  Singapore: Springer.  doi:10.1007/978-981-15-4624-2_54.
  * Fulcher & Jones [2014] Fulcher, B. D., & Jones, N. S. (2014).  Highly comparative feature-based time-series classification.  IEEE Transactions on Knowledge and Data Engineering, 26, 3026–3037. doi:10.1109/TKDE.2014.2316504.
  * Hernandez-Lopez et al. [2021] Hernandez-Lopez, F. J., Legarda-Sáenz, R., & Brito-Loeza, C. (2021).  Parallel algorithm for fringe pattern demodulation.  Journal of Real-Time Image Processing, 18, 2441–2451. URL: https://link.springer.com/10.1007/s11554-021-01129-4. doi:10.1007/s11554-021-01129-4.
  * Kiefer & Wolfowitz [1952] Kiefer, J., & Wolfowitz, J. (1952).  Stochastic estimation of the maximum of a regression function.  The Annals of Mathematical Statistics, 23, 462–466. URL: https://doi.org/10.1214/aoms/1177729392. doi:10.1214/aoms/1177729392.
  * Kingma & Ba [2014] Kingma, D. P., & Ba, J. (2014).  Adam: A method for stochastic optimization.  arXiv:arXiv:1412.6980.
  * Klein et al. [2020] Klein, M., Dudek, M., Clauss, G. F., Ehlers, S., Behrendt, J., Hoffmann, N., & Onorato, M. (2020).  On the deterministic prediction of water waves.  Fluids, 5. URL: https://www.mdpi.com/2311-5521/5/1/9. doi:10.3390/fluids5010009.
  * Klein et al. [2021] Klein, M., Hartmann, M., & von Bock und Polach, F. (2021).  Note on the application of transient wave packets for wave–ice interaction experiments.  Water, 13. URL: https://www.mdpi.com/2073-4441/13/12/1699. doi:10.3390/w13121699.
  * Pathak et al. [2018] Pathak, J., Hunt, B., Girvan, M., Lu, Z., & Ott, E. (2018).  Model-free prediction of large spatiotemporally chaotic systems from data: A reservoir computing approach.  Physical review letters, 120, 024102. doi:10.1103/PhysRevLett.120.024102.
  * Pathak et al. [2017] Pathak, J., Lu, Z., Hunt, B. R., Girvan, M., & Ott, E. (2017).  Using machine learning to replicate chaotic attractors and calculate lyapunov exponents from data.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 27, 121102\. URL: https://doi.org/10.1063/1.5010300. doi:10.1063/1.5010300.
  * Perlin & Bustamante [2014] Perlin, M., & Bustamante, M. (2014).  A robust quantitative comparison criterion of two signals based on the sobolev norm of their difference.  Journal of Engineering Mathematics, 101, 115--124.
  * Raissi & Karniadakis [2018] Raissi, M., & Karniadakis, G. E. (2018).  Hidden physics models: Machine learning of nonlinear partial differential equations.  Journal of Computational Physics, 357, 125--141. doi:10.1016/j.jcp.2017.11.039.
  * Robbins & Monro [1951] Robbins, H., & Monro, S. (1951).  A stochastic approximation method.  The Annals of Mathematical Statistics, 22, 400--407. URL: https://doi.org/10.1214/aoms/1177729586. doi:10.1214/aoms/1177729586.
  * Sharie et al. [2020] Sharie, M., Mosavi, M. R., & Rahemi, N. (2020).  Determination of an appropriate mother wavelet for de-noising of weak GPS correlation signals based on similarity measurements.  Engineering Science and Technology, an International Journal, 23, 281--288. URL: https://www.sciencedirect.com/science/article/pii/S2215098618322018. doi:https://doi.org/10.1016/j.jestch.2019.05.006.
  * Smaoui & Al-Enezi [2004] Smaoui, N., & Al-Enezi, S. (2004).  Modelling the dynamics of nonlinear partial differential equations using neural networks.  Journal of Computational and Applied Mathematics, 170, 27--58. URL: https://doi.org/10.1016/j.cam.2003.12.045. doi:10.1016/j.cam.2003.12.045.
  * Süli & Mayers [2003] Süli, E., & Mayers, D. F. (2003).  An Introduction to Numerical Analysis.  Cambridge: Cambridge University Press.
  * Vlachas et al. [2020] Vlachas, P., Pathak, J., Hunt, B., Sapsis, T., Girvan, M., Ott, E., & Koumoutsakos, P. (2020).  Backpropagation algorithms and reservoir computing in recurrent neural networks for the forecasting of complex spatiotemporal dynamics.  Neural Networks, 126, 191--217. URL: https://doi.org/10.1016/j.neunet.2020.02.016. doi:10.1016/j.neunet.2020.02.016.
  * Vlachas et al. [2018] Vlachas, P. R., Byeon, W., Wan, Z. Y., Sapsis, T. P., & Koumoutsakos, P. (2018).  Data-driven forecasting of high-dimensional chaotic systems with long short-term memory networks.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 474, 20170844. URL: https://doi.org/10.1098/rspa.2017.0844. doi:10.1098/rspa.2017.0844.