# An End-to-End Vehicle Trajcetory Prediction Framework
††thanks: Funded by The Hong Kong Polytechnic University College of
Undergraduate Researchers & Innovators (PolyU CURI)’s Undergraduate Research &
Innovation Scheme (URIS)

Fuad Hasan  Department of Aeronautical & Aviation Engineering
The Hong Kong Polytechnic University
Kowloon, Hong Kong SAR
<EMAIL_ADDRESS>Hailong Huang  Department of Aeronautical&
Aviation Engineering
The Hong Kong Polytechnic University
Kowloon, Hong Kong SAR
<EMAIL_ADDRESS>
###### Abstract

Anticipating the motion of neighboring vehicles is crucial for autonomous
driving, especially on congested highways where even slight motion variations
can result in catastrophic collisions. An accurate prediction of a future
trajectory does not just rely on the previous trajectory, but also, more
importantly, a simulation of the complex interactions between other vehicles
nearby. Most state-of-the-art networks built to tackle the problem assume
readily available past trajectory points, hence lacking a full end-to-end
pipeline with direct video-to-output mechanism. In this article, we thus
propose a novel end-to-end architecture that takes raw video inputs and
outputs future trajectory predictions. It first extracts and tracks the 3D
location of the nearby vehicles via multi-head attention-based regression
networks as well as non-linear optimization. This provides the past trajectory
points which then feeds into the trajectory prediction algorithm consisting of
an attention-based LSTM encoder-decoder architecture, which allows it to model
the complicated interdependence between the vehicles and make an accurate
prediction of the future trajectory points of the surrounding vehicles. The
proposed model is evaluated on the large-scale BLVD dataset, and has also been
implemented on CARLA. The experimental results demonstrate that our approach
outperforms various state-of-the-art models.

###### Index Terms:

autonomous driving, vehicle trajectory prediction, multi-head attention, LSTM,
3D multi-object tracking

##  I Introduction

Recently, there has been a great deal of interest in autonomous driving
systems, as it promises safer commute on the road. Among the numerous
fundamental functions of autonomous driving, the prediction of the future
trajectory of nearby vehicles has received special attention. When driving on
a congested highway where cars frequently travel at high velocities, even a
minor error can cause a chain reaction that results in a serious accident. In
order to assess the danger of its own maneuvers, a self-driving vehicle must
be able to predict the future paths of nearby vehicles.

In practice, however, it is exceedingly difficult to accurately predict the
future trajectories of surrounding vehicles because it is dependent not only
on their historical trajectories but also on the socio-temporal
interdependence [1] of the dense web of vehicular traffic surrounding the car.
A model of the interaction between the participants must be used as one of the
prediction parameters in order to make accurate predictions.

Initially, the majority of deep learning models concentrated on extracting the
temporal correlation using RNN architectures such as LSTM and GRU [2, 3],
which only considered the historical trajectory to make predictions for the
future [4]. Due to their encoder-decoder sequence-to-sequence architecture,
RNNs were incapable of modeling social or spatial interactions. To address
this issue, recent research [5, 6] has explicitly concentrated on coupling
RNNs with social feature extraction architectures to model vehicle
interaction. Attention-based networks are a recent innovation due to their
capacity to swiftly extract vital data from historical data. Two attention
layers were stacked to extract both spatial and temporal attention in [7]
predicting the trajectory of pedestrians. Due to the problem of accumulating
errors resulting from the autoregressive decoding procedure of Transformers
[8], the multi-headed attention mechanism, originally proposed in the
traditional Transformer [9], has not been extensively explored in the highway
trajectory prediction problem despite these advancements.

In addition, the NGSIM dataset [10] has been widely used to evaluate the
validity of the vast majority of highway trajectory prediction studies. The
limitation of implementing the NGSIM dataset is the fact that it contains
readily available past trajectories. Most of the literature thus has
exclusively focused on future trajectory prediction based on these past
trajectories. To the knowledge of the authors, there has not been much work on
a more robust, end-to-end pipeline that can both extract the 3D information
from video sequences and also produce readily available past trajectories in
order to extract the socio-temporal features and consequently predict future
trajectories. An end-to-end pipeline makes it easier to deploy the solution on
an actual autonomous vehicle as well, negating the dependence on 3D pose
estimation architectures.

Figure 1: Overview of Proposed architecture

Therefore, in this research, we propose an end-to-end vehicle trajectory
prediction framework. We exploit the multi-head attention (MHA) mechanism of
the Transformer to regress the 3D information from raw images and then fuse it
with an LSTM sequence-to-sequence architecture which solves the decoder’s
autoregressive accumulative error problem.

The following is a summary of the contribution and academic novelty of our
original work:

  * •

We propose a novel end-to-end architecture in order to produce predictions
directly from video/image sequences.

  * •

We utilize the MHA mechanism to extract 3D information from images, encode the
trajectory in order to extract socio-temporal features and also extract future
socio-temporal features.

  * •

We evaluate our method on BLVD, a large-scale dataset of vehicle trajectories
as well as video sequences. The experimental results demonstrate that our
model strongly performs against state-of-the-art methodologies implemented on
both the NGSIM and BLVD datasets, despite using video sequences as inputs.

##  II Methodology

###  II-A Problem Formulation

In this article, the problem of vehicle trajectory prediction has been
formulated as a non-linear regression task where the inputs to the model are
arbitrary video clip with frames
$\left\\{\mathcal{I}_{1},\mathcal{I}_{2},...\mathcal{I}_{n}\right\\}$ where
$\mathcal{I}_{i}$ is a single image frame of the video clip. Based on the
video, the task is to first extract the ego-centric coordinates of each
vehicle $i$ in the frame as $p^{i}=(x^{i},y^{i})$. This can then be
concatenated in order to produce the past trajectories of each vehicle
represented as

$X=\left\\{p_{1},p_{2},\ldots,p_{T}\right\\},$ (1)

where

$p_{t}=\left\\{(x_{t}^{0},y_{t}^{0}),(x_{t}^{1},y_{t}^{1}),\ldots,(x_{t}^{N},y_{t}^{N})\right\\}.$
(2)

$p_{t}^{i}=(x_{t}^{i},y_{t}^{i})$ are the coordinates $x$ and $y$ of vehicle
$i$ at timestep $t$, where $t\in[1,T]$ and $i\in[1,N]$, $T$ stands for the
total length of the observed trajectory and $N$ is the total number of
observed neighbouring vehicles within 30 metres of the target vehicle, i.e.
$\left|p_{t}^{target}-p_{t}^{i}\right|\leq 30$.

Based on the past trajectories, the objective of the proposed model is to
learn a non-linear regression function that predicts the coordinates of all
observed vehicles in the prediction horizon $F$, such that the predicted
trajectories of the neighbouring vehicles are represented as $\hat{Y}$, as
follows, which approximates the ground truth (GT) trajectory,
$Y=\left\\{{q}_{T+1},{q}_{T+2},\ldots,{q}_{T+F}\right\\}$

$\hat{Y}=\left\\{\hat{q}_{T+1},\hat{q}_{T+2},\ldots,\hat{q}_{T+F}\right\\},$
(3)

where

$\hat{q}_{t}=\left\\{(\hat{x}_{t}^{0},\hat{y}_{t}^{0}),(\hat{x}_{t}^{1},\hat{y}_{t}^{1}),\ldots,(\hat{x}_{t}^{N},\hat{y}_{t}^{N})\right\\}.$
(4)

The frame of reference of the dataset used in this paper is by default ego-
centric as the trajectories were extracted from the onboard sensors. The
longitudinal y-axis indicates the direction forward to the ego vehicle and the
lateral x-axis direction indicates the axis perpendicular to it. The right-
hand side is considered positive according to the dataset [11]. As a result,
the model is independent of the curvature of the road which conveniently
allows it to be used in the highway as long as an object-detection and lane
estimation algorithm is built on the target vehicle [12].

###  II-B Network Overview

An overview of the model we propose is shown in Figure 1. It is divided into
the following three modules: Trajectory Processor, Encoder and Decoder. We
treat the surrounding vehicles as 3D bounding boxes. The Trajectory Processor
first extracts the location of the centroid of these bounding boxes. We
implemented DeepSORT[13] in order to solve the Multi-Object Tracking (MOT)
problem [14] first and get the 2D bounding boxes of surrounding vehicles. From
the 2D bounding box patch provided by DeepSORT, we regress the dimension and
orientations via two MHA layers to constraint and solve the optimization of
the 3D real-world, ego-centric position of the 3D bounding boxes. This is then
concatenated across frames to provide the historical tracks of the vehicles.
The encoder module then encodes the socio-temporal attention data into the
historical tracks via two distinct multi-head-attention layers and then also
encodes temporal memory via the LSTM layer. The decoder then decodes the
encoded information into future predictions. The hidden features are passed
onto successive decoding steps via further MHA layers to improve further
future predictions which also prevents the traditional transformer decoding
accumulation error by enriching the hidden features.

###  II-C Trajectory Processor

Figure 2: Skeleton of the IMHA layers

####  II-C1 Constraint Design

The 3D bounding box for surrounding vehicles is described by its
translation/position $T=(t_{x},t_{y},t_{z})$, dimensions
$D=(d_{x},d_{y},d_{z})$ and orientation $R=(\theta,\phi,\psi)$. The projection
of a 3D point $X=(X,Y,Z,1)^{T}$ into the image coordinate frame
$x=(x,y,1)^{T}$ with Camera Intrinsic $K$ is described as follows:

$\mathbf{x}=K\left[\begin{array}[]{ll}R&T\end{array}\right]\mathbf{X}$ (5)

Taking the assumption that the center of the object coincides with the center
of the 3D bounding box with known dimensions D, the vertex coordinates of the
3D bounding box is as follows:

$\mathbf{X_{1}}=(\frac{d_{x}}{2},\frac{d_{y}}{2},\frac{d_{z}}{2}),\ldots\mathbf{X_{8}}=(-\frac{d_{x}}{2},-\frac{d_{y}}{2},-\frac{d_{z}}{2})$
(6)

In line with existing literature [13, 15], we further assume that the
projection of the 3D bounding box fits tightly inside the 2D bounding box.
This assumption requires that every side of the 2D box must contain at least
one of the vertices of the projection of the 3D box. As an example, we can say
that vertex $\mathbf{X_{1}}=(\frac{d_{x}}{2},\frac{d_{y}}{2},\frac{d_{z}}{2})$
lies on the left side of the 2D box with coordinate $x_{max}$. This
correspondence of the point to the side of the 2D box can be represented in
the following equation:

$x_{\min}=\left(\left[\begin{array}[]{ll}R&T\end{array}\right]\left[\begin{array}[]{c}d_{x}/2\\\
-d_{y}/2\\\ d_{z}/2\\\ 1\end{array}\right]\right)_{x}$ (7)

Similarly, three other constraint equations can be derived for
$x_{min},y_{max},y_{min}$, totaling to four constraint equations, so far.
Further we make the assumption that the vehicles are upright, the relative
roll of the vehicles are zero, the pitch of the vehicles are also zero and the
driving scene is flat (highway), which reduces the possible configurations to
only 64. Therefore, regressing the dimensions and orientations enables us to
optimize the configurations in order to solve for translation $T$, from which
past trajectories can be derived.

####  II-C2 Image Multi-Head Attention

We implemented two Multi-Head Attention layers in the module Image Multi-Head
Attention (IMHA) in order to regress the Dimension $d_{x},d_{y},d_{z}$ and the
orientation $R=(\theta,\phi,\psi)$. As we assume roll $\phi$ and pitch $\psi$
to be zero, we only regress $\theta$ which is the yaw or azimuth angle of the
vehicles. However, the visual appearance of the images correlates strongly
with local orientation $\theta_{l}$ not the global orientation $\theta$, the
two being related according to Equation 8. $\theta_{ray}$, which is the
relative angle of the ray that intersects with a particular pixel, can be
easily computed via the camera intrinsics.

$\theta=\theta_{ray}+\theta_{l}$ (8)

