# A dual-space classification scheme of spatial heterogeneities in 2D
monatomic supercooled liquids

Viet Nguyen  Xueyu Song<EMAIL_ADDRESS><EMAIL_ADDRESS>Ames
Laboratory and Department of Chemistry, Iowa State University, Ames, IA, USA

###### Abstract

Understanding the physics of supercooled liquids near glassy transition
remains one of the major challenges in condensed matter science. There has
been long recognized that supercooled liquids have spatially dynamical
heterogeneity Cubuk  _et al._ (2015); Kob  _et al._ (1997); Sastry  _et al._
(1998); Ediger (2000) whose dynamics in some regions of the sample could
potentially be orders of magnitude faster than the dynamics in other regions
only a few nanometers away. However, to identify such domain structures both
structurally and configurationally in a consistent fashion and the connection
between structures and dynamics remains elusive. A new approach to classify
these spatial heterogeneities of supercooled liquids is developed. Using
average weighted coordination numbers (WCNs) of particles as features for the
Principle Component Analysis (PCA) and Gaussian Mixture (GM) clustering, a
representation of the feature space to perform GM clustering is constructed
after applying PCA to WCNs. Nano-domains or aggregated clusters are found in
configurational (real) space by direct mapping the identities of clusters in
the feature (structural) space with some classification uncertainties. These
uncertainties can be improved by a co-learning strategy to transmit
information between the structural space and configurational space iteratively
until convergence. The domains are indeed found to be remarkably consistent in
both spaces over long times as well as measured to have heterogeneous
dynamics. Using such a classification scheme the order parameter time
correlation function follows the non-conserved order parameter scaling law,
which indicates that such domains are the result of liquid-liquid phase
separation after quenching.

##  I Introduction

Phenomena of spatial heterogeneity in supercooled liquids near glass
transition have been studied intensively and widely accepted over past few
decades due to considerable experimental Vidal Russell and Israeloff (2000);
Kegel and van Blaaderen (2000); Marcus  _et al._ (2000) and theoretical Kob
_et al._ (1997); Perera and Harrowell (1996); Doliwa and Heuer (1998) efforts.
Observation of heterogeneous dynamics is closely related to the plateau region
of the relaxation dynamics and the onset of caging process, where particles
with the similar mobility tend to move in cooperatively manner that form
dynamically correlated mesoscopic domains or clusters. Royall and Williams
(2015); Cavagna (2009). Several theories of glass transition have been
developed for a fundamental description of these spatial domains such as the
energy landscape picture Goldstein (1969), the Adam-Gibbs theory Adam and
Gibbs (1965) and the random first-order transition theory Kirkpatrick  _et
al._ (1989). These theories are linked explicitly to thermodynamics of the
system to characterize the formation of domain structures. Despite of their
simple and intuitive description to the spatial domains of liquids, it lacks
reliable methods to characterize them. Meanwhile, some classification
approaches like Voronoi polyhedra, bond-orientational order parameters
Boattini  _et al._ (2020); Spellings and Glotzer (2018); Boattini  _et al._
(2019), the common-neighbour analysis Reinhart  _et al._ (2017), and
topological clusters require some specific structures $\it a$ $\it priori$ of
the system which is unknown in general while some other ‘order-agnostic”
approaches utilize mutual information based on Shannon entropy lacks
microscopic structural details of the system Royall and Williams (2015);
Dunleavy  _et al._ (2015).

In this study, a new classification approach is developed to identify these
spatially heterogeneous domains showing improvement over previous methods. In
our approach, radial distribution function $\it g(r)$ does not play a role in
qualitatively delivering signatures of structural heterogeneities Sastry  _et
al._ (1998) but becomes an order-parameter to unlock the picture of
heterogeneous domains by a collection of shell coordination numbers (CNs). The
collection of CNs is highlighted by two points: collecting from a simulation
frame to frame to avoid the statistically average of structures and employing
a Gaussian weight to overcome a conventional cut-off distance criterion to
scale the contribution of neighboring particles continuously within a shell.
With aids of some machine learning (ML) algorithms such as Principle Component
Analysis (PCA), K-means clustering and Gaussian Mixture (GM), the information
of CNs are fully extracted out to reveal the formation of nano-domains both in
structural and configurational spaces consistently. These nano-domains from
our classification scheme indeed reflect widely accepted spatial
heterogeneities phenomena.

By taking advantage of the classification scheme, nature of these domains is
clearly clarified and the spatially correlated heterogeneous dynamics are
characterized: the intra-domain relaxation is primarily due to the same
mobility of particles inside the domain and inter-domain relaxation which
influences the morphological evolution of nano-domains structures.
Furthermore, the slow inter-domain relaxation is related to the coarsening
kinetics of liquid-liquid phase separation.

In this paper, the method is implemented for a 2D simple monatomic system with
Lennard-Jones Gauss (LJG) interatomic potential Mizuguchi and Odagaki (2009);
Engel and Trebin (2007); Hoang and Odagaki (2011). The reason to select this
2D model is because this LJG potential leads to glass transition which
normally occurs with more complex systems Mizuguchi and Odagaki (2009) and the
formation of nano-domains in 2D can be visualized easily to provide a clearer
picture of the phenomenon than in 3D Nguyen and Song (2022, 2023).

The paper is organized as follows. In Sec. II, our MD simulated model and
proposed methods are described in detail. This is followed by an extensive
results presentation with discussions in Sec. III. Some concluding remarks are
given in Sec. IV.

##  II Methodology

###  II.1 Simulation Details

Fig. 1: The plot of LJG potential as a function of distance between two atoms.

We consider a system of particles which interact isotropically through the LJG
potential whose form is given by:

$V(r)=\epsilon\left[\left(\frac{\sigma}{r}\right)^{12}-2\left(\frac{\sigma}{r}\right)^{6}\right]-1.5\epsilon\left[\frac{\left(r-1.47\sigma\right)^{2}}{0.04\sigma^{2}}\right],$
(1)

where the parameter $\epsilon$ is the potential well depth and $\sigma$ is the
characteristic atomic diameter. The LJG potential is a double-well potential
where the first well is the LJ potential and the second well is represented by
a Gaussian function as shown in Fig. 1. The second minimum is added at a
typical range of second-nearest neighbor coordination distance in close-packed
crystals, which suppress the crystallization of the system. The values of
parameters for the position, depth and width of the Gaussian part are selected
so that glassy states are formed. It was found that the existence of the
second well with chosen parameters favors the formation of a pentagonal local
order and pentagons are essential for stability of a glassy state Mizuguchi
and Odagaki (2009). Hence, the LJG potential is a suitable model for our study
of simple supercooled and glass-forming liquids.