The IMHA module thus directly outputs the dimension $d_{x},d_{y},d_{z}$ and
local orientation $\theta_{l}$ (which consequently gives $R$) of the bounding
boxes from the 2D image patch, which serves as inputs to Equation 5. For
example, let the vertical side of the 2D bounding box be $\mathbf{X^{i}}$,
which is the $i^{t}h$ corner of the 3D Box constrained in the 2D box, then the
constraint equation is the following:

$\left(K\left[\begin{array}[]{cc}I&R\times\mathbf{X}^{\mathbf{k}}\\\
0&1\end{array}\right]\left[\begin{array}[]{c}T_{x}\\\ T_{y}\\\ T_{z}\\\
1\end{array}\right]\right)_{x}=x_{\text{min}}$ (9)

where $I$ is a $3\times 3$ identity matrix and $i$ corresponds to possible
configurations of the 3D box corned projected on the 2D box as explained in 5.
The only unknowns remain to be the 3DoF relative position of the vehicles,
$T=\left[T_{x},T_{y},T_{z}\right]$. As it is an over-determined set of linear
equations with only one global minimum, it is solved by simple least square
fitting optimization with a zero initial guess. The position coordinates are
then concatenated to produce trajectory arrays.

###  II-D Encoder-Decoder

We have taken inspiration from [16] in order to design the Encoder-Decoder
module in this section. To input the trajectory data into our proposed model,
the coordinates are scaled first to a range of (-1, 1) to aid the model to
reach faster convergence [17]. This is then run through an embedding layer, a
multi-layer perceptron (MLP) which, as shown in Eq. 10, maps raw input
features to a higher dimension and outputs $X_{\lambda}$.

${\Lambda}_{t}^{i}=\operatorname{MLP}\left(p_{t}^{i},W_{\lambda}\right)$ (10)

The input is then fed into the Social-Multi-Head Attention Layer (SMHA). In
this layer, the features of the input are first split into $n_{heads}$
transformer heads and fed into the Multi-Head Attention (MHA) layer. The MHA
first computes the query, key, and value matrices at timestep $t$, as shown in
Eq. 11

$\left\\{Q_{t}=\operatorname{MLP}\left(\left\\{x_{t}^{i}\right\\}_{i=1}^{N},w_{q}\right),\quad
K_{t},\quad V_{t}\right\\}$ (11)

The inter-vehicular interaction is represented by an undirected graph
$G_{t}=\left\\{V_{t},E_{t}\right\\}$ where the nodes $V_{t}$ represent the
observed vehicles $V_{t}=\left\\{v\mid i=1,2,...,N\right\\}$ and the edges
represent the interaction as binary values between vehicles $i$ and $j$,
$E_{t}[i][j]=1$ if $p_{t}^{i}-p_{t}^{j}\leq d_{near}$, otherwise the value is
zero. The interaction graph is used to mask social attention.

The masked multi-head attention of head $h$ at timestep $t$ is then computed
as

$\operatorname{A}_{h}\left(Q_{t},K_{t},V_{t}\right)=\frac{\operatorname{Softmax}\left(\left[\left(q_{t}^{i}\right)^{T}k_{t}^{j}\right]_{\Gamma}\right)}{\sqrt{d_{k}}}\left[v_{t}^{i}\right]_{\Gamma},$
(12)
$\operatorname{MHA}\left(Q_{t},K_{t},V_{t}\right)=\operatorname{Concat}\left(A_{1},...,A_{n_{heads}}\right).$
(13)

where $\Gamma=\left\\{j\mid G_{t}[i,j]=1,j\in[1,N]\right\\}$ is the attention
mask. The output from the Multi-Head Attention layer of size $O_{MHA}$ is then
fed through an intra-layer MLP which helps boost training speed. This is then
followed by an Add & Norm layer. The Add layer is simply a residual connection
of the input $X_{\lambda}$ which adds the current output with the input. This
connection greatly helps models such as transformers reach convergence faster
by resolving the vanishing gradient problem so it is also used extensively in
our model [18]. The output of the Social Multi-Head Attention Layer is thus
computed as follows:

$\operatorname{Norm}(\operatorname{MLP}(O_{MHA})+X_{\lambda}).$ (14)

Positional Encoding, via a simple MLP is then added to the output which is
then fed through a Temporal Multi-Head Attention (TMHA) layer. The TMHA layer
has a similar architecture as SMHA along with a time mask that prevents the
current steps from accessing features from the relative future. For example,
if the current timestep is $t$ then the features from timestep $t+1$ to $T$
will be masked to zero. This prevents the model from overfitting and making
exclusive correlations on the training data. The output is also fed through
another multi-layer perceptron network (FFN). FFN is made up of two feed-
forward networks that successively map the feature vectors to a higher
dimension and then back to the original dimension. This is done to add more
model parameters so that the temporal attention vectors can be fed through
further layers such as the LSTM. Thus the output from the TMHA layer has both
social and temporal correlation encoded into it. We refer to the SMHA, TMHA,
and FFN layers compounded together as Socio-Temporal Multi-Head Attention
(STMHA) layer, as illustrated in Figure 1. We then stack $M$ number of STMHA
layers successively to extract more complex socio-temporal dependence from the
past trajectory information.

We propose a traditional LSTM encoder after the first STMHA layer stack to
extract the hidden memory information from the output tensor. This is done to
facilitate the hidden tensor passing to the LSTM decoder module. We propose
this sequence-to-sequence architecture as a replacement for the traditional
transformer architecture with an MHA-based transformer decoder. As mentioned
before, due to the autoregressive nature of the transformer decoder, it
carries the problem of accumulated errors. To decode the hidden-state tensor
and at the same time extract the rich socio-temporal dependence information
encoded into the hidden tensor from the decoded information, we propose an
STMHA-LSTM decoder. The hidden tensor at timestep $T+1$ is then fed into an
STMHA layer to extract and update the socio-temporal dependence at the next
decoding step. Such is repeated for every following decoding step. This is
done to further improve the accuracy of successive decoding steps, as well as
decode the future trajectories from the hidden information of the transformer
encoder with lower auto-regressive accumulated errors.

##  III Experiments

###  III-A Implementation Details

####  III-A1 Dataset

The proposed model was trained, validated, and tested on the BLVD dataset. It
consists of a total of 654 high-resolution videos which were extracted from
Changshu city, Jiangsu province. Each frame consists of the ID, 3D
coordinates, vehicle direction, and interactive behavior information of all
observed vehicles by the ego. These frame annotations were used to train, test
and evaluate our coordinate estimation algorithm and produce trajectories
which were then run through the encoder-decoder to make future predictions.
The prediction was then evaluated against the trajectory ground truth of the
BLVD dataset. In line with the existing literature [1, 12, 7] we have adopted
the root mean squared error (RMSE) between the prediction and the ground truth
as our evaluation metric. RMSE at prediction time
$t^{\prime},t\in[T+1,...,T+F]$ can be calculated as illustrated in Eq. 15.

$\operatorname{RMSE}_{t}=\sqrt{\frac{1}{L}\sum_{l=1}^{L}\left(\hat{Y}_{t}^{l}-Y_{t}^{l}\right)^{2}},$
(15)

where $\hat{Y}$ is the predicted positions and $Y$ is the Ground Truth
Position of the l-th testing sample at timestep $\textit{t}^{\prime}$ and L is
the total length of the test set.

####  III-A2 Hyperparameter Settings

For the observable region, we set it to be a 30m radial region. We set the
threshold $d_{near}$ to be 15m. We set the number of STMHA layers $M=2$ and
$n_{heads}$ of all the MHA modules to be 4. For the trajectory, we have used
the past timesteps $T$ to be 3 s and future timestamps F to be 5 s. For model
training, we have used the Adam [19] optimizer with $\eta$ = 0.001, $\beta$ =
0.999. The learning rate used is 0.0001 and the batch size of 32. The teacher-
forcing ratio used is 0.5.

###  III-B Ablative Tests

To defend the effectiveness of our architecture, we perform three ablative
experiments with the following three variants of our architecture

  1. 1.

TP:  Excludes the Trajectory Processor and uses trajectories extracted from
the dataset to test and validate.

  2. 2.

EST:  Excludes the STMHA in the Encoder.

  3. 3.

DST:  Excludes the STMHA in the Decoder.

Following Table I are the results for the above ablative setup.

TABLE I: Ablative Analysis Excluded | RMSE | Average
---|---|---
Module | 1s | 2s | 3s | 4s | 5s | Improvement
TP | 0.48 | 1.01 | 1.60 | 2.31 | 3.36 | -29.6%
EST | 0.67 | 1.48 | 2.10 | 3.03 | 4.70 | 34.8%
DST | 0.64 | 1.61 | 2.36 | 3.35 | 5.66 | 67.6%
Control | 0.57 | 1.32 | 1.86 | 2.66 | 3.82 |

At first, the results for the improvements in the EST and DST imply the
effectiveness of the STMHA layer. The Encoder STMHA layer proved to be crucial
in encoding socio-temporal interaction, as evident by the increase in error
through all timesteps resulting from its exclusion (EST). The Decoder STMHA
layer seems to play a vital role in reducing errors as well especially in
future timesteps as evident by the exponential rise in error due to its
exclusion in DST. Secondly, the exclusion of TP and directly importing the
annotated trajectories of the dataset seems to positively impact the model,
however, the error between the control (with TP) and the variant without TP
does not seem that significant, with an average improvement of only -29.6%.
This suggests that although TP cannot compare to manually constructed
annotations as in the BLVD dataset, it does not produce 3D coordinate
estimates that are too far off.

We also studied the mean distance error (MDE) (in meters) of the position
estimation and the Intersection over Union (IoU) of the bounding box estimated
vs the distance of the target from the ego. Figure 3 summarizes the findings.

Figure 3: MDE & IoU of Target vs Distance from Ego

###  III-C Comparison with Other Models

Table II summarizes our method’s performance with the existing literature. It
demonstrates our model’s strength, especially in enhancing the prediction
accuracy in future timesteps.

TABLE II: Comparative Analysis Model | RMSE
---|---
| 1s | 2s | 3s | 4s | 5s
CV | 0.73 | 1.78 | 3.13 | 4.78 | 6.68
V-LSTM | 0.68 | 1.65 | 2.91 | 4.46 | 6.27
S-LSTM | 0.65 | 1.31 | 2.16 | 3.25 | 4.55
CS-LSTM | 0.61 | 1.27 | 2.09 | 3.10 | 4.37
DSCAN | 0.58 | 1.26 | 2.03 | 2.98 | 4.13
SGAN | 0.57 | 1.32 | 2.22 | 3.26 | 4.40
HMNet | 0.50 | 1.13 | 1.89 | 2.85 | 4.04
Ours | 0.57 | 1.32 | 1.86 | 2.66 | 3.82
Figure 4: Visualization of Real-time Trajectory Prediction on CARLA

###  III-D Visualization

The trained model was also implemented on CARLA’s [20] Town10 in order to
qualitatively validate the strengths and weaknesses of the model through
visualized scenarios. Figure 4 demonstrates four distinct scenes each with
different levels of surrounding traffic, with the trajectories marked in
different colors. The performance of the model seems similar for the first
three scenes, with very low absolute path error. Our model was designed to
deal with congested traffic by extracting socio-temporal dependencies and
making a prediction, and this is further evident by these three scenarios.
However, as can be seen in Figure 4(d), the absolute path increasingly
deviates from GT in lane change scenarios, even though the model predicts an
accurate direction of the lane change. We believe this is due to the highway
dataset of BLVD not having enough lane change cases, which causes a sort of
underfitting of the model’s prediction in those scenarios.

Visualization tests were also carried out to check the performance of the 3D
bounding box prediction, as illustrated in Figure 5. This demonstrates our
technique of converting from 2D to 3D bounding box via the Trajectory
Processor module. A wrong 2D prediction resulting in a faulty bounding box can
be observed in this scenario (on the building to the left). Our method relies
strongly on the 2D box prediction in order to estimate its 3D counterpart,
which can be regarded as one of its limitations. So, if the 2D box is
inaccurate, this would also propagate the errors to the 3D box.

Figure 5: Visualization of Bounding Box Reconstruction on CARLA

##  IV Conclusion and Future Work

This paper has introduced an end-to-end trajectory prediction framework,
enabling trajectories to be extracted and forecasted at the same time from raw
video inputs. The framework is able to extract 3D position and concurrently
make accurate future predictions by modeling the complex socio-temporal
interaction. The strengths of the model have been proven to be in congested
highway driving scenarios and longer predictions, while the limitations lie in
the position estimation and lane change scenarios. Further work can focus on
implementing more efficient 3D position estimation techniques and possibly
also implementing driving style into the prediction in order to improve the
performance on lane change scenarios. Further suggestions may include focusing
on more diverse driving scenarios such as urban driving with other road
participants such as pedestrians and cyclists.

## References

  * [1] K. Messaoud, I. Yahiaoui, A. Verroust-Blondet, and F. Nashashibi, “Attention based Vehicle Trajectory Prediction,” _IEEE Trans. Intell. Vehicles_ , vol. 6, pp. 175–185, 2021.
  * [2] S. Hochreiter and J. Schmidhuber, “Long Short-Term Memory,” _Neural Comput._ , vol. 9, pp. 1735–1780, 1997.
  * [3] J. Chung, C. Gulcehre, K. Cho, and Y. Bengio, “Empirical Evaluation of Gated Recurrent Neural Networks on Sequence Modeling,” _arXiv_ , arXiv:1412.3555, 2016.
  * [4] M. Althoff and A. Mergel, “Comparison of Markov Chain Abstraction and Monte Carlo Simulation for the Safety Assessment of Autonomous Cars,” _IEEE Trans. Intell. Transp. Syst._ , vol. 12, pp. 1237–1247, 2011.
  * [5] Y. Xing, C. Lv, X. Mo, Z. Hu, C. Huang, P. Hang, and D. Cao, “Toward Safe and Smart Mobility: Energy-Aware Deep Learning for Driving Behavior Analysis and Prediction of Connected Vehicles,” _IEEE Trans. Intell. Transp. Syst._ , vol. 22, pp. 4267–4280, 2021.
  * [6] Y. Xing, C. Huang, C. Lv, L. Yahui, W. Hong, and D. Cao, “A Personalized Deep Learning Approach for Trajectory Prediction of Connected Vehicles,” _SAE International Journal of Advances and Current Practices in Mobility_ , 2020-01-0759, 2020.
  * [7] C. Yu, X. Ma, J. Ren, H. Zhao, and S. Yi, “Spatio-Temporal Graph Transformer Networks for Pedestrian Trajectory Prediction,” _Computer Vision–ECCV_ , vol. 12357, pp. 507–523, 2020.
  * [8] K. Chen, G. Chen, D. Xu, L. Zhang, Y. Huang, and A. Knoll, “NAST: Non-autoregressive Spatial-Temporal Transformer for Time Series Forecasting,” _arXiv_ , arXiv:2102.05624, 2021.
  * [9] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, A. N. Gomez, L. Kaiser, and I. Polosukhin, “Attention Is All You Need,” _IEEE Adv. Neural Inf. Process. Syst._ , pp. 5998–6008, 2017.
  * [10] Next Generation Simulation (NGSIM). Available online: http://ops.fhwa.dot.gov/trafficanalysistools/ngsim.html (accessed on 20 September 2022).
  * [11] J. Xue, J. Fang, T. Li, B. Zhang, P. Zhang, Z. Ye, and J. Dou, “Blvd: Building a Large-Scale 5D Semantics Benchmark for Autonomous Driving,” _arXiv_ , vol. 1903.06405, 2019.
  * [12] N. Deo and M. M. Trivedi, “Convolutional Social Pooling for Vehicle Trajectory Prediction,” _IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops (CVPRW)_ , pp. 1468–1476, 2018.
  * [13] N. Wojke, A. Bewley, and D. Paulus, “Simple online and realtime tracking with a deep association metric,” _Proceedings of the 2017 IEEE International Conference on Image Processing (ICIP)_ , pp. 3645–3649, 2017.
  * [14] W. Luo, J. Xing, A. Milan, X. Zhang,W. Liu, T. Kim “Multiple Object Tracking: A Literature Review,” _arXiv_ , arXiv:1409.7618v5, 2022.
  * [15] A. Mousavian, D. Anguelov, J. Flynn, and J. Kosecka, “3D Bounding Box Estimation Using Deep Learning and Geometry,” _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , pp. 7074–7082, 2017.
  * [16] F. Hasan, H. Huang, “MALS-Net: A Multi-Head Attention-Based LSTM Sequence-to-Sequence Network for Socio-Temporal Interaction Modelling and Trajectory Prediction,” _Sensors_. vol. 23(1), pp. 530, 2023.
  * [17] S. Santurkar, D. Tsipras, A. Ilyas, and A. Madry, “How Does Batch Normalization Help Optimization?,” _arXiv_ , vol. 1805.11604, 2018.
  * [18] R. Pascanu, T. Mikolov, and Y. Bengio, “On the difficulty of training recurrent neural networks,” _arXiv_ , vol. 1211.5063, 2012.
  * [19] D. P. Kingma and J. Ba, “Adam: A Method for Stochastic Optimization,” _arXiv_ , vol. 1412.6980, 2014.
  * [20] A. Dosovitskiy, G. Ros, F. Codevilla, A. Lopez, and V. Koltun, “CARLA: An Open Urban Driving Si


# Reconstruction of Stochastic Dynamics from Large Streamed Datasets

William Davis<EMAIL_ADDRESS>Cecil H. and Ida M. Green Institute of
Geophysics and Planetary Physics, Scripps Institution of Oceanography,
University of California, San Diego, La Jolla, California 92037, USA

###### Abstract

The complex dynamics of physical systems can often be modeled with stochastic
differential equations. However, computational constraints inhibit the
estimation of dynamics from large time-series datasets. I present a method for
estimating drift and diffusion functions from inordinately large datasets
through the use of incremental, online, updating statistics. I demonstrate the
validity and utility of this method by analyzing three large, varied synthetic
datasets, as well as an empirical turbulence dataset. This method will
hopefully facilitate the analysis of complex systems from exceedingly large,
“big data” scientific datasets, as well as real-time streamed data.

††preprint: APS/123-QED

##  I Introduction

The dynamics of complex systems with many degrees of freedom can often be
modeled as continuous-time stochastic processes [1]. When a system is modeled
by a stochastically forced, scalar, first-order differential equation, the
temporal evolution of a quantity $X(t)$ is described by a Langevin-type
equation [2]:

$\frac{d}{dt}X(t)=f(X)+g(X)\Gamma(t),$ (1)

where a separation of scales partitions the dynamics of $X(t)$ into slow
changes modulated by $f(X)$, and rapidly-varying changes modulated by $g(X)$.
Fluctuations are driven by Gaussian white noise $\Gamma(t)$, with
$\langle\Gamma(t)\rangle=0$ and
$\langle\Gamma(t)\Gamma(t^{\prime})\rangle=\delta(t-t^{\prime})$. Here and
throughout, the Itô interpretation is adopted. If $X(t)$ contains no
discontinuous jumps, then the evolution of the probability density function
can be described by the Fokker-Planck equation [2]

$\frac{\partial}{\partial
t}p(x,t|x^{\prime},t^{\prime})=\Bigg{[}-\frac{\partial}{\partial
x}D^{(1)}(x)\\\ +\frac{\partial^{2}}{\partial
x^{2}}D^{(2)}(x)\Bigg{]}p(x,t|x^{\prime},t^{\prime})$ (2)

where $p(\circ|\circ)$ is the transition probability, and $x$ and $x^{\prime}$
are state variables of $X$. The Fokker-Planck equation contains the Kramers-
Moyal (KM) coefficients

$D^{(k)}(x)=\lim_{\tau\rightarrow
0}\frac{1}{k!\tau}\int_{-\infty}^{\infty}\big{[}x^{\prime}-x\big{]}^{k}p(x^{\prime},t+\tau|x,t)\
dx^{\prime}.$ (3)

The $k=1$ and $k=2$ KM coefficients are called the drift and diffusion
functions, respectively, and they correspond to terms in the dynamical
equation (1), with $f(x)=D^{(1)}(x)$ and $g(x)=\sqrt{2D^{(2)}(x)}$.

It has been shown that KM coefficients—and hence drift and diffusion
functions—can be estimated from empirical samples of $X(t)$, using a
conditional averaging technique called “direct estimation” [3, 4]. Direct
estimation and descendant methods [e.g., 5, 6, 7, 8, 9, 10, 11] have been
applied to time-series data in various fields of science [12, 13], including
turbulence [14, 15, 16], wind energy [17], climate data [18, 6], and
geomagnetic field variations [19, 20].

Although the calculation of KM coefficients is conceptually simple [3, 21],
estimations are prone to bias, especially in areas of rarely sampled state
space [22]. Inaccuracies are particularly apparent for processes with heavy
tails, or for systems that exhibit rare, transient dynamics. Attempts to
resolve KM coefficients in rarely sampled regions by reducing the resolution
of conditioning also result in biased drift and diffusion estimates [21].

A rudimentary but effective solution to the sampling problem is to perform
analyses on datasets that are as large as possible. This approach is effective
because the estimation bias of KM coefficients scales as $1/\sqrt{N\Delta t}$,
where $N$ is the number of samples and $\Delta t$ is the sampling interval
[23]. Indeed, in the era of “big data,” there is growing interest in
estimating drift and diffusion functions for increasingly large scientific
datasets [24, 13]. However, existing KM procedures calculate KM coefficients
using offline methods [25, 22, 26], requiring the complete dataset to be
available at once, and with memory requirements that scale with the number of
data points. Large datasets are often incompatible with offline methods,
either because the data cannot fit into computer memory, or because the data
originates from arbitrarily large data streams [27]. An alternative approach
is to use online methods, which incrementally update statistical estimates
from streamed data, arriving one data point at a time [28, 29]. In this paper
I present an online method of computing KM coefficients from streamed time-
series data, enabling the estimation of drift and diffusion functions from
large time-series datasets which are unreachable with previous methods.

##  II Estimation of Conditional Moments

Consider a finite sample of $N$ points in $X(t)$ from process (1), denoted as

$\mathscr{S}_{N}:=\\{(t_{1},X_{1}),(t_{2},X_{2}),\dots,(t_{N},X_{N})\\}.$ (4)

Here I assume a regular sampling interval $\Delta t$. The aim is to use these
data to construct non-parametric estimates of drift and diffusion coefficients
of the Langevin-type equation that generated $X(t)$. Estimation of drift and
diffusion coefficients is conducted at a set of $N_{x}$ evaluation points in
$x$, represented by the vector

$\bm{\mathcal{X}}:=\left[x_{1},x_{2},\dots,x_{N_{x}}\right].$ (5)

Drift and diffusion estimates at these points will be denoted by the vector
$\mathbf{\hat{D}}^{(k)}$, with
$\hat{D}^{(k)}_{j}:=\hat{D}^{(k)}(\mathcal{X}_{j})$. Estimation of
$\mathbf{\hat{D}}^{(k)}$ requires evaluation of the conditional process
increments—or “conditional moments” [9]—in (3), namely

$M^{(k)}(\tau,x)=\int_{-\infty}^{\infty}[x^{\prime}-x\big{]}^{k}p(x^{\prime},t+\tau|x,t)\
dx^{\prime},$ (6)

for $k=1,2$. As the $\tau\rightarrow 0$ limit in (3) cannot be performed for
empirical data, (6) is estimated at a set of $N_{\tau}$ evaluation points in
$\tau$ values, represented by the vector

$\bm{\mathcal{T}}:=\left[\Delta t,2\Delta t,\dots,N_{\tau}\Delta
t\right]^{T}.$ (7)