All MD simulations of the system are performed using LAMMPS Thompson  _et al._
(2022) with $\it{NVT}$ ensemble. The cut-off distance is set to 2.5$\sigma$
and periodic boundary conditions are applied to all directions. All units are
in LJ-reduced units in line with parameters of solid Ar Montero de Hijes  _et
al._ (2020): the time unit in term of $\tau=t\sqrt{m\sigma^{2}\over\epsilon}$
and reduced temperature is defined by $T^{*}=T(\frac{\epsilon}{k_{B}})$where
$k_{B}$ is the Boltzmann constant. The Verlet algorithm is employed and MD
time step is $0.005\tau$ or about 10 fs for Ar parameters. Initial atomic
configuration of 2D simple square lattice of 6400 particles with a fixed
density $\rho(N/A)=1.0$ is heated up to a high temperature to obtain a liquid
state and relaxed for 2 x $10^{6}$ time steps. The system is then quenched to
a target temperature $T^{*}$ = 0.35 or 0.2, and is equilibrated again before
collecting configurations at every 100 time steps or $0.05\tau$ in the
production run. The cooling process has been done by linearly decreasing
temperature via re-scaling atomic velocity: $T=T_{0}-\gamma n$, where $\gamma$
is the cooling rate (9.67 x $10^{10}$ K/s if taking Ar parameters) and $n$ is
the number of MD steps. The cooling rate we used is similar to the cooling
rate in the reference Spellings and Glotzer (2018). At this cooling rate which
is effectively a quenching, the model shows a metastable supercooled liquid
state for a long time, while the crystallization is suppressed. For each
temperature, three independent runs are generated to improve statistics and
consistency of the results. These temperatures are selected because: $T^{*}$ =
0.35 is just above the glass transition temperature ($T^{*}_{g}=0.31$ )
predicted by Van Hoang  _et al._ (2015) to be able to observe substantial
change of dynamics while $T^{*}$ = $0.2$ is below the $T^{*}_{g}$ to study the
trend of structural heterogeneity for temperature dependence.

###  II.2 Radial Distribution Function (rdf) and Weighted Coordination Numbers
(WCNs)

To detect hidden heterogeneous structures of a disordered system, order
parameters are required to accurately characterize its spatial local
environment. In one extreme, the radial distribution function g(r) has been
commonly used to describe the normal liquid structure by means of aggregating
coordination numbers (CNs) in all solvation shells, i.e., the relative number
of neighboring particles in a particular surrounding spherical shell averaging
over many configurations of all particles, hence a highly averaged function
that lacks the details to capture the heterogeneous structural features of a
super-cooled liquid. At another extreme, for a particular configuration of the
system either a snapshot of simulations or an experimental image, local
structures of an M particles system can be given by M particles’ local
coordination shell structure. Naturally, a middle ground for a useful
structural description of the system will be to classify these local
structures into a few meso-states, which represent the overall structural
heterogeneity of the system. Then the particles in the same meso-state should
form domains in the configurational space, and these domains from different
meso-states together can tile up the whole system, which provide a
classification scheme both structurally and configurationally. Furthermore,
lifetime of the meso-states in the structural space and domains in the
configurational space is in the order of nano-second timescale which is long
enough to make further analysis possible. The formation of meso-states (PC-
space) or nano-domains (real space) are attributed to the plateau region or
caging process. In the 2(b), the plateau region in diffusion at various
temperatures shows the salient features of the onset of caging processes
observed in supercooled liquids near glass transition Kob and Andersen (1995)
and ranges from 5x$10^{1}$ to 2x$10^{4}$ MD units or converted to 0.1 to 40
(ns) for Ar parameters Schrøder and Dyre (2020). Thus, the lifetime of nano-
domains is indeed longer than or at the order of nanoseconds in the
supercooled liquids.

(a) WCNs based on rdf

(b) MSD trajectories

Fig. 2:  (a) Radial distribution function (rdf) employed by weighted Gaussian
functions for the system at $T^{*}=0.2$. The color and superscript indicate
the WCN’s position and number along the rdf. b) All-particle MSDs as a
function of time at various temperatures

From a snapshot of MD simulations, the CNs of each particle in that
configuration can be calculated. However, the major drawback of conventional
CN-based features is the discontinuity of counting neighboring particles to be
in or out of a particular shell. To alleviate this strict assignment, weighted
coordination numbers (WCNs) Rudzinski  _et al._ (2019), which utilize
normalized Gaussian distributions on the shell structure (2(a)), are employed
to weight the contribution of each surrounding particle based on the
particle’s distance to the central one. WCNs are proved to be a better choice
because it could count the fractional contribution of particles based on the
Gaussian distribution in each shell, hence it smoothes out the discontinuity
among solvation shells. After identifying relevant solvation shells at various
maxima and minima along its rdf, normalized Gaussian weighting functions are
placed at the center of those shells. The width of the Gaussian function is
determined by an area of associated solvation feature and intersecting
tolerance of 0 to 0.25 with neighboring Gaussians is chosen in this study, but
other reasonable choices of width and the size of the overlapping areas of the
Gaussians do not change the consistency of the final results. The dimension of
WCNs vector in each configuration is equivalent to number of shells in the
g(r), N = 19 is shown in the 2(a), where the value of each WCN is the
summation of all weighted surrounding particles within that shell in a single
configuration; other reasonable numbers of shells tested yield the same
results. For each configuration, employing this WCNs implementation, the
features data for all particles in a particular configuration is represented
by a matrix ${\widetilde{\bf X}}$ of MxN, which consists of N WCNs for each of
the M particles. Such collected WCNs are still noisy in a disordered system,
hence a smoothing is achieved by using an average of WCNs:
$\overline{WCN}_{i}=\frac{1}{N_{b}}{\sum_{j}^{N_{b}}WCN_{j}}$, where $N_{b}$
is the number of neighboring particles in each shell plus the particle $i$
itself. Each particle’s N features in the $\overline{WCNs}$ matrix is
constructed independently to describe its own local environment with respect
to its surrounding particles, thus doesn’t correlate to each other to provide
order parameters that reveal the hidden structures.

(a)

(b)

(c)

(d)

Fig. 3:  A 2D snapshot of 2 meso-states of 6400 particles system at
$T^{*}=0.35$, for both direct mapping and co-learning method. Panels a) and c)
are for the PC-space and configurational space from direct mapping while
panels b) and d) present clusters formed in the PC-space and configurational
space after convergence from co-learning. The hat symbol for labelling axes of
the PC space represents the inner product of a particle’s $\overline{WCN}$s
with the PCs basis, in this case for the first three PC components. Blue,
orange-colored particles corresponds to meso-state 1 and meso-state 2,
respectively. X,Y are the atomic coordinates in reduced units.

(a) direct mapping

(b) converged iteration

(c) direct mapping

(d) converged iteration

Fig. 4:  A 2D snapshot of 2 meso-states of 16900 particles system at
$T^{*}=0.35$, for both direct mapping and co-learning method. Panels a) and c)
are for the PC-space and configurational space from direct mapping while
panels b) and d) present clusters formed in the PC-space and configurational
space after convergence from co-learning.

(a) direct mapping

(b) converged iteration

(c) direct mapping

(d) converged iteration

Fig. 5:  A 2D snapshot of 2 meso-states of 48400 particles system at
$T^{*}=0.35$, for both direct mapping and co-learning method. Panels a) and c)
are for the PC-space and configurational space from direct mapping while
panels b) and d) present clusters formed in the PC-space and configurational
space after convergence from co-learning.