Estimates of conditional moments (6) are performed at all points in
$\bm{\mathcal{T}}$ and $\bm{\mathcal{X}}$, and will be denoted as
${N_{\tau}\times N_{x}}$ matrices $\mathbf{\hat{M}}^{(k)}$, with
$\hat{M}^{(k)}_{ij}:=\hat{M}^{(k)}(\mathcal{T}_{i},\mathcal{X}_{j})$. I now
outline an existing estimation procedure for conditional moments, before
proposing online updating formulas.

###  II.1 Offline calculation

One method of estimating conditional moments is kernel-based regression (KBR)
[21]. A chosen kernel function $K(\cdot)$ applies conditioning on the state
variable, $x$, and, assuming ergodicity, the estimators for (6) can be written
as

$\hat{M}^{(k)}_{ij}=\frac{\sum\limits_{n=1}^{N-i}K_{h}(\mathcal{X}_{j}-X_{n})\big{[}X_{n+i}-X_{n}\big{]}^{k}}{\sum\limits_{n=1}^{N-i}K_{h}(\mathcal{X}_{j}-X_{n})},$
(8)

for $k=1,2$, where $K_{h}(\cdot)=K(\cdot/h)/h$ is a scaling of the kernel with
bandwidth $h$. The methods presented here are independent of the chosen
kernel. Here I use the Epanechnikov kernel [30]

$K(x)=\begin{cases}\frac{3}{4}(1-x^{2})&\text{if~{}$x^{2}<1$},\\\
0&\text{otherwise},\end{cases}$ (9)

which has computationally favorable properties [31]. If kernel conditioning is
replaced with bin counting, the estimation becomes histogram-based regression
(HBR) [3]. Some studies also analyze the variance of the conditional process
increments (e.g., [32, 33, 11]). I will refer to this quantity as the
“conditional variance,” and denote it as

$\hat{M}^{(2^{*})}_{ij}=\frac{\sum\limits_{n=1}^{N-i}K_{h}(\mathcal{X}_{j}-X_{n})\left(\big{[}X_{n+i}-X_{n}\big{]}-\hat{M}^{(1)}_{ij}\right)^{2}}{\sum\limits_{n=1}^{N-i}K_{h}(\mathcal{X}_{j}-X_{n})}.$
(10)

Both HBR and KBR are implemented in modern software libraries [25, 22, 34,
26], and can be extended to irregularly sampled time-series data [35].
However, these offline methods require the entire input $\mathscr{S}_{N}$ to
be available at once: the entire calculation must be repeated if more data is
appended to $\mathscr{S}_{N}$.

###  II.2 Online calculation

I now present formulas for updating sample conditional moments (8), previously
calculated from $\mathscr{S}_{N-1}$, with a single new observation
$(t_{N},X_{N})$. I refer to this approach as “online kernel-based regression”
(OKBR). To facilitate indexing, the subscript notation $[\dots]\big{|}_{n}$
denotes a quantity calculated from the first $n$ observations. The updating
formulas are written as (see Appendix)

$\hat{M}^{(k)}_{ij}\big{|}_{N}=\hat{M}^{(k)}_{ij}\big{|}_{N-1}+K_{h}(\mathcal{X}_{j}-X_{N-i})\\\
\times\left(\left[X_{N}-X_{N-i}\right]^{k}-\hat{M}^{(k)}_{ij}\big{|}_{N-1}\right)\Big{/}W_{ij}\big{|}_{N},$
(11)

for $k=1,2$, where $W_{ij}\big{|}_{N}$ are cumulative weights

$W_{ij}\big{|}_{N}=W_{ij}\big{|}_{N-1}+K_{h}\left(\mathcal{X}_{j}-X_{N-i}\right).$
(12)

To define a corresponding updating formula for (10), I introduce the
intermediate quantity $S_{ij}\big{|}_{N}$, which corresponds to the weighted
sum of squares of differences from the current mean,

$S_{ij}\big{|}_{N}:=\sum\limits_{n=1}^{N-i}K_{h}(\mathcal{X}_{j}-X_{n})\Big{(}\big{[}X_{n+i}-X_{n}\big{]}-\hat{M}^{(1)}_{ij}\big{|}_{N}\Big{)}^{2},$
(13)

and is related to (10) by

$\hat{M}^{(2^{*})}_{ij}\big{|}_{N}=S_{ij}\big{|}_{N}\bigg{/}W_{ij}\big{|}_{N}.$
(14)

The corresponding online formula is (see Appendix .2)

$S_{ij}\big{|}_{N}=S_{ij}\big{|}_{N-1}+K_{h}(\mathcal{X}_{j}-X_{N-i})\\\
\times\left(\left[X_{N}-X_{N-i}\right]-\hat{M}^{(1)}_{ij}\big{|}_{N-1}\right)\\\
\times\left(\left[X_{N}-X_{N-i}\right]-\hat{M}^{(1)}_{ij}\big{|}_{N}\right),$
(15)

These formulas have been constructed to avoid numerical instability and loss
of precision [36, 37]. In the next section, I validate the presented methods
on three synthetic datasets.

##  III Numerical Examples

###  III.1 Ornstein-Uhlenbeck process

I examine a simple example where the drift and diffusion functions are set as

$\displaystyle D^{(1)}(x)$ $\displaystyle=-x,$ (16a) $\displaystyle
D^{(2)}(x)$ $\displaystyle=1.$ (16b)

I numerically integrate [38] this process using a sampling interval of $\Delta
t=10^{-3}$ for $N=10^{7}$ data-points. I estimate conditional moments at 26
equally-spaced points in the range $[-5,5]$ using a bandwidth of $h=0.4$, and
perform time sampling at a single time-step $\bm{\mathcal{T}}=[\Delta t]$. I
conduct estimation using both the KBR formulas (8) and OKBR formulas (11). To
illustrate the ability of OKBR to conduct analysis on an inordinately large
dataset, I also repeat the OKBR estimation for a simulated time-series with
$N=10^{10}$ data-points. For all three cases, I estimate drift and diffusion
coefficients from the conditional moments using direct estimation [3]

$\mathbf{\hat{D}}^{(k)}=\frac{1}{k!\Delta t}\mathbf{\hat{M}}^{(k)}.$ (17)

Results are shown in Fig. 1.

Figure 1: Drift (upper) and diffusion (lower) estimates for example III.1. The
true drift and diffusion functions are indicated with dashed black lines.
Estimates from $N=10^{7}$ data-points conducted with KBR and OKBR are shown
with red “$\times$” crosses and green “$+$” crosses, respectively. Estimates
from $N=10^{9}$ data points conducted with OKBR are shown with open blue
circles. For the $N=10^{7}$ simulation, estimates of drift and diffusion
coefficients at $|x|\gtrsim 2$ are noisy or non-existent, due to sparse
sampling.

I find that for the $N=10^{7}$ case, KBR and OKBR give identical estimates for
the drift and diffusion coefficients, and the coefficients $|x|\lesssim 2$ are
estimated fairly. However, at the rarely sampled edges, either large errors
are present or there are no samples available to make an estimate. For the
$N=10^{10}$ case, OKBR accurately recovers the drift and diffusion
coefficients over the entire estimation range. It is not possible to use KBR
on the $N=10^{10}$ dataset, as the data does not fit within computer memory.

####  III.1.1 Empirical performance of estimation procedures

To empirically benchmark the time and space requirements for KBR and OKBR, I
repeat the estimations in Sec. III.1, varying the number of data-points $N$
and leaving all other parameters unchanged. The number of data-points
considered is $N\in(10^{4},10^{5},\dots,10^{10})$, however estimation using
the largest dataset using KBR is not possible due to memory requirements.
Table 1 shows the benchmark results. Both methods show linear scaling in time.
KBR shows linear scaling in space, whereas the space requirements of OKBR
scale to a constant, independent of $N$.

Table 1: Time and space requirements for KBR and OKBR, varying the number of data points $N$ in the time-series dataset. Estimation using KBR is not possible for the $N=10^{10}$ dataset due to memory restrictions. Similar scalings were found for examples III.2 and III.3. For KBR, the memory estimate is dominated by the size of the time series itself.  | KBR | OKBR
---|---|---
N | Time (s) | Space (GB) | Time (s) | Space (GB)
${10}^{4}$ |  $4.49\times 10^{-2}$ |  $1.28\times 10^{-4}$ |  $4.08\times 10^{-2}$ |  $7.44\times 10^{-5}$
${10}^{5}$ |  $3.89\times 10^{-1}$ |  $8.48\times 10^{-4}$ |  $3.65\times 10^{-1}$ |  $1.03\times 10^{-4}$
${10}^{6}$ |  $4.05\times 10^{0}$ |  $8.05\times 10^{-3}$ |  $3.65\times 10^{0}$ |  $1.63\times 10^{-4}$
${10}^{7}$ |  $4.05\times 10^{1}$ |  $8.00\times 10^{-2}$ |  $3.64\times 10^{1}$ |  $1.63\times 10^{-4}$
${10}^{8}$ |  $4.08\times 10^{2}$ |  $8.00\times 10^{-1}$ |  $3.64\times 10^{2}$ |  $1.64\times 10^{-4}$
${10}^{9}$ |  $4.06\times 10^{3}$ |  $8.00\times 10^{0}$ |  $3.64\times 10^{3}$ |  $1.64\times 10^{-4}$
${10}^{10}$ |  |  |  $3.64\times 10^{4}$ |  $1.64\times 10^{-4}$

###  III.2 Tri stable system

I consider a system which exhibits poorly sampled regions of state space,
arising from fast, transient dynamics through unstable (or metastable) states.
One natural example of such a system is the time-variability of the axial
dipole moment of Earth’s geomagnetic field, which shows two prominently stable
states at positive and negative polarity, and an unstable (or possibly
metastable) “weak” state during polarity transitions [39, 40, 41]. The
qualitative dynamics of this system can be represented by the toy model

$\displaystyle D^{(1)}(x)$ $\displaystyle=-x+27x^{3}-26x^{5},$ (18a)
$\displaystyle D^{(2)}(x)$ $\displaystyle=\frac{7}{10}.$ (18b)

This system is characterized by two strong attractors at $x=\pm 1$, and one
weaker, rarely sampled attractor at $x=0$. Here, one might aim to determine
the stability of the middle state from empirical data.

I integrate system (18) with a sampling interval of $\Delta t=10^{-4}$, using
$N=5\times 10^{7}$ and $N=10^{10}$ points for KBR (8) and OKBR (11),
respectively. I estimate conditional moments at 45 equally spaced points in
the interval $[-1.4,1.4]$ using a bandwidth of $h=0.03$, and perform sampling
in $\tau$ at a series of time steps $\bm{\mathcal{T}}=\left[\Delta t,2\Delta
t,3\Delta t,4\Delta t\right]^{T}$. I estimate drift and diffusion coefficients
in the $\tau\rightarrow 0$ limit in (3) by minimizing

$V\left(\mathbf{\hat{D}}^{(k)}\right)=\left|\left|\mathbf{\hat{M}}^{(k)}-\bm{\mathcal{T}}\mathbf{\hat{D}}^{(k)}\right|\right|^{2},$
(19)

using ordinary least squares. Results are shown in Fig. 2.

Figure 2: Drift (upper) and diffusion (lower) estimates for example III.2. The
inset in the upper plot shows a zoomed in section close to the origin. The
inset in the lower plot shows a section of the time-series data. The true
drift and diffusion functions are indicated with dashed black lines. Estimates
using KBR with $N=5\times 10^{7}$ data points are shown with red “$\times$”
crosses. Estimates using OKBR with $N=10^{10}$ data-points are shown with
green “$+$” crosses. A repeated estimation using OKBR and $N=5\times 10^{7}$
data points gives identical results to KBR, and is not plotted for
conciseness.

I find that for the $N=5\times 10^{7}$ case, KBR is able to reasonably recover
the drift and diffusion coefficients close to the attractors at $x=\pm 1$.
However, poor estimates are made for the rarely sampled transitions, for
$x\in[-0.5,0.5]$, and the details of stability at $x=0$ are unresolvable. For
the $N=10^{10}$ case, OKBR accurately recovers the drift and diffusion
coefficients across the entire sampling domain, revealing the presence of the
weak attractor at $x=0$.

###  III.3 Multiplicative and correlated noise

I consider a system with a multiplicative diffusion term and an exponentially
correlated noise source $\eta(t)$,

$\displaystyle\frac{d}{dt}X$
$\displaystyle=D^{(1)}(X)+\sqrt{2D^{(2)}(X)}\eta(t),$ (20a)
$\displaystyle\frac{d}{dt}\eta$
$\displaystyle=-\frac{1}{\theta}\eta+\frac{1}{\theta}\Gamma(t),$ (20b)

where

$\displaystyle D^{(1)}(x)$
$\displaystyle=-\frac{1}{8}-\frac{9}{4}x-\frac{4}{15}x^{3},$ (21a)
$\displaystyle D^{(2)}(x)$
$\displaystyle=1+\frac{1}{50}x^{2}+\frac{1}{40}x^{4},$ (21b)

and $\theta=0.01$ is the correlation time of the noise $\eta(t)$, and
$\Gamma(t)$ is internal Gaussian white noise. Only the time-series of $X(t)$
is observed.

I analyze process (20) and (21) using the non parametric inversion method of
Lehle and Peinke [11], assuming that the timescale $\theta$ has already been
estimated (e.g., using [29]). This method requires estimation of the sample
conditional mean—$k=1$ in (8) and (11)—as well as the conditional variance,
(10) and (14). I integrate the process with a sampling interval of $\Delta
t=5\times 10^{-3}$, using $N=10^{7}$ and $N=5\times 10^{9}$ points for KBR and
OKBR, respectively. I estimate the conditional quantities
$\mathbf{\hat{M}}^{(k)}$ at 100 equally spaced points in the interval
$[-2.5,2.5]$ using a bandwidth of $h=0.01$, and perform sampling in $\tau$
using 25 time steps, $\bm{\mathcal{T}}=\left[\Delta t,\dots,25\Delta
t\right]^{T}$. To estimate the drift and diffusion coefficients using the
method of Lehle and Peinke [11], I decompose the sample conditional mean and
variance into basis functions $r_{i}(\tau,\theta)$ and coefficients
$\lambda_{i}^{(k)}(x)$, given by

$M^{(k)}(x,\tau)\approx\sum_{i=1}^{3}\lambda_{i}^{(k)}(x)r_{i}(\tau,\theta).$
(22)

Here, the basis functions are

$\displaystyle r_{1}(\tau;\theta)=$
$\displaystyle\tau-\theta(1-e^{-\tau/\theta}),$ (23) $\displaystyle
r_{2}(\tau;\theta)=$ $\displaystyle\tau^{2}/2-\theta r_{1}(\tau;\theta),$
$\displaystyle r_{3}(\tau;\theta)=$ $\displaystyle\tau^{3}/6-\theta
r_{2}(\tau;\theta),$

and are expressed in matrix form with elements

$R_{ij}:=r_{j}(\mathcal{T}_{i}).$ (24)

I solve for the coefficients by minimizing

$V\left(\bm{\lambda}^{(k)}\right)=\left|\left|\mathbf{\hat{M}}^{(k)}-\bm{R}\bm{\lambda}^{(k)}\right|\right|^{2},$
(25)

using ordinary least squares. Finally, I use the $i=1$ components of the
coefficients to solve differential algebraic equations for estimates of the
drift and diffusion coefficients $\mathbf{\hat{D}}^{(k)}$; see Lehle and
Peinke [11] for details. The estimated drift and diffusion coefficients are
shown in Fig. 3.

Figure 3: Drift (upper) and diffusion (lower) estimates, for example III.3.
The true drift and diffusion functions are indicated with dashed black lines.
Estimates using KBR with $N=10^{7}$ data points are shown with red “$\times$”
crosses. Estimates using OKBR with $N=5\times 10^{9}$ data points are shown
with green “$+$” crosses. A repeated estimation from $N=10^{7}$ data points
using OKBR gives identical results to KBR, and is not plotted for conciseness.

I find that for the $N=10^{7}$ case, KBR is able to recover the drift and
diffusion coefficients in the range $x\in[-1,1]$, but poor estimates are made
in the rarely sampled tails. For the $N=5\times 10^{9}$ case, OKBR is able to
accurately recover the drift and diffusion coefficients over a much larger
range. To illustrate the consequences of poorly resolved tails, I use
$\mathbf{\hat{D}}^{(k)}$ to estimate the parametric coefficients of the
diffusion function,

$D^{(2)}(x)=A+Bx^{2}+Cx^{4}.$ (26)

Parameter estimates in Table 2 show that both the quadratic and quartic
coefficients are poorly resolved for the KBR case, with uncertainty intervals
overlapping zero. However, the increased resolution that OKBR enables results
in accurate parameter estimation.

Table 2: Fit results, with $2\sigma$ uncertainties and $R^{2}$ values. N | A | B | C | $R^{2}$
---|---|---|---|---
True | $1.000$ | $0.020$ | $0.025$ |
${10}^{7}$ | $0.946\pm 0.133$ | $0.113\pm 0.133$ | $0.002\pm 0.024$ |  $0.301$
$5\times{10}^{9}$ | $1.001\pm 0.005$ | $0.023\pm 0.005$ | $0.024\pm 0.001$ |  $0.998$

##  IV Application to Turbulence Data

To illustrate one possible application of OKBR, I examine a turbulence dataset
from Fuchs  _et al._ [34]. This dataset—originally published by Renner  _et
al._ [15]—comes from a turbulent air jet experiment, where time-variable
observations of local air velocity were made using hot-wire measurements. The
dataset comprises $N=1.25\times 10^{7}$ points sampled at 8 kHz, although
other turbulence datasets can be orders of magnitude larger [42].

The data can be used to investigate a statistical description of a turbulent
cascade [43]. The measurements, under the assumption of Taylor’s hypothesis of
frozen turbulence, reflect spatial velocity variations $u(x)$. Increments of
these velocity variations

$\xi_{n,i}:=\xi(x_{n},r_{i})=u(x_{n})-u(x_{n}-r_{i}),$ (27)

define a “zooming-in” process in $\xi$ for decreasing $r$. Following the
phenomenological model of Friedrich and Peinke [14], velocity increments
evolve as a Markov process in scale $r$. From this, the turbulent cascade is
interpreted as a stochastic process described by a Fokker-Planck equation
evolving through a sequence of velocity increments
$\xi_{n,0},\xi_{n,1},\xi_{n,2},\dots,$ at increasingly smaller
$r_{0}>r_{1}>r_{2}>\dots$ scales. One can use the empirical velocity
measurements to not only verify the Markov property of $\xi(r)$, but also to
estimate the corresponding drift and diffusion coefficients [43].

The conditional moments for two increment scales separated by $\delta$ are
defined as

$M^{(k)}(\delta,\xi,r,u_{N})=\\\
\int_{-\infty}^{\infty}[\xi^{\prime}(r-\delta,u_{N})-\xi(r,u_{N})\big{]}^{k}p(\xi^{\prime}|\xi,u_{N})\
d\xi^{\prime},$ (28)

for $k=1,2$. Then, the KM coefficients are given by [43]

$D^{(k)}(\xi,r,u_{N})=\frac{r}{k!}\lim_{\delta\rightarrow
0}\frac{1}{\delta}M^{(k)}(\delta,\xi,r,u_{N})$ (29)

Analogously to (11), the online formulas for the estimator of (28) can be
written as

$\hat{M}^{(k)}_{ij}\big{|}_{N}=\hat{M}^{(k)}_{ij}\big{|}_{N-1}+K_{h}(\mathcal{X}_{j}-\xi_{N,0})\\\
\times\left(\left[\xi_{N,i}-\xi_{N,0}\right]^{k}-\hat{M}^{(k)}_{ij}\big{|}_{N-1}\right)\Big{/}W_{ij}\big{|}_{N},$
(30)

where

$W_{ij}\big{|}_{N}=W_{ij}\big{|}_{N-1}+K_{h}\left(\mathcal{X}_{j}-\xi_{N,0}\right).$
(31)

I analyze the turbulence dataset comparably to Fuchs  _et al._ [34] by
normalizing the velocity by its variance, $\sigma$, and estimating conditional
moments using the same parameters described in Fuchs  _et al._ [34], their
Fig. 23. I use OKBR with a boxcar kernel and a bandwidth of $h=0.038$ to
estimate conditional moments at a range of scales separated by $\delta$, from
$\Delta_{EM}<\delta<2\Delta_{EM}$, where $\Delta_{EM}$ is the Einstein-Markov
length. KM coefficients are estimated in the $\delta\rightarrow 0$ limit
through linear extrapolation. The estimated drift and diffusion coefficients
are shown in Fig. 4, reproducing the previously determined results of Fuchs
_et al._ [34].

Figure 4: Drift (upper) and diffusion (lower) estimates, for example IV. The
coefficients $D^{(k)}(\xi,r)$ are estimated with respect to velocity increment
$\xi$ for a fixed scale $r=2.7\lambda=3.2\Delta_{EM}$, where $\lambda$ is the
Taylor length scale. Estimates from Fuchs  _et al._ [34] are shown with red
“$\times$” crosses. Estimates using OKBR are shown with green “$+$” crosses.

##  V Discussion and Conclusion

I present online updating formulas for estimating conditional moments and
variance from time-series data. These formulas enable the non-parametric
estimation of drift and diffusion functions from arbitrarily large datasets,
without requiring the entire set of input data to be available at once. I
demonstrate this with three numerical examples. Even for datasets that far
exceed the working memory of most computers, OKBR is able to generate accurate
estimates of drift and diffusion functions. OKBR is also applied to a
turbulence dataset, and the estimated drift and diffusion functions reproduce
previously determined results. Although OKBR’s sequential calculations do
inhibit vectorized optimizations [22], the constant memory usage enables the
analysis of exceedingly large scientific datasets. This method could thus be
incorporated into existing software packages (e.g., [25, 22, 34, 26]).
Furthermore, if the updating step is faster than the data sampling interval—as
is the case in example IV—OKBR may be used for real-time analysis of high-
frequency streamed data.

The method presented here is demonstrated in one dimension; however,
extensions to higher dimensions are straightforward. Extensions cannot be
assumed for higher-order conditional moments ($k>2$ in $\hat{M}^{(k)}_{ij}$),
as updating formulas for skewness, kurtosis, and other moments are nontrivial
[44]. Further work should seek to extend the online framework to higher-order
conditional moments, which would aid the analysis of jump-diffusion processes
[45]. Similar online approaches might also be used to estimate Einstein-Markov
length-scales from real-time turbulence experiments [15, 46, 47].

Although OKBR reduces the memory complexity to calculate conditional moments
from $\mathcal{O}(N)$ to $\mathcal{O}(1)$, the time complexity remains at
$\mathcal{O}(N)$: competitive with traditional offline methods, as well as
recent polynomial-time approaches [48]. However, as detailed by Chan  _et al._
[49], online formulas can sometimes be altered for calculation by multiple
processing units in parallel. It may thus be possible to estimate conditional
moments in sub linear time [50].

## Data Availability

A prototype Julia implementation of the estimation procedure is available at
[51]. The dataset from Fuchs  _et al._ [34] in Sec. IV is used under the GNU
General Public License (GPL) version 3.

###### Acknowledgements.

I thank Matthias Morzfeld, Catherine Constable, Katherine Armstrong, and three
anonymous reviewers for helpful discussions and comments which benefited this
research. This work is supported by the Cecil H. and Ida M. Green Foundation’s
John W. Miles postdoctoral fellowship in theoretical and computational
geophysics.

## Appendix: Derivation of Incremental Quantities

###  .1 Weights and conditional moments

First I define the cumulative weights,

$W_{ij}\big{|}_{N}:=\sum\limits_{n=1}^{N-i}K_{h}\left(\mathcal{X}_{j}-X_{n}\right).$
(A1)

This is rearranged to permit incremental updates

$W_{ij}\big{|}_{N}=W_{ij}\big{|}_{N-1}+K_{h}\left(\mathcal{X}_{j}-X_{N-i}\right).$
(A2)