###  II.3 Principle Component Analysis and clustering

Principal Component Analysis (PCA) Bishop (2006) is used for dimension
reduction, namely to reduce N particles’ features to a few correlated ones.
The new coordinates (PC’s representation) which are generated from the inner
product of original $\overline{WCNs}$ matrix with PC’s basis (PC-space), are
used for K-means clustering to detect meso-states in the structural space. PCA
is performed by the following steps:

  * •

Obtaining the mean-free data
$\mathbb{X=\widetilde{X}-\langle\widetilde{X}\rangle}$ where the average is
over ${M}$ particles for each component of WCNs.

  * •

Forming the correlation matrix $\mathbb{C=X^{\intercal}X}$, which is $N\times
N$.

  * •

The principle components $\mathbb{u_{i}}$ are obtained after solving the
eigenvalue problem: $\mathbb{Cu_{i}={\bm{\sigma_{i}}}^{2}u_{i}}$. The
eigenvalue $\mathbb{\bm{\sigma_{i}^{2}}}$ measures the variance of the data
along each principle component(PC) $i$. PCA is optimal in term of seeking
small numbers of PCs but maximizing cumulative proportion of variance
explained (PVE) $\mathbb{\bm{\sigma_{i}^{2}}}$ by each principle component. In
other words, the numbers of retained PCs depend on their total PVE such that
the total PVE explains $\geq$ 95$\%$ of total variances presented in
${\widetilde{\bf X}}$.

The PC’s basis is $\mathbb{U=[u_{1},u_{2},..u_{N}]}$ where each
$\mathbb{u_{i}}$ is a collective coordinate with $N$ components corresponding
to the number of features in the data input and the new PC representation
mathematically represented as $\mathbb{Y=U^{\intercal}\widetilde{\bf X}}$.

Despite K-means is well-known for unsupervised classification Scherer  _et
al._ (2015), it requires a good initial guess of possible numbers of clusters
K, hence an implementation of an Elbow convergent test provides a range of
number of clusters K effectively where K = 2 to 4 is sufficient enough to
serve our purpose. An alternative solution could be an implementation of
Bayesian algorithms Bishop (2006). After a careful trial-and-error process, K
= 2 is selected. K-means clustering is implemented to classify PC-
representation (3(a)) into 2 distinct meso-states, then clusters in the
configurational space can also be obtained by mapping identities of particles
in each meso-state directly shown in 3(c). Although a direct mapping of
particles associated with a meso-state in the PC-space also results in a
domain in real space, the classification of interfacial particles among
domains remains uncertain. Initially, minimizing euclidean distance to the
centroid is a measure to assign particles to a meso-state in K-means
clustering, hence K-means works well to determine membership of particles
whose distance are well-separated like core particles of each meso-state.
Meanwhile, assignment of interfacial particles could be inclusive due to small
differences in the distance of interfacial particels to either states. In
other words, the membership of interfacial particles are linearly weighted
combination of both states, so a hard assignment using K-means could lead to
misclassifications. Therefore, an alternative soft clustering algorithm like
Gaussian Mixture is used to improve classification of interfacial particles.
To transfer information between the structural space and the configurational
space consistently, a co-learning strategy is implemented:

  1. 1.

Perform K-means clustering in the PC-space.

  2. 2.

Use the initial knowledge of the clustering from the PC-space (direct mapping)
to perform Gaussian Mixture (GM) classification to soften the hard assignment
of K-means in the configurational space:

    1. (a)

Directly map all particles corresponding to each of PC’s states to Cartesian
coordinate to identify distinct aggregated domains in the real space.

    2. (b)

Build a mixture model of multivariate Gaussian distributions of domains, then
maximizing Gaussian probability to assign the membership of each particle to a
domain .

  3. 3.

Repeat the step number 2 but in the PC space: perform GM in the PC-space from
the clustering knowledge in the configurational space.

  4. 4.

Perform GM classification iteratively in both spaces until converged .

This co-learning strategy is highly non-linear because it gains missed
information from a direct transfer between two spaces and consistently re-
assigned interfacial particles by utilizing a probabilistic clustering. After
a few iterations, all meso-states in both spaces are formed and converged as
shown in 3(b) for PC-space and in 3(d) for configurational space which present
clear evidences of spatial heterogeneity of the system at a supercooled
condition. The consistency of the classification protocol is confirmed by
increasing the system size: nano-domains structure are quickly converged and
formed in both spaces after few iterations in the (Fig. 4 a-d and Fig. 5 a-d.

##  III Results and Discussion

###  III.1 Nature of nano-domains: Statics

(a) partial S(q)

(b) partial g(r)

Fig. 6:  (a) Weighted partial static structure factors S(q). (b) Weighted
partial g(r) decomposition of each meso-state at $T^{*}=0.35$.

To demonstrate non-crystallinity of meso-states, the static structure factors
of two randomly selected domains of the same meso-states for both states and
the total system has been calculated. From the observation of all the peaks
along the scattering pattern S(q) in 6(a), there is no evidence of crystal
structure inside nano-domains.

The meso-states observed in the PC space show bimodality of $\overline{WCNs}$
distribution along each solvation shell in Fig. 7. Fig. 7 explores the bimodal
distribution of $\overline{WCNs}$ along the first four shells; details of more
shells are illustrated in A1 a-f of Appendix A. The presence of a bimodal
distribution reflects a coexistence of two meso-states with different local
structures.

From the bimodality of $\overline{WCNs}$, g(r) for each meso-states are
constructed accordingly which is shown in 6(b). The total $\it g(r)$ (black)
of the whole system can be decomposed into two partial $\it g(r)$ (blue and
orange curves) representing two distinct meso-states smoothed by each
solvation shell. The weighted decomposition of $\it g(r)$ is indeed similar to
the weighted distribution in $\overline{WCNs}$, thus explains how
classification works by classifying unique local structures of the system.

(a) ss1

(b) ss2

(c) ss3

(d) ss4

Fig. 7:  Weighted $\overline{WCNs}$ distributions for the first four solvation
shells of each meso-state and whole system at $T^{*}$ = 0.35, respectively.
(ss stands for the solvation shell).

(a) Density profile

(b) Pressure profile

Fig. 8:  (a) Density profile varying with distance from the center of a patch
of meso-state1 at $T^{*}=0.35$. (b) Normal and tangential components of the
pressure tensor varying with the distance from the center of a patch of meso-
state1. The green line and black line are normal ($P_{N}^{sim}$) and
tangential ($P_{T}^{sim}$) pressures obtained in the simulation, respectively.
The red line is the fit to the normal pressure ($P_{N}$) and the blue line is
the tangential pressure ($P_{T}$) obtained to verify mechanical equilibrium.

Another way to quantitatively interpret these distinguished meso-states is to
compute their density and pressure profile in the configurational space.
Because the domain shapes of meso-state1 are irregular, the thermodynamic
properties of the two meso-states such as density and pressure are calculated
using a circular region inside a domain of meso-state1 and a shell in the
outermost meso-state 2 region as shown in Fig. 8 along one direction, while
other directions give similar results. 8(a) shows the radial atomic density
distribution from the center of a domain of meso-state 1. Meanwhile, the
pressure distribution of the system is calculated as described below: The four
components of the pressure tensor ($\it p_{xx}$, $\it p_{yy}$, $\it p_{xy}$
and $\it p_{yx}$) for each atom are computed in the Cartesian coordinate. The
pressure tensor is then transformed into polar coordinates whose corresponding
components will be ($\it p_{rr}$, $\it p_{\theta\theta}$, $\it p_{r\theta}$
and $\it p_{\theta r}$). It is noted that the order of magnitude of the off-
diagonal terms is the same as the diagonal terms. By default, LAMMPS treats 2D
particles as finite-size 3D spheres, hence actually simulated a pseudo-2D
system whose 3D-sphered particles are enforced to only move in the $\it x$ and
$\it y$ direction while the $\it z$-components of velocities and forces are
zeroed out at every time step. This pseudo-2D model physically explains why
the off-diagonal terms have the same magnitude order comparing to the diagonal
ones. In an exact-2D system, the pressure tensor should have negligible
contribution of off-diagonal terms:

$P(r)=P_{N}(r)\textbf{e}_{r}\textbf{e}_{r}+P_{T}(r)\textbf{e}_{\theta}\textbf{e}_{\theta}$
(2)

where $\textbf{e}_{r}$, $\textbf{e}_{\theta}$ are unit vectors and $\it P_{N}$
and $\it P_{T}$ are the radial or normal and transverse components of the
pressure tensor, respectively. The radial profiles of the components $\it
P_{N}(r)$ and $\it P_{T}(r)$ are obtained by integrating out the angular
degrees of freedom over thin cicular shells extending outwards from the
origin. Fig. 8(b) shows the normal ($\it P_{N}$) and tangential ($\it P_{T}$)
pressure profiles for the circular interface. It is verified that the normal
and tangential profiles statisfy the mechanical
equilibrium,$\bm{\nabla}\cdot\bm{P}=0$, which in polar coordinates is given
by:

$P_{T}(r)=P_{N}(r)+r\frac{dP_{N}(r)}{dr},$ (3)

where the second term is the derivative of the normal pressure with respect to
distance from the center of the patch of meso-state1.

Fig. 8 indicates the formation of interfaces which signifies the co-existence
of two local distinct states. It will be interesting to further characterize
such co-existence by computing the chemical potential in both meso-states.

###  III.2  Nature of the nano-domains: Dynamics

(a) $T^{*}$ = 0.35, Snapshot 1

(b) $T^{*}$ = 0.35, Snapshot 2

(c) $T^{*}$ = 0.35, Snapshot 3

(d) $T^{*}$ = 0.2, Snapshot 1

(e) $T^{*}$ = 0.2, Snapshot 2

(f) $T^{*}$ = 0.2, Snapshot 3

Fig. 9:  2D snapshots of meso-states at different temperatures. The snapshots
are took at 10ns lag time.

(a) $T^{*}=0.2$

(b) $T^{*}=0.35$

(c) D distribution

(d) $C(r,t)$

Fig. 10:  a-b) Particles fluctuations of meso-states at different temperatures
over 6ns. c) Diffusion constants distribution at $T^{*}=0.35$. The dash and
solid lines represent measurements of diffusion constants distributions for
two different domains or patches of a corresponding meso-state. d) Equal time
correlation function $C(r,t)$, where the abscissa is the scaling variable
$r/t^{1/2}$ at $T^{*}$= 0.2.