Next I derive incremental formulas for conditional moments (8). Identifying
the denominator of (8) as (A1) and rearranging gives

$\hat{M}^{(k)}_{ij}\big{|}_{N}\cdot
W_{ij}\big{|}_{N}=\sum\limits_{n=1}^{N-i}K_{h}(\mathcal{X}_{j}-X_{n})\big{[}X_{n+i}-X_{n}\big{]}^{k}.$
(A3)

Separating the last term in the sum and substituting (A2) gives

$\hat{M}^{(k)}_{ij}\big{|}_{N}\cdot W_{ij}\big{|}_{N}=\\\
\hat{M}^{(k)}_{ij}\big{|}_{N-1}\cdot\Big{(}W_{ij}\big{|}_{N}-K_{h}(\mathcal{X}_{j}-X_{N-i})\Big{)}\\\
+K_{h}(\mathcal{X}_{j}-X_{N-i})\big{[}X_{N}-X_{N-i}\big{]}^{k}.$ (A4)

Finally, dividing by $W_{ij}\big{|}_{N}$ and rearranging gives

$\hat{M}^{(k)}_{ij}\big{|}_{N}=\hat{M}^{(k)}_{ij}\big{|}_{N-1}+K_{h}(\mathcal{X}_{j}-X_{N-i})\\\
\times\left(\left[X_{N}-X_{N-i}\right]^{k}-\hat{M}^{(k)}_{ij}\big{|}_{N-1}\right)\Big{/}W_{ij}\big{|}_{N},$
(A5)

as required by (11).

###  .2 Conditional variance

An online calculation of the conditional variance (10) is achieved through
incremental updating of the quantity $S_{ij}\big{|}_{N}$, the weighted sum of
squares of differences from the current mean

$S_{ij}\big{|}_{N}:=\sum\limits_{n=1}^{N-i}K_{h}(\mathcal{X}_{j}-X_{n})\Big{(}\big{[}X_{n+i}-X_{n}\big{]}-\hat{M}^{(1)}_{ij}\big{|}_{N}\Big{)}^{2}.$
(A6)

Derivation of an incremental formula for this expression uses (A2) and (A5),
and follows in a similar fashion to Sec. A.1:

$\displaystyle\begin{aligned}
\mathllap{S_{ij}\big{|}_{N}}&=\Bigg{[}\sum\limits_{n=1}^{N-i}K_{h}(\mathcal{X}_{j}-X_{n})\big{[}X_{n+i}-X_{n}\big{]}^{2}\Bigg{]}-\Big{(}\hat{M}^{(1)}_{ij}\big{|}_{N}\Big{)}^{2}\cdot
W_{ij}\big{|}_{N},\end{aligned}$ (A7) $\displaystyle\begin{aligned}
&=S_{ij}\big{|}_{N-1}+K_{h}(\mathcal{X}_{j}-X_{N-i})\big{[}X_{N}-X_{N-i}\big{]}^{2}+\Big{(}\hat{M}^{(1)}_{ij}\big{|}_{N-1}\Big{)}^{2}\cdot\left(W_{ij}\big{|}_{N}-K_{h}(\mathcal{X}_{j}-X_{N-i})\right)-\Big{(}\hat{M}^{(1)}_{ij}\big{|}_{N}\Big{)}^{2}\cdot
W_{ij}\big{|}_{N},\\\ \end{aligned}$ (A8) $\displaystyle\begin{aligned}
&=S_{ij}\big{|}_{N-1}+K_{h}(\mathcal{X}_{j}-X_{N-i})\bigg{\\{}\big{[}X_{N}-X_{N-i}\big{]}^{2}-\Big{(}\hat{M}^{(1)}_{ij}\big{|}_{N-1}\Big{)}^{2}\bigg{\\}}\\\
&\qquad\qquad\qquad\qquad-
W_{ij}\big{|}_{N-1}\cdot\left(\hat{M}^{(1)}_{ij}\big{|}_{N}-\hat{M}^{(1)}_{ij}\big{|}_{N-1}\right)\left(\hat{M}^{(1)}_{ij}\big{|}_{N}+\hat{M}^{(1)}_{ij}\big{|}_{N-1}\right),\\\
\end{aligned}$ (A9) $\displaystyle\begin{aligned}
&=S_{ij}\big{|}_{N-1}+K_{h}(\mathcal{X}_{j}-X_{N-i})\bigg{\\{}\big{[}X_{N}-X_{N-i}\big{]}^{2}-\Big{(}\hat{M}^{(1)}_{ij}\big{|}_{N-1}\Big{)}^{2}\\\
&\qquad\qquad\qquad\qquad-\Big{(}\big{[}X_{N}-X_{N-i}\big{]}-\hat{M}^{(1)}_{ij}\big{|}_{N-1}\Big{)}\Big{(}\hat{M}^{(1)}_{ij}\big{|}_{N}+\hat{M}^{(1)}_{ij}\big{|}_{N-1}\Big{)}\bigg{\\}},\\\
\end{aligned}$ (A10) $\displaystyle\begin{aligned}
&=S_{ij}\big{|}_{N-1}+K_{h}(\mathcal{X}_{j}-X_{N-i})\left(\left(X_{N}-X_{N-i}\right)-\hat{M}^{(1)}_{ij}\big{|}_{N-1}\right)\left(\left(X_{N}-X_{N-i}\right)-\hat{M}^{(1)}_{ij}\big{|}_{N}\right).\end{aligned}$
(A11)

## References

  * Haken [2004] H. Haken,  _Synergetics: Intoduction and advanced topics_ (Spinger-Verlag, 2004).
  * Risken [1996] H. Risken, in  _The Fokker-Planck Equation_ (Springer, 1996) pp. 63–95.
  * Siegert  _et al._ [1998] S. Siegert, R. Friedrich, and J. Peinke, Physics Letters A 243, 275 (1998).
  * Gottschall and Peinke [2008] J. Gottschall and J. Peinke, New Journal of Physics 10, 083034 (2008).
  * Böttcher  _et al._ [2006] F. Böttcher, J. Peinke, D. Kleinhans, R. Friedrich, P. G. Lind, and M. Haase, Physical Review Letters 97, 090603 (2006).
  * Lind  _et al._ [2010] P. G. Lind, M. Haase, F. Böttcher, J. Peinke, D. Kleinhans, and R. Friedrich, Physical Review E 81, 041125 (2010).
  * Scholz  _et al._ [2017] T. Scholz, F. Raischel, V. V. Lopes, B. Lehle, M. Wächter, J. Peinke, and P. G. Lind, Physics Letters A 381, 194 (2017).
  * Lade [2009] S. J. Lade, Physics Letters A 373, 3705 (2009).
  * Honisch and Friedrich [2011] C. Honisch and R. Friedrich, Physical Review E 83, 066701 (2011).
  * Rydin Gorjão  _et al._ [2021] L. Rydin Gorjão, D. Witthaut, K. Lehnertz, and P. G. Lind, Entropy 23, 517 (2021).
  * Lehle and Peinke [2018] B. Lehle and J. Peinke, Physical Review E 97, 012113 (2018).
  * Friedrich  _et al._ [2011] R. Friedrich, J. Peinke, M. Sahimi, and M. R. R. Tabar, Physics Reports 506, 87 (2011).
  * Tabar [2019] M. R. R. Tabar,  _Analysis and data-based reconstruction of complex nonlinear dynamical systems_ , Vol. 730 (Springer, 2019).
  * Friedrich and Peinke [1997] R. Friedrich and J. Peinke, Physical Review Letters 78, 863 (1997).
  * Renner  _et al._ [2001] C. Renner, J. Peinke, and R. Friedrich, Journal of Fluid Mechanics 433, 383 (2001).
  * Friedrich and Grauer [2020] J. Friedrich and R. Grauer, Atmosphere 11, 1003 (2020).
  * Milan  _et al._ [2013] P. Milan, M. Wächter, and J. Peinke, Physical Review Letters 110, 138701 (2013).
  * Lind  _et al._ [2005] P. G. Lind, A. Mora, J. A. C. Gallas, and M. Haase, Physical Review E 72, 056706 (2005).
  * Buffett  _et al._ [2013] B. A. Buffett, L. Ziegler, and C. G. Constable, Geophysical Journal International 195, 86 (2013).
  * Davis and Buffett [2021] W. Davis and B. Buffett, Geophysical Journal International 228, 1478 (2021).
  * Lamouroux and Lehnertz [2009] D. Lamouroux and K. Lehnertz, Physics Letters A 373, 3507 (2009).
  * Rydin Gorjão and Meirinhos [2019] L. Rydin Gorjão and F. Meirinhos, Journal of Open Source Software 4, 1693 (2019).
  * Kleinhans and Friedrich [2007] D. Kleinhans and R. Friedrich, in  _Wind Energy_ (Springer, 2007) pp. 129–133.
  * Raischel  _et al._ [2014] F. Raischel, A. Moreira, and P. G. Lind, The European Physical Journal Special Topics 223, 2107 (2014).
  * Rinn  _et al._ [2016] P. Rinn, P. G. Lind, M. Wächter, and J. Peinke, Journal of Open Research Software 4 (2016), 10.5334/jors.123.
  * Rydin Gorjão  _et al._ [2023] L. Rydin Gorjão, D. Witthaut, and P. G. Lind, Journal of Statistical Software 105, 1 (2023).
  * Wang  _et al._ [2016] C. Wang, M.-H. Chen, E. Schifano, J. Wu, and J. Yan, Statistics and Its Interface 9, 399 (2016).
  * Karp [1992] R. M. Karp, in  _Algorithms, Software, Architecture: Information Processing 92: Proceedings of the IFIP 12th World Computer Congress_ , Vol. 1 (International Computer Science Institute, Berkeley, 1992) p. 416.
  * Day and Zhou [2020] J. Day and H. Zhou, Journal of Open Source Software 5 (2020), 10.21105/joss.01816.
  * Epanechnikov [1969] V. A. Epanechnikov, Theory of Probability & Its Applications 14, 153 (1969).
  * Härdle  _et al._ [2004] W. Härdle, M. Müller, S. Sperlich, and A. Werwatz,  _Nonparametric and Semiparametric Models_ , Vol. 1 (Springer, 2004).
  * Ragwitz and Kantz [2001] M. Ragwitz and H. Kantz, Physical Review Letters 87, 254501 (2001).
  * Siefert  _et al._ [2003] M. Siefert, A. Kittel, R. Friedrich, and J. Peinke, EPL (Europhysics Letters) 61, 466 (2003).
  * Fuchs  _et al._ [2022] A. Fuchs, S. Kharche, A. Patil, J. Friedrich, M. Wächter, and J. Peinke, Physics of Fluids 34, 101801 (2022).
  * Davis and Buffett [2022] W. Davis and B. Buffett, Physical Review E 106, 014140 (2022).
  * Welford [1962] B. Welford, Technometrics 4, 419 (1962).
  * West [1979] D. West, Communications of the ACM 22, 532 (1979).
  * Mil’shtejn [1975] G. Mil’shtejn, Theory of Probability & Its Applications 19, 557 (1975).
  * Constable and Parker [1988] C. Constable and R. Parker, Journal of Geophysical Research: Solid Earth 93, 11569 (1988).
  * Lhuillier  _et al._ [2013] F. Lhuillier, G. Hulot, and Y. Gallet, Physics of the Earth and Planetary Interiors 220, 19 (2013).
  * Wicht and Meduri [2016] J. Wicht and D. G. Meduri, Physics of the Earth and Planetary Interiors 259, 45 (2016).
  * Fuchs  _et al._ [2017] A. Fuchs, A. Girard, J. Peinke, P. Diribarne, J. Moro, and G. Guelker, in  _16th European Turbulence Conference, 21-24 August, 2017, Stockholm, Sweden_ (European Mechanics Society, 2017).
  * Peinke  _et al._ [2019] J. Peinke, M. R. R. Tabar, and M. Wächter, Annual Review of Condensed Matter Physics 10, 107 (2019).
  * Pèbay [2008] P. P. Pèbay,  _Formulas for robust, one-pass parallel computation of covariances and arbitrary-order statistical moments._ , Tech. Rep. SAND2008-6212 (Sandia National Laboratories (SNL), Albuquerque, NM, and Livermore, CA, 2008).
  * Anvari  _et al._ [2016] M. Anvari, M. R. R. Tabar, J. Peinke, and K. Lehnertz, Scientific Reports 6, 35435 (2016).
  * Lück  _et al._ [2006] St. Lück, C. Renner, J. Peinke, and R. Friedrich, Physics Letters A 359, 335 (2006).
  * Stresing  _et al._ [2011] R. Stresing, D. Kleinhans, R. Friedrich, and J. Peinke, Physical Review E 83, 046319 (2011).
  * Nikakhtar  _et al._ [2023] F. Nikakhtar, L. Parkavousi, M. Sahimi, M. R. R. Tabar, U. Feudel, and K. Lehnertz, New Journal of Physics 25, 083025 (2023).
  * Chan  _et al._ [1982] T. F. Chan, G. H. Golub, and R. J. LeVeque, in  _COMPSTAT 1982 5th Symposium held at Toulouse 1982_ (Springer, 1982) pp. 30–41.
  * Schubert and Gertz [2018] E. Schubert and M. Gertz, in  _Proceedings of the 30th International Conference on Scientific and Statistical Database Management_ (Association for Computing Machinery, New York, 2018) pp. 1–12.
  * Davis [2023] W. Davis, “OnlineMoments.jl,”  (2023), DOI:10.5281/zenodo.8104832.


# Gender mobility in the labor market with skills-based matching models

Ajaya Adhikari1, Steven Vethman1, Daan Vos1, Marc Lenz1, Ioana Cocu1, Ioannis
Tolios1, Cor J. Veenman1,2

###### Abstract

Skills-based matching promises mobility of workers between different sectors
and occupations in the labor market. In this case, job seekers can look for
jobs they do not yet have experience in, but for which they do have relevant
skills. Currently, there are multiple occupations with a skewed gender
distribution. For skills-based matching, it is unclear if and how a shift in
the gender distribution, which we call _gender mobility_ , between occupations
will be effected. It is expected that the skills-based matching approach will
likely be data-driven, including computational language models and supervised
learning methods.

This work, first, shows the presence of gender segregation in language model-
based skills representation of occupations. Second, we assess the use of these
representations in a potential application based on simulated data, and show
that the gender segregation is propagated by various data-driven skills-based
matching models. These models are based on different language representations
(bag of words, word2vec, and BERT), and distance metrics (static and machine
learning-based). Accordingly, we show how skills-based matching approaches can
be evaluated and compared on matching _performance_ as well as on the risk of
_gender segregation_. Making the gender segregation bias of models more
explicit can help in generating healthy trust in the use of these models in
practice.

##  1 Introduction

The skills-based matching approach aims to create a better fit between those
looking for a job (candidates) and those offering a job (employers) by means
of focusing on the skills necessary to be successful at the job. Previous
research has shown that job satisfaction of employees increases, when one’s
skills are well matched with the job activities (Vieira 2005). Next, matching
based on skills has the potential to allow more labor mobility of job seekers.
The focus on skills stimulates the formulation of which tasks the job seekers
can perform, instead of diplomas obtained or professional experiences in
similar functions in a specific sector. The emphasis on what job seekers can
do, relative to what they have previously experienced, allows for more job
opportunities (possibly in a new sector) for job seekers as well as a bigger
pool of applicants for employers, especially in bottleneck occupations. The
potential of this approach has been recognized by various organizations
through different initiatives (Prevoo et al. 2022; Rentzsch and Staneva 2020).

Gender segregation is largely present in the labor market (Bureau of Labor
Statistics 2021; Bettio, Verashchagina, and Camilleri-Cassar 2009; Tomaskovic-
Devey et al. 2006). Previous research has shown that the increase in diversity
in teams is associated with increased innovation (Nathan and Lee 2013; Díaz-
García, González-Moreno, and Jose Sáez-Martínez 2013), effectiveness
(Phillips, Liljenquist, and Neale 2009), productivity (Herring 2009), and fair
decision making (Sommers 2006). Moreover, the increase of diversity in
occupations with a good representation of different groups, can serve social
justice goals for those groups that were historically banned or discouraged to
practice certain occupations (Eagly 2016). Although labor mobility of workers
may increase with skills-based matching, it is unclear whether a shift in
gender distribution between occupations, which we call _gender mobility_ ,
will be effected.

Recruitment using skills-based matching is most probably to be data-driven
(Prevoo et al. 2022) in which computational language models are used to
interpret and represent skills formulation in machine readable vectors. Using
these vectors, machine learning has the opportunity to learn from historical
candidate-job matches which skill sets of candidates fit well with the
required skill sets of job openings.

This work investigates the effect of data-driven skills-based matching models
on the gender segregation in occupations. We first examine whether gender
segregation exists in different skills representations (bag of words (Li et
al. 2020), word2vec (Mikolov et al. 2013), and Bert (Reimers and Gurevych
2019)) of occupations. Second, we measure the impact of biases in skills
representation on gender segregation for the application of skills-based
matching based on simulated data. Typically, these biases are only studied in
the language representations (Delobelle et al. 2022; Goldfarb-Tarrant et al.
2021) without taking the effect on the possible down stream tasks into
account. We examine various skills-based matching models based on the above
mentioned type of representations and different distance metrics (static and
machine learning-based). This allows the evaluation of a model not only
according to its performance but also on its effect on the propagation of
gender segregation.

This paper is structured as follows. Section 2 elaborates on the related work.
In section 3, we describe the dataset and the preprocessing steps that were
done. Section 4 describes the experiments and the outcomes. At last, we
conclude in section 5.

##  2 Related work

To put our contribution into perspective, we first highlight related work on
gender segregation in the context of skills. Since we consider skills
formulations as one of the sources of the segregation as well as the
historical presence of genders in the respective jobs, we also review the
literature on bias measurement with language models.

We follow related work by defining gender segregation as under (over)
representation of a given gender in occupations (Bettio, Verashchagina, and
Camilleri-Cassar 2009). One study showed that the female-dominated occupations
are often associated with lower salary than male-dominated occupations, which
partly stems from low valuation and visibility of female-associated skills
(Grimshaw and Rubery 2007). Academics also point towards the fact that the
type of skills matter. In countries with more focus on general skills rather
than role- or firm-specific skills less occupational gender segregation is
present (Estévez-Abe 2005). Here, general skills are argued to be more gender-
neutral as women have higher concern for skill portability given potential
career interruptions due to family responsibilities. In the U.S., patterns of
skills demand signal an increase in gender segregation that reinforces gender
income inequality (Bettio, Verashchagina, and Camilleri-Cassar 2009). This is
given that the shifts in demand towards lower paid work concern female-
associated skills and shifts towards higher paid work concern male associated
skills. In Britain, evidence was found that skills commonly acquired by women
are increasingly in demand in highly-skilled occupations, such that gender
segregation in those occupations may decrease (Brynin and Perales 2016). All
in all, the interrelation of skills and gender segregation have been
established in related work, where the level and type of skills one has may
impact gender segregation. We contribute by investigating the interplay of
skills and gender segregation within the scope of how risks for this bias can
be measured within AI and language model applications that may foster skills-
based matching.

These language models are based on our written texts which reflect the
stereotypes and human biases that are embedded in our language (Menegatti and
Rubini 2017; Caliskan, Bryson, and Narayanan 2017). The field of trustworthy
AI (Thuraisingham 2022) has therefore established a large knowledge base on
measuring the bias in these models (Hovy and Spruit 2016). We acknowledge the
recent criticism that relevant bias measurement requires the context of the
downstream task in which the language model is applied, i.e. a shift from
measuring language bias in the models themselves towards measuring the impact
language bias has when language models are used in e.g. decision support tools
(Delobelle et al. 2022; Goldfarb-Tarrant et al. 2021). In our contribution, we
therefore measure the risk of language bias exacerbating gender segregation in
the context of the application of skills-based matching.

##  3 Dataset

We position our work in the framework of skills of the O*NET database (O*NET
2022). The O*NET database provides a structured hierarchy of sectors to
occupations and their corresponding required skills. The skills described as
_detailed work activities_ are chosen for the experiments, as they have a
similar level of tangibility in describing the required skills as is observed
in job vacancies.

To assess the gender segregation we also need the gender distribution of the
O*NET occupations. The average gender distribution per O*NET occupation is
retrieved from the U.S. Bureau of Labor Statistics (BLS) (Bureau of Labor
Statistics 2021).

We would, ideally, need real historic matching data between the skills of job
seekers and skills in job ads such that we can train and test skills-based
matching models. As skills-matching is not yet widespread and known
initiatives protect their data for privacy reasons, this data is not
available. For this reason, we gauge the potential impact of skills matching
models on simulated matches of candidates and job openings based on skill
profiles of O*net occupations. Two random sub-sets of skills sampled from the
same occupation are regarded as a good match while two random sub-set of
skills sampled from two different occupations are regarded as a bad match. For
the experiments, we choose to sample sub-sets of size 5. The ability of a
candidate is present in the combination of the different skills as a whole.
One possible way to materialize this in the experiment is by concatenating the
skill sets into one string. This results in a dataset of good and bad matches
between two skills descriptions. 3940 amount of pairs were created with equal
amount of good and bad matches. We choose for the (possible) occurrence of
overlap in sampling of two sub-sets of skills from an occupation; this is also
likely to be the case in reality as people/job openings will overlap in skill.
This set is further equally spit into a train and a test set. For both the
train and test we draw skill profiles from all occupations while the specific
skills of different occupations do not overlap between the train and test set.
This is also the realistic scenario, as the model in practice would also train
and test on skill sets of candidates and job openings related to the same
occupations. To create a fair test set, some sub-sets of skills of a baker,
for example, is trained upon, while the model is tested on other subsets of
the the total skill set of a baker.

##  4 Experiments

This section describes the setup and results of two experiments: (1) presence
of gender segregation in occupations according to their skills representation,
and (2) the evaluation of skills-based matching models according to their
performance and the risk of propagating gender segregation.

For both experiments, we consider three types of vectorizers with different
levels of complexity to convert a string of skills descriptions into one
vector, namely Bag of Words (BoW) (Li et al. 2020), word2vec (Mikolov et al.
2013) and Sentence-BERT (Reimers and Gurevych 2019). BoW counts the presence
of single or combination of words, while word2vec and Sentence-BERT are
language models which encode the semantic meaning of a word in a vector.
Word2vec encodes each word of a sentence independently while Sentence-BERT
also encodes the specific sentence context of a word.

All experiments were run in python 3.9.12. The BoW vectorization is
implemented using 1-gram and 2-gram words in Scikit-learn 1.1.1. library
(Pedregosa et al. 2011). For word2vec, the vectors of different words are
averaged to create one vector using the Spacy 3.4.1. (Honnibal and Montani
2017) library. At last, the sentence transformers 1.2.1. library
(implementation of Sentence-BERT) is used, which can transform a description
of skills into one vector.

### Gender segregation in occupations based on skills

(a) BoW

(b) Word2vec

(c) Sentence-BERT

Figure 1: These figures visualize the current gender segregation in
occupations based on their skills (represented by 3 types of vectorizers)1.

We investigate the current gender segregation in occupations according to
their underlying skills. The list of skills of each O*NET occupation is
concatenated together to create one large string per occupation. These strings
are fed to the three above mentioned vectorizers. The resulting 3 vectors per
occupation are further separately visualized in Figure 1111The interactive web
version of the results of both experiments can be found in
https://fate2022-demo.tnodatalab.nl/#/model by mapping them into two
dimensions using PCA, which encodes the variance in the data as much as
possible. Each point represents an occupation, and the color indicates the
average female gender ratio of that occupation according to BLS.