Armed by the identities of these domains, heterogeneous dynamics is related to
relaxation dynamics of nano-domains which can be measured by their diffusion
constant distributions. Diffusion constants distribution is different for
domains belonging to different meso-states while is the same for domains
associated with the same meso-state. After characterizing nano-domains in the
configurational space, the core particles, which remain static during the
simulation time, are collected in each domain. Those core particles are then
used to obtain mean-square displacements to compute diffusion constants
through the Einstein relation. Final results (10(c) support the picture that
the same diffusion constant distribution of core particles for different
domains associated with the same meso-state and different distribution in
different meso-states, but the long time scale heterogeneous dynamics
reflecting the activated cage escaping need much longer simulations to study.

The stability to form nano-domains is closely related to the onset of cage
processes which is signified by the plateau region of diffusion dynamics
through MSD 2(b). The plateau physically corresponds to the trapping of
particles within the cage. These particles take long times but eventually
mange to escape out of the cage and returns to a normal diffusive motion. The
timescale of the cage process depends on temperatures, thus a quantitative
characterization will need further studies, but some qualitative description
can be made by the observation of different configurational snapshots. Fig. 9
shows three different configurations with 10ns lag time at two different
temperatures to qualitatively examine the timescale of nano-domains. At
$T^{*}$ = 0.35 which is above the $T^{*}_{g}$, Fig. 9 a-c show that the
particles have sufficient thermal energy to be able to escape out the cage and
self-rearrange more freely, thus particles of one domain will likely have
chance to hop to other domains whose shapes are frequently changing and become
unstable. Meanwhile, the shapes of nano-domains are more static due to a
freeze of most of particles at their local cage or domain at $T^{*}$ = 0.2
(below the $T^{*}_{g}$) as illustrated in the Fig. 9 d-f. These phenomena are
supported by quantitatively assessing particles fluctuation of the domains in
Fig. 10 a-b. The particles fluctuate more as the temperature increases because
of higher number of mobile particles.

Furthermore, as the system size increases we would expect that a meso-state in
the PC-space should be able to bifurcate into more domains or patches in the
configurational space to tile up the whole space. 4(d) and 5(d) demonstrate
our prediction. The meso-state 1 splits into more domains in correspondence to
the growing system size, thus the larger the system is, the higher number of
patches each meso-state owns.

Finally, to further confirm the nature of these domains reflects a liquid-
liquid phase separation after quenching, the kinetics of domain growth
following the quenching from a high temperature equilibrium state (normal
liquid) state to a deeply supercooled state is calculated. At the quenched
temperature, the morphology of the system evolved either by nucleation or
spinodal decomposition with two coexisting thermodynamic states represented by
the two meso-states shown in 3(d) (we also observed nucleation scenario by
change the initial density of the system Brickley  _et al._ (2023) at the same
quenched temperature). It has been known that the non-equilibrium domain
growth laws of the system follows a growth law of characteristic domain size
$L(t)\sim t^{1/2}$ for a non-conserved scalar order parameter, which can be
measured by an equal-time correlation function $C(r,t)$ Humayun and Bray
(1991); Gunton  _et al._ (1983); Mazenko (1990). The $C(r,t)$ has form:
$C(r,t)=N^{-1}\sum_{i}\psi_{i}(t)\psi_{i+r}(t)$ where N is the number of
particles, $\psi_{i}$ is the scalar order-parameter which $\psi_{i}$ = +1 for
particles identities of state 1 and $\psi_{i}$ = -1 for particles identities
of state 2. $i+r$ indicates a neighboring particle displaced by a distance $r$
relative to the reference particle $i$, hence the product of
$\psi_{i}(t)\psi_{i+r}(t)$ will be +1 between pair of particles from the same
state and will be -1 otherwise. Since the domain identities of the particles
are known from our classification 10(d) indeed shows the results of the scaing
law as $L(t)\sim t^{1/2}$ using the autocorrelation function.

##  IV Conclusions

The results reported in this study indicate that a general approach can be
developed to understand the existence of spatial heterogeneities both
dynamically and structurally in a supercooled liquid. Such heterogeneity
results from a liquid-liquid phase separation when the system is cooled
rapidly. In our approach, weighted coordinate numbers (WCNs) are collected at
various frames. The collection of WCNs this way provides crucial information
to avoid statistical homogeneity, hence reveals the unique difference in
structure of supercooled liquids in contrast with other works that the
structure of supercooled liquids is only investigated by from the features
along the g(r). These features only reflect the highly averaged WCNs
representation of the system. Local environment of all particles is
statistically the same and each particle will have the same CNs which lacks
enough details to provide any realistic description of the spatial
heterogeneities.

PCA with K-means and GM clustering provides a robust method of classifying
particles into different thermodynamic states where one state forms domains
surrounded by another state in the configurational space. Representations of
partial rdf, bimodal distribution in $\overline{WCNs}$, density and pressure
profiles, and domain growth scaling law clearly support our physical
interpretation of the liquid-liquid phase separation picture. Given the
spatial classification scheme, dynamics of such meso-states shows difference
in the distribution of diffusion constants, hence the dynamical
heterogeneities may be the results of different thermodynamic states after
phase separation. Therefore, our classification scheme paves the way for
further statistical mechanics analysis of supercooled liquids.

##  V Acknowledgement

This work is supported by the Division of Chemical and Biological Sciences,
Office of Basic Energy Sciences, U.S. Department of Energy, under Contact No.
DE-AC02-07CH11358 with Iowa State University.

##  Appendix A Bimodality of $\overline{WCNs}$ distribution

Bimodality of $\overline{WCNs}$ distribution is described in the main text. A1
present examples of bimodal distribution along the fifth to tenth solvation
shells at $T^{*}=0.35$.

(a) ss5

(b) ss6

(c) ss7

(d) ss8

(e) ss9

(f) ss10

A1:  Weighted $\overline{WCNs}$ distribution for other six solvation shells of
each meso-state and whole system at $T^{*}$ = 0.35, respectively. (ss stands
for the solvation shell).

## References

  * Cubuk  _et al._ (2015) E. D. Cubuk, S. S. Schoenholz, J. M. Rieser, B. D. Malone, J. Rottler, D. J. Durian, E. Kaxiras, and A. J. Liu, Physical Review Letters 114, 108001 (2015).
  * Kob  _et al._ (1997) W. Kob, C. Donati, S. J. Plimpton, P. H. Poole, and S. C. Glotzer, Physical Review Letters 79, 2827 (1997).
  * Sastry  _et al._ (1998) S. Sastry, P. G. Debenedetti, and F. H. Stillinger, Nature 393, 554 (1998).
  * Ediger (2000) M. D. Ediger, Annual Review of Physical Chemistry 51, 99 (2000), pMID: 11031277, https://doi.org/10.1146/annurev.physchem.51.1.99 .
  * Vidal Russell and Israeloff (2000) E. Vidal Russell and N. E. Israeloff, Nature 408, 695 (2000).
  * Kegel and van Blaaderen (2000) W. Kegel and A. van Blaaderen, Science (New York, N.Y.) 287, 290 (2000).
  * Marcus  _et al._ (2000) A. H. Marcus, J. Schofield, and S. Rice, Physical Review E 61, 7260 (2000).
  * Perera and Harrowell (1996) D. N. Perera and P. Harrowell, Physical Review E 54, 1652 (1996).
  * Doliwa and Heuer (1998) B. Doliwa and A. Heuer, Physical Review Letters 80, 4915 (1998).
  * Royall and Williams (2015) C. P. Royall and S. R. Williams, Physics Reports 560, 1 (2015), the role of local structure in dynamical arrest.
  * Cavagna (2009) A. Cavagna, Physics Reports 476, 51 (2009).
  * Goldstein (1969) M. Goldstein, The Journal of Chemical Physics 51, 3728 (1969), https://doi.org/10.1063/1.1672587 .
  * Adam and Gibbs (1965) G. Adam and J. H. Gibbs, The Journal of Chemical Physics 43, 139 (1965), https://doi.org/10.1063/1.1696442 .
  * Kirkpatrick  _et al._ (1989) T. R. Kirkpatrick, D. Thirumalai, and P. G. Wolynes, Physical Review A 40, 1045 (1989).
  * Boattini  _et al._ (2020) E. Boattini, S. Marín-Aguilar, S. Mitra, G. Foffi, F. Smallenburg, and L. Filion, Nature Communications 11, 5479 (2020).
  * Spellings and Glotzer (2018) M. Spellings and S. C. Glotzer, AIChE Journal 64, 2198 (2018), https://aiche.onlinelibrary.wiley.com/doi/pdf/10.1002/aic.16157 .
  * Boattini  _et al._ (2019) E. Boattini, M. Dijkstra, and L. Filion, J. Chem. Phys 151, 154901 (2019).
  * Reinhart  _et al._ (2017) W. F. Reinhart, A. W. Long, M. P. Howard, A. L. Ferguson, and A. Z. Panagiotopoulos, Soft Matter 13, 4733 (2017).
  * Dunleavy  _et al._ (2015) A. J. Dunleavy, K. Wiesner, R. Yamamoto, and C. P. Royall, Nature Communications 6, 6089 (2015).
  * Mizuguchi and Odagaki (2009) T. Mizuguchi and T. Odagaki, Physical Review E 79, 051501 (2009).
  * Engel and Trebin (2007) M. Engel and H.-R. Trebin, Physical Review Letters 98, 225505 (2007).
  * Hoang and Odagaki (2011) V. V. Hoang and T. Odagaki,  _The Journal of Physical Chemistry B_ , The Journal of Physical Chemistry B 115, 6946 (2011).
  * Nguyen and Song (2022) V. Nguyen and X. Song, arXiv.org (2022), https://arxiv.org/abs/2203.11989.
  * Nguyen and Song (2023) V. Nguyen and X. Song, J. Phys. Con. Matt. (in press, 2023), https://arxiv.org/abs/2304.03469.
  * Thompson  _et al._ (2022) A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton, Computer Physics Communications 271, 108171 (2022).
  * Montero de Hijes  _et al._ (2020) P. Montero de Hijes, J. R. Espinosa, V. Bianco, E. Sanz, and C. Vega,  _The Journal of Physical Chemistry C_ , The Journal of Physical Chemistry C 124, 8795 (2020).
  * Van Hoang  _et al._ (2015) V. Van Hoang, V. Teboul, and T. Odagaki,  _The Journal of Physical Chemistry B_ , The Journal of Physical Chemistry B 119, 15752 (2015).
  * Kob and Andersen (1995) W. Kob and H. C. Andersen, Physical Review E 51, 4626 (1995).
  * Schrøder and Dyre (2020) T. B. Schrøder and J. C. Dyre, The Journal of Chemical Physics 152, 141101 (2020), https://doi.org/10.1063/5.0004093 .
  * Rudzinski  _et al._ (2019) J. F. Rudzinski, M. Radu, and T. Bereau, The Journal of Chemical Physics 150, 024102 (2019), https://doi.org/10.1063/1.5064808 .
  * Bishop (2006) C. M. Bishop,  _Pattern recognition and machine learning_ (Springer, 2006).
  * Scherer  _et al._ (2015) M. K. Scherer, B. Trendelkamp-Schroer, F. Paul, G. Pérez-Hernández, M. Hoffmann, N. Plattner, C. Wehmeyer, J.-H. Prinz, and F. Noé,  _Journal of Chemical Theory and Computation_ , Journal of Chemical Theory and Computation 11, 5525 (2015).
  * Brickley  _et al._ (2023) J. Brickley, V. Nguyen, and X. Song, “unpublished results,” (2023).
  * Humayun and Bray (1991) K. Humayun and A. J. Bray, J. Phys. A: Math. Gen. 24, 1915 (1991).
  * Gunton  _et al._ (1983) J. D. Gunton, M. S. Miguel, and P. Sahni,  _The Dynamics of First Order Phase Transitions_ , edited by E. C. Domb and J. Lebowitz, Phase Transitions and Critical Phenomena, Vol. 8 (Academic Press, London, 1983).
  * Mazenko (1990) G. F. Mazenko, Phys. Rev. B 42, 4487 (1990).


# Neuro-symbolic Meta Reinforcement Learning for Trading

S I Harini,1 Gautam Shroff,2 Ashwin Srinivasan,1 Prayushi Faldu,3 Lovekesh Vig
2

###### Abstract

We model short-duration (e.g. day) trading in financial markets as a
sequential decision-making problem under uncertainty, with the added
complication of continual concept-drift. We, therefore, employ meta
reinforcement learning via the RL2 algorithm. It is also known that human
traders often rely on frequently occurring symbolic patterns in price series.
We employ logical program induction to discover symbolic patterns that occur
frequently as well as recently, and explore whether using such features
improves the performance of our meta reinforcement learning algorithm. We
report experiments on real data indicating that meta-RL is better than vanilla
RL and also benefits from learned symbolic features.

## Introduction

Deep learning techniques have achieved human-like and even superhuman
performance in a number of arenas: video games, strategy games, robotics, etc.
In many of these arenas, the spectrum of human performance varies widely, from
average to expert. Human traders in financial markets also differ greatly in
skill and performance. The consistent success of expert traders is unlikely to
be due to chance alone; it is more likely that such traders are explicitly or
implicitly relying on patterns in the data they see.

External events in the world clearly affect prices in financial markets. Thus,
accurately forecasting financial prices over medium term, e.g., months, weeks
or even days is a challenging if not hopeless task, even for deep-learning
techniques, see (Shah and Shroff 2021). The impact of external events is
reduced when trading over a shorter time frame, such as sub-second high-
frequency trading as well as, to an extent, intra-day trading. However, intra-
day price variations are the result of complex feedback between buyers and
sellers, as modeled in (Wang 2014). Such feedback results in near chaotic, but
not random, behavior, i.e., the auto-correlation in the price series is not
zero: the price series has ‘memory’. Thus, short-duration price series, while
near chaotic, do seem to exhibit patterns albeit not consistently, which is
what human traders exploit.

It is known that deep-learning techniques can indeed track chaotic systems,
such as the Lorenz equations (Madondo and Gibbons 2018); thus, it is not
unreasonable to expect deep learning could be effective in near chaotic
environments such as financial markets. Indeed, there have been many attempts
at applying machine learning, and, more recently, deep learning to trading,
some of which are mentioned in the next section.

Even so, applying deep learning in financial markets faces challenges: First,
data is scarcer than one may expect - there are only so many days of past
history, and so many financial instruments, even if one assumes all
instruments behave alike. This is far fewer data than, say, text or images on
the internet. Second, price series exhibit patterns, these change over time,
and it is well known that markets continually change their ‘regime’, e.g. from
‘trending’ to ‘mean-reverting’, as well as their volatility.

In this paper we (a) model trading as one of sequential decision-making under
uncertainty, (b) apply deep meta reinforcement learning to make trading
decisions and (c) investigate whether the incorporation of features based on
hand-crafted patterns that are often used by human traders improves
performance. Further, we observe a meta-pattern in such hand-crafted patterns
which we use to automatically learn a large number of similar features using
techniques borrowed from inductive logic programming, and investigate whether
these add to the effectiveness of our meta-RL based trading agent. We present
_preliminary_ results on real data that indicate that both meta reinforcement
learning and logical features, both hand-crafted and learned, are more
effective than vanilla RL or primary price features alone. We conclude with
ideas for future exploration.

## Background/Related Work

### Meta-learning

Meta-learning approaches seek to learn in situations training data is scarce,
either inherently or due to rapid distribution shifts that render older data
less relevant to the current task, as is the case for short-duration trading.
Meta-learning techniques are‘ learning to learn’ by training on many related
tasks, so that performance on similar future tasks is improved. Optimization-
based meta-learning, exemplified by the MAML algorithm (Finn, Abbeel, and
Levine 2017) attempts to learn a good parameter initialization such that a few
steps of gradient descent starting from there are sufficient to adapt rapidly
to a new task even with very few data samples. While MAML and related meta-
learning techniques do apply in reinforcement learning setting, they still
require training on a new task, albeit with limited data. In the case of
trading, where an episode is ideally an entire day, it is of limited use since
a model so adapted could only be used the next day, and without further
adaptation as the day progresses. (Metric-based meta-learning techniques, such
as matching networks (Vinyals et al. 2016) do not easily apply in
reinforcement learning.) The third class of traditional meta-learning
techniques are model-based, wherein adaption on new data takes place within
the activations of a network’s hidden states rather than via any gradient-
based updates. The RL2 algorithm (Duan et al. 2016) which we use here is also
a model-based meta-learning technique. In such a technique, adaption on new
data, in our case new rewards, takes place within the network activations as
rewards arrive, making such a technique most applicable in our scenario.

### Machine-learning for Trading

Recent works applying machine learning to trading based on price signals alone
have also used deep neural networks as well as reinforcement learning: ‘Deep
Momentum Networks’ (Lim, Zohren, and Roberts 2019) as well as ‘Momentum
Transformer’ (Wood et al. 2021) formulate the trading task as one of
suggesting the position to take, e.g., 1 for a long (i.e., buy) position, -1
for a short (i.e., sell) position, and 0 for no position. Exiting a buy/sell
position takes place when a 0 action follows the previous 1 actions, etc.
Neural networks are trained to directly optimize volatility-adjusted expected
returns, adjusted for transaction costs, over a trading period. Transaction
costs are computed by tracking when positions are entered and exited, i.e.,
when actions change from one time step to the next. The former paper uses MLPs
and LSTMs, while the latter uses transformers. Both works are essentially
applying vanilla REINFORCE to the MDP formulation of the trading problem. Deep
Reinforcement Learning in Trading (Zhang, Zohren, and Roberts 2020) uses the
same formulation as the above two works, but applies more refined
reinforcement learning techniques, e.g., policy-gradient, actor-critic, and
deep-Q-learning algorithms. In contrast to the formulation used in all the
above three works, we also model actions as buy (1), sell(-1), or do nothing
(0), but these can _only_ be taken when no position (i.e., zero shares) are
held. As soon as a position is taken our training environment computer when
this position exits due to the pre-defined stop-loss/target being met or the
day end being reached, at which point a reward is returned to the RL agent. We
postulate that such a formulation makes for easier learning since the agent
only needs to deal with one kind of situation, i.e. when it holds no position.
The downside is that exit conditions (i.e., stop-loss/targets/end-of-day) are
fixed in advance, rather than determined based on price movements. Note that
however, such exit conditions can also be outputs of the policy, i.e., varying
stop-loss/targets for each buy or sell, based on current volatility, or
whatever the agent finds useful; nevertheless, we have not reported
experiments with this enhancement here.

### Inductive Logic Programming

Inductive logic programming (ILP) (Muggleton and Raedt 1994) investigates the
inductive construction of first-order clausal theories from examples and
background knowledge. ILP is ordinarily employed in a supervised learning
setting, where positive (and usually negative) examples of a target concept
are given in terms of base features. Also supplied is background knowledge in
the form of facts as well as logical rules, typically in a logic programming
language, i.e., Prolog. The ILP process involves constructing a theory that
explains the examples provided with the desired accuracy, support, and
confidence. At each stage in this process, possible theories are tested
against target examples via resolution using background knowledge.

In our case we do not have target examples or concepts, instead, price data is
translated into Prolog facts, and feature _templates_ , or ‘meta-rules’, are
added as background knowledge using techniques introduced in (Muggleton, Lin,
and Tamaddoni-Nezhad 2015). Thereafter, starting from instances selected
randomly as in (Muggleton and Raedt 1994), features with high support in the
data are discovered as in (Dehaspe and Raedt 1997), via (SLD)
resolution111https://en.wikipedia.org/wiki/SLD_resolution using the supplied
background facts and meta-rules.

## Methodology

### Task Formulation and Learning Environment

Each task represents a trading day for a particular symbol (i.e., stock). Data
arrives each minute with the open, high, low, and close prices for the past
minute along with technical analysis indicators (as will be detailed below in
a subsequent section). The agent issues buy/sell/do-nothing actions based on
the data seen so far for the day. Once an order (i.e., buy or sell) is placed,
the agent does not see any data (or reward) until the price moves by an amount
determined by pre-defined stop-loss or target values; e.g. if these are each
1%, the agent receives a reward (positive or negative) when the close price
changes by 1% from the point at which the order was placed, At this point, the
agent resumes receiving data every minute until it places another order.
Alternatively, if the day ends, the agents receive a reward based on the final
closing price of the day. As noted earlier, this formulation differs from that
in prior works (Lim, Zohren, and Roberts 2019; Wood et al. 2021; Zhang,
Zohren, and Roberts 2020).

### Meta-reinforcement Learning: RL2

To deal with continual distribution shift, we employ the meta reinforcement
learning approach RL2 from (Duan et al. 2016). In a standard RL formulation
the agent predicts the next action based on the current state (or history of
states, in the case of a recurrent network) and subsequently receives a
reward. In RL2, the previous action and reward are also input to the network,
and a recurrent network is used. As a result, changes in the
{state,action,reward} distribution are visible to the agent as it takes
actions. Note that the agent is trained over many trials where it encounters
sequences of tasks with possibly different {state, action, reward}
distributions. Thus, in principle, the agent can learn to adapt to a new
distribution when encountering a new task.

The RL2 agent is trained on past data comprising of trials, where each trial
is a sequence of tasks. In the trading context, this entails training the
agent over many day-symbol combinations, and then testing for a new day (for
one of the symbols already seen; though unseen symbols could also be used -
here we use seen symbols). In principle, the meta-RL agent should rapidly
adapt to the reward pattern it experiences in its first few orders even if
these differ from the recent past. We use PPO to train the meta-RL agent and
an LSTM agent (whereas (Duan et al. 2016) used TRPO and a GRU agent), since
PPO is known to be more stable while training and LSTMs are more expressive
than GRUs.

### Hand-crafted Features

Human traders use price patterns as signals on which to base their trading
decisions. An example of two such patterns are depicted in Figure 1. The
‘three crows’ pattern involves three successive observations over which both
opening prices and closing prices decrease sequentially. Similarly, the ’four
horsemen’ pattern involves a sequence of four rising open/close prices.

Figure 1: Example of hand-crafted features

In order to explore whether such features add additional value we detect the
presence or absence of such patterns and append two Boolean features to the
state at each time step to indicate if and which of these holds: thus these
features would be $(0,0)$ for all time steps except at time step 17 where a
$(0,1)$ would indicate the presence of ‘four horsemen’, and at 6,7, and 8
where $(1,0)$ would indicate the presence of ‘three crows’ over three recent
steps.

### Learned Logical Features

The two handcrafted features above are based on patterns involving
increasing/decreasing sequences of primary features, viz., open and close
prices respectively. We postulate that increasing and decreasing sequences of
other primary features might also form useful features. For example,
increasing/decreasing sequences of highs, lows, or even moving averages or
other technical indicators.

Further, it is also possible that increasing/decreasing sequences of _derived_
features, in particular, differences between primary features may be useful.
For example, whether the difference between open and close is narrowing or
widening may be indicative of decreasing or increasing volatility, which in
turn should be useful for deciding an appropriate action. The same can be said
for differences in highs and lows, opens and highs, etc. It is also reasonable
to consider differences between technical indicators as well, e.g.,
differences between moving averages of different lengths; the crossing of such
moving averages is known to be used by traders, so narrowing differences would
point to a possible impending crossing.

Of course, exhaustively enumerating all such differences for each time step
would be both computationally inefficient and likely to confuse the neural
network model. Instead, we use techniques borrowed from inductive logic
programming to enumerate such features in a principled manner and filter these
first based on the frequency of occurrence in the training data, and then by
importance using standard feature-importance determinants.

Meta-patterns are defined in Prolog to capture the concept of ‘runs’, i.e., a
sequence of continually increasing/decreasing values (of one or more
primary/derived features) of arbitrary length. Thereafter randomized search
followed by resolution in Prolog is used to enumerate frequently occurring
patterns (i.e., those that occur in the training data more often than a given
support value). Background clauses are included to define the concept of a
‘derived’ feature as the difference between two primary features. Randomness
is used in the search process to select features to test for. Search proceeds
until the number of high support patterns found reaches a pre-determined
maximum limit.

The above procedure yields a very large number of patterns, which are used to
augment the price data with pattern-based features determined by the presence
or absence of one or more patterns at each time step. These features are used
to train a random forest regression model to predict reward; as a side-effect,
the random forest model returns the importance of each feature towards
predicting reward. These importance values are used to rank pattern-based
features. Finally, a small number of top-ranking features are used to augment
the meta-RL agent’s neural network model.

## Evaluation

### Data

Data received by the agent at each time step is price data, i.e., open, high,
low, close, and volume. These are normalized by dividing the volume column by
the first non zero volume at the beginning of the episode and the rest by
dividing by the close price at beginning of the episode. These normalized
values are then used to compute technical analysis indicators such as simple
moving averages, relative strength indicator,s etc. We use a total of 15 such
technical analysis features in addition to open, high, low, close prices and
volume.

### Results

The agent is trained on $n$ symbols for $m$ days and test scores on the day
$m+1$ are calculated. The results averaged over multiple non-overlapping
subsets of symbols and days are as shown in Table 1 and Table Results. Each
table entry indicates the % average daily return achieved on the test day when
trained using the data of the given number of symbols and previous days. Table
1 shows the performance of vanilla reinforcement learning (i.e., wherein the
previous action and reward is _not_ fed back into the neural network, vs meta
reinforcement learning, and Table Results shows the performance when using
different sets of features as inputs to the RL-agent’s neural network.

|  | 1 symbol | 3 symbols | 6 symbols
---|---|---|---|---
Meta RL | 5 days | 0.26 | -0.23 | 0.36
10 days | -0.06 | 0.11 | 0.27
15 days | 0.11 | 0.04 | 0.09
Vanilla RL | 5 days | -0.12 | -0.33 | -0.45
10 days | -0.17 | -0.31 | -0.52
15 days | -0.27 | -0.22 | -0.38
Table 1: Avg. % returns on test day using vanilla/meta RL.

|  | 1 symbol | 3 symbols | 6 symbols
---|---|---|---|---
Technical Indicator  features | 5 days | 0.26 | -0.23 | 0.36
10 days | -0.06 | 0.11 | 0.27
15 days | 0.11 | 0.04 | 0.09
+Handcrafted features | 5 days | 0.32 | -0.16 | 0.41
10 days | 0.07 | 0.21 | 0.19
15 days | 0.15 | 0.12 | 0.14
+Learned logical features | 5 days | 0.31 | 0.27 | 0.45
10 days | 0.17 | 0.39 | 0.21
15 days | -0.1 | 0.23 | 0.11
+Handcrafted +logical features | 5 days | 0.33 | 0.36 | 0.29
10 days | 0.19 | 0.29 | 0.25
15 days | -0.04 | 0.27 | 0.05

Table 2: Avg. % returns on test day for different feature sets.

## Discussion

We draw the following indicative conclusions from the results presented above:

  1. 1.

Meta reinforcement learning improves over vanilla reinforcement learning,
indicating that distribution shift may be impacting the former.

  2. 2.

Logical features, both hand-crafted as well as learned, improve performance vs
using primary features alone.

  3. 3.

Learned logical features add value over and above hand-crafted features alone.

  4. 4.

Training using too much past data (15 days) is inferior to training using a
moderate amount of past data (5 or 10 days). This may be further evidence of
distribution shift over longer time periods.

  5. 5.

Training on more symbols is better, indicating that distribution shifts take
place more over days rather than across different symbols.

## Conclusions and Future Work

We submit that meta reinforcement learning is a promising direction to explore
for building trading agents using deep learning. Also, logical features
learned using meta-patterns inspired by hand-crafted features may be useful.

Many recent advances in deep learning are worthwhile exploring in the context
of building trading agents: Language models have proven to be few-shot
learners even for numerical data expressed symbolically (Hegselmann et al.
2022); could logical features as we have used here form the basis for
exploring whether these systems could be applied in the trading arena?
Recently de-noising diffusion models have also been used for planning (Janner
et al. 2022); trying such approaches in trading may be worth exploring as
well.

## References

  * Dehaspe and Raedt (1997) Dehaspe, L.; and Raedt, L. D. 1997.  Mining association rules in multiple relations.  In _International Conference on Inductive Logic Programming_ , 125–132. Springer.
  * Duan et al. (2016) Duan, Y.; Schulman, J.; Chen, X.; Bartlett, P. L.; Sutskever, I.; and Abbeel, P. 2016.  Rl 2: Fast reinforcement learning via slow reinforcement learning.  _arXiv preprint arXiv:1611.02779_.
  * Finn, Abbeel, and Levine (2017) Finn, C.; Abbeel, P.; and Levine, S. 2017.  Model-agnostic meta-learning for fast adaptation of deep networks.  In _International conference on machine learning_ , 1126–1135. PMLR.
  * Hegselmann et al. (2022) Hegselmann, S.; Buendia, A.; Lang, H.; Agrawal, M.; Jiang, X.; and Sontag, D. 2022\.  TabLLM: Few-shot Classification of Tabular Data with Large Language Models.  _arXiv preprint arXiv:2210.10723_.
  * Janner et al. (2022) Janner, M.; Du, Y.; Tenenbaum, J. B.; and Levine, S. 2022.  Planning with Diffusion for Flexible Behavior Synthesis.  _arXiv preprint arXiv:2205.09991_.
  * Lim, Zohren, and Roberts (2019) Lim, B.; Zohren, S.; and Roberts, S. 2019.  Enhancing time-series momentum strategies using deep neural networks.  _The Journal of Financial Data Science_ , 1(4): 19–38.
  * Madondo and Gibbons (2018) Madondo, M.; and Gibbons, T. 2018.  Learning and modeling chaos using lstm recurrent neural networks.  In _MICS 2018 proceedings_.
  * Muggleton and Raedt (1994) Muggleton, S.; and Raedt, L. D. 1994.  Inductive Logic Programming: Theory and Methods.  _J. Log. Program._ , 19/20: 629–679.
  * Muggleton, Lin, and Tamaddoni-Nezhad (2015) Muggleton, S. H.; Lin, D.; and Tamaddoni-Nezhad, A. 2015.  Meta-interpretive learning of higher-order dyadic datalog: Predicate invention revisited.  _Machine Learning_ , 100(1): 49–73.
  * Shah and Shroff (2021) Shah, V.; and Shroff, G. 2021.  Forecasting Market Prices using DL with Data Augmentation and Meta-learning: ARIMA still wins!  In _I (Still) Can’t Believe It’s Not Better! NeurIPS 2021 Workshop_.
  * Vinyals et al. (2016) Vinyals, O.; Blundell, C.; Lillicrap, T.; Wierstra, D.; et al. 2016.  Matching networks for one shot learning.  _Advances in neural information processing systems_ , 29.
  * Wang (2014) Wang, L.-X. 2014.  Dynamical Models of Stock Prices Based on Technical Trading Rules—Part II: Analysis of the Model.  _IEEE Transactions on Fuzzy Systems_ , 23(4): 1127–1141.
  * Wood et al. (2021) Wood, K.; Giegerich, S.; Roberts, S.; and Zohren, S. 2021.  Trading with the Momentum Transformer: An Intelligent and Interpretable Architecture.  _arXiv preprint arXiv:2112.08534_.
  * Zhang, Zohren, and Roberts (2020) Zhang, Z.; Zohren, S.; and Roberts, S. 2020.  Deep reinforcement learning for trading.  _The Journal of Financial Data Science_ , 2(2): 25–40.