For all three types of vectorizers, we see clusters of occupations according
to their gender ratio. For example, Sentence-BERT shows a cluster of male
majority occupations on the middle right and a female majority occupations on
top left. The male majority cluster contains mostly technical or labor
intensive occupations such as machinist, dishwasher and roofers, while the
female majority cluster contains mostly healthcare occupations such as dental
assistants, registered nurse and respiratory therapists. Thus, this experiment
indicates that there is a risk that language models interpreting the skills
sets of occupations are influenced by gender segregation.

Vectorizer | Similarity/Distance Metr. | AUC | GSR
---|---|---|---
BoW | Metric learning | 0,90 | 0,60
Cosine | 0,94 | 0,72
Euclidean | 0,88 | 0,59
Word2Vec | Metric learning | 0,91 | 0,68
Cosine | 0,82 | 0,57
Euclidean | 0,83 | 0,60
Sentence-BERT | Metric learning | 0,91 | 0,69
Cosine | 0,94 | 0,70
Euclidean | 0,90 | 0,67
Table 1: This table shows the trade-off between AUC performance and GSR for
three types of vectorizers (BoW, Word2vec and Sentence-BERT embedding) in
combination with three types of similarity/distance metrics (Euclidean
distance, Cosine similarity and Metric learning)1.

### Evaluation of skills-based matching models

In the second experiment, we put the perceived risk for gender segregation in
the language representation of skills from the first experiment into the
relevant context. In particular, we demonstrate how skills-matching models can
be evaluated not only for their matching performance but also for their risk
of propagating gender segregation. We demonstrate this setting with three
types of the above mentioned vectorizers and three types of similarity
measures namely Euclidean distance, cosine similarity and metric learning.
This results in 9 versions of skills-based matching models. While Euclidean
and Cosine metrics are pre-defined, metric learning is a supervised learning
method which learns a task-specific similarity measure from the training data.
The _Information Theoretic Metric Learning_ model with default parameters of
the metric-learn 0.6.2 library(de Vazelhes et al. 2020) is used in this
experiment.

The dataset containing skills description pairs of good and bad matches
described in section 3 is used to train and test the skills-based matching
models. Each pair of skills descriptions of the test set is first vectorized
and the resulting two vectors are further provided as input to a distance or
similarity metric resulting in a matching score. The performance evaluation of
the 9 versions according to the AUC metric can be found in Table 11. We didn’t
further optimize the metric learning algorithm as it is not the focus of this
work.

Figure 2: The figure illustrates the female gender ratio of an occupation and
the average female ratio of 10 occupations with the highest matching score
according to the BoW and cosine similarity version.

To assess the risk of propagating gender segregation by the above mentioned
skills-based matching models, we create a new test-set. Similar to the
training data, a random sub-set of 5 skills of each O*NET occupation is
extracted and concatenated together. Per O*NET occupation the matching score
is computed with the rest of the occupations. For each occupation, we find the
10 occupations with the highest matching score with respect to their skill
set. Then, the overall _Gender Segregation Risk (GSR)_ of a skills-based
matching model is measured by the Pearson correlation between the following
variables: the female gender ratio of an occupation and the average female
ratio of 10 occupations with the highest matching score. An example of the
correlation between these two variables for the BoW and cosine similarity
version is shown in Figure 2. This score provides an indication of how likely
job seekers will end up in a new occupation with similar gender distribution
as their previous occupation if they follow one of the top 10 suggestions made
by the matching model. Table 11 shows the GSR of the different matching
models. Note that in this experiment, it is not an issue that there might be
overlap between the train-set and this new test-set, because the model is not
optimized to minimize the GSR.

We see a positive correlation of higher than 0.5 for all models, indicating
that gender segregation is propagated by these models. We also see that this
risk is correlated with the performance of the model. This is expected
because, for example, a bad model which randomly assigns a matching score,
will randomly suggest the top 10 occupations. The gender ratio of this
suggested occupations will not be correlated with the query occupation.

This allows one to evaluate and compare different models according to its
performance and the GSR. In our simulated application, the combination of BoW
and cosine has the best performance. While the BoW with metric learning
performs slightly less, it has the best trade-off between high performance and
low GSR.

##  5 Conclusion

We measure the risk for gender segregation in skills descriptions and next to
that demonstrate with simulated data the risk of propagating gender
segregation within a potential skills-matching application that uses language
models to interpret skills descriptions. The measured language bias in skills
has shown the need to consider the risk for gender segregation when language
models are put into use for skills-based matching. To facilitate this, our
work provides a first exploration on how quantifying the measured risk for
gender segregation can aid design choices in practice by considering both the
performance and the risk of propagating gender segregation when choosing a
skills-based matching model. The quantification of this bias can contribute
towards generating healthy trust in the use of skills-based matching models.

Future work is needed with real data to validate to what extent this risk is
also present in applications and to establish how this risk can be integrated
in design choices. Moreover, future research is advised to include more
diversity aspects such as ethnicity and educational background.

## Acknowledgements

This work was part of FATE and Skills-matching projects which were funded by
the Appl.AI program within TNO.

## References

  * Bettio, Verashchagina, and Camilleri-Cassar (2009) Bettio, F.; Verashchagina, A.; and Camilleri-Cassar, F. 2009.  Gender segregation in the labour market: Root causes, implications and policy responses in the EU.
  * Brynin and Perales (2016) Brynin, M.; and Perales, F. 2016.  Gender wage inequality: The de-gendering of the occupational structure.  _European Sociological Review_ , 32(1): 162–174.
  * Bureau of Labor Statistics (2021) Bureau of Labor Statistics. 2021.  Labor Force Statistics from the Current Population Survey.  https://www.bls.gov/cps/cpsaat11.htm.  Accessed: 2022-06-10.
  * Caliskan, Bryson, and Narayanan (2017) Caliskan, A.; Bryson, J. J.; and Narayanan, A. 2017.  Semantics derived automatically from language corpora contain human-like biases.  _Science_ , 356(6334): 183–186.
  * de Vazelhes et al. (2020) de Vazelhes, W.; Carey, C.; Tang, Y.; Vauquier, N.; and Bellet, A. 2020\.  metric-learn: Metric Learning Algorithms in Python.  _Journal of Machine Learning Research_ , 21(138): 1–6.
  * Delobelle et al. (2022) Delobelle, P.; Tokpo, E. K.; Calders, T.; and Berendt, B. 2022.  Measuring fairness with biased rulers: A comparative study on bias metrics for pre-trained language models.  In _NAACL 2022: the 2022 Conference of the North American chapter of the Association for Computational Linguistics: human language technologies_ , 1693–1706.
  * Díaz-García, González-Moreno, and Jose Sáez-Martínez (2013) Díaz-García, C.; González-Moreno, A.; and Jose Sáez-Martínez, F. 2013.  Gender diversity within R&D teams: Its impact on radicalness of innovation.  _Innovation_ , 15(2): 149–160.
  * Eagly (2016) Eagly, A. H. 2016.  When passionate advocates meet research on diversity, does the honest broker stand a chance?  _Journal of Social Issues_ , 72(1): 199–222.
  * Estévez-Abe (2005) Estévez-Abe, M. 2005.  Gender bias in skills and social policies: The varieties of capitalism perspective on sex segregation.  _Social Politics: International Studies in Gender, State & Society_, 12(2): 180–215.
  * Goldfarb-Tarrant et al. (2021) Goldfarb-Tarrant, S.; Marchant, R.; Muñoz Sánchez, R.; Pandya, M.; and Lopez, A. 2021.  Intrinsic Bias Metrics Do Not Correlate with Application Bias.  In _Proceedings of the 59th Annual Meeting of the Association for Computational Linguistics and the 11th International Joint Conference on Natural Language Processing (Volume 1: Long Papers)_ , 1926–1940. Online: Association for Computational Linguistics.
  * Grimshaw and Rubery (2007) Grimshaw, D.; and Rubery, J. 2007.  Undervaluing women’s work.  _European Work and Employment Research Centre - University of Manchester_.
  * Herring (2009) Herring, C. 2009.  Does diversity pay?: Race, gender, and the business case for diversity.  _American sociological review_ , 74(2): 208–224.
  * Honnibal and Montani (2017) Honnibal, M.; and Montani, I. 2017.  spaCy 2: Natural language understanding with Bloom embeddings, convolutional neural networks and incremental parsing.  To appear.
  * Hovy and Spruit (2016) Hovy, D.; and Spruit, S. L. 2016.  The social impact of natural language processing.  In _Proceedings of the 54th Annual Meeting of the Association for Computational Linguistics (Volume 2: Short Papers)_ , 591–598.
  * Li et al. (2020) Li, Q.; Peng, H.; Li, J.; Xia, C.; Yang, R.; Sun, L.; Yu, P. S.; and He, L. 2020\.  A survey on text classification: From shallow to deep learning.  _arXiv preprint arXiv:2008.00364_.
  * Menegatti and Rubini (2017) Menegatti, M.; and Rubini, M. 2017.  Gender bias and sexism in language.  In _Oxford Research Encyclopedia of Communication_.
  * Mikolov et al. (2013) Mikolov, T.; Sutskever, I.; Chen, K.; Corrado, G. S.; and Dean, J. 2013.  Distributed representations of words and phrases and their compositionality.  In _Adv. in neural information processing systems_ , 3111–3119.
  * Nathan and Lee (2013) Nathan, M.; and Lee, N. 2013.  Cultural Diversity, Innovation, and Entrepreneurship: Firm-level Evidence from L ondon.  _Economic geography_ , 89(4): 367–394.
  * O*NET (2022) O*NET. 2022.  National Center for O*NET Development.  https://www.onetonline.org/.  Accessed: 2022-06-10.
  * Pedregosa et al. (2011) Pedregosa, F.; Varoquaux, G.; Gramfort, A.; Michel, V.; Thirion, B.; Grisel, O.; Blondel, M.; Prettenhofer, P.; Weiss, R.; Dubourg, V.; Vanderplas, J.; Passos, A.; Cournapeau, D.; Brucher, M.; Perrot, M.; and Duchesnay, E. 2011.  Scikit-learn: Machine Learning in Python.  _Journal of Machine Learning Research_ , 12: 2825–2830.
  * Phillips, Liljenquist, and Neale (2009) Phillips, K. W.; Liljenquist, K. A.; and Neale, M. A. 2009.  Is the pain worth the gain? The advantages and liabilities of agreeing with socially distinct newcomers.  _Personality and Social Psychology Bulletin_ , 35(3): 336–350.
  * Prevoo et al. (2022) Prevoo, T.; Heyma, A.; Künn-Nelen, A.; Montizaan, R.; Steens, S.; and Goedhart, R. 2022.  Inzicht in skills: een verkenning van Nederlandse initiatieven.  _Tijdschrift voor Arbeidsvraagstukken_ , 38(2).
  * Reimers and Gurevych (2019) Reimers, N.; and Gurevych, I. 2019.  Sentence-BERT: Sentence Embeddings using Siamese BERT-Networks.  In _Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing_. Association for Computational Linguistics.
  * Rentzsch and Staneva (2020) Rentzsch, R.; and Staneva, M. 2020.  Skills-Matching and Skills Intelligence through curated and data-driven ontologies.  In _Proceedings of the DELFI Workshops_.
  * Sommers (2006) Sommers, S. R. 2006.  On racial diversity and group decision making: identifying multiple effects of racial composition on jury deliberations.  _Journal of personality and social psychology_ , 90(4): 597.
  * Thuraisingham (2022) Thuraisingham, B. 2022.  Trustworthy Machine Learning.  _IEEE Intelligent Systems_ , 37(1): 21–24.
  * Tomaskovic-Devey et al. (2006) Tomaskovic-Devey, D.; Zimmer, C.; Stainback, K.; Robinson, C.; Taylor, T.; and McTague, T. 2006.  Documenting desegregation: Segregation in American workplaces by race, ethnicity, and sex, 1966–2003.  _American sociological review_ , 71(4): 565–588.
  * Vieira (2005) Vieira, J. A. C. 2005.  Skill mismatches and job satisfaction.  _Economics letters_ , 89(1): 39–47.