# Data-driven intrinsic localized mode detection and classification in one-
dimensional crystal lattice model††thanks: Preprint.

Jānis Bajārs
Faculty of Physics, Mathematics and Optometry
University of Latvia
Jelgavas Street 3
Riga, LV-1004, Latvia
<EMAIL_ADDRESS>&Filips Kozirevs
Faculty of Physics, Mathematics and Optometry
University of Latvia
Jelgavas Street 3
Riga, LV-1004, Latvia
<EMAIL_ADDRESS>
###### Abstract

In this work we propose Support Vector Machine classification algorithms to
classify one-dimensional crystal lattice waves from locally sampled data.
Different learning datasets of particle displacements, momenta and energy
density values are considered. Efficiency of the classification algorithms is
further improved by two dimensionality reduction techniques: Principal
Component Analysis and Locally Linear Embedding. Robustness of classifiers is
investigated and demonstrated. Developed algorithms are successfully applied
to detect localized intrinsic modes in three numerical simulations considering
a case of two localized stationary breather solutions, a single stationary
breather solution in noisy background and two mobile breather collision.

_K_ eywords intrinsic localized modes, discrete breathers, crystal lattice
models, data-driven methods, classification, localization

##  1 Introduction

With constant increase of computational power data-driven methods have gained
a significant popularity in the past few decades with applications spanning
all fields of science [1]. The study of intrinsic localized modes (ILMs) in
crystal lattice models, also known as discrete breathers (DBs), is an active
field of research in mathematics and physics [2, 3, 4, 5] which relies heavily
on obtained data from numerical simulations. Thus, emerging data-driven
methods provide great opportunity to explore new techniques for understanding
nonlinear wave phenomena in crystal lattice models.

Such localized wave solutions are of particular interest in material science
where transport of charge in silicates by moving nonlinear localized modes is
experimentally confirmed, the phenomenon called hyperconductivity [5, 6]. In
addition, localized modes provide mechanisms for the formation of long
decorated dark lines in muscovite mica [4] and can have a functional role in
proteins [7, 8]. Strong evidence of localized phenomena in crystal lattice
models is provided by numerical experiments [9, 10, 11, 12, 13, 14], to
mention but a few.

There are already well established theoretical and numerical understandings of
the existence of one-dimensional stationary DBs [2, 3, 15], i.e., spatially
localized time-periodic excitations. The same may not be said about
propagating DBs [16] and, in particular, in higher dimensions. There are still
open theoretical questions regarding the existence and properties of
propagating DBs in general nonintegrable lattices. Mobile discrete breathers
are not expected to be long-lived, since the lattice models considered are
likely to be nonintegrable and their lifespan is subject to interactions with
lattice defects and the phonon background, which may be viewed as thermal
noise. Thus, an obvious challenge is to study the existence and interactions
of long-lived propagating discrete breather solutions [11, 14] to understand
their role in the transport of energy in sputtering experiments [17]. Such
demanding computations require good data-based numerical tools to detect and
track breather solutions in space and time which serves as a motivation for
our work.

In this work we develop classification algorithms to differentiate between
nonlinear localized waves and nonlocalized linear waves in numerical
simulations of one-dimensional crystal lattice model. Obtained classifiers
rely only on locally sampled data at a given time in opposed to nonlocal time
series analysis methods, such as discrete Fourier and wavelet transforms. For
instance, discrete Fourier transform in space and time provides frequency-
momentum representation of the lattice waves and allows characterization of
exact ILMs [18, 19] but eliminates all space-time domain information and does
not provide any means to locate the localization regions in space at a given
time. On the other hand, localization of frequency and time of a single
particle dynamics can be studied with the Gábor transformation, also known as
short-time Fourier transform, with localization accuracy subject to the
Heisenberg uncertainty relationship, i.e., there are trade-offs between time
and frequency resolutions. Wavelet transforms, which are extensions of the
Gábor transformation, subsequently shorten the scaling window to extract
higher frequencies at better time resolution. Wavelet analysis of ILMs in the
time-frequency domain has been studied in [20, 21], while authors in [22]
considered wavelet analysis in the space-wavenumber domain. Fourier, Gábor and
wavelet transforms provide valuable spectral information and its evolution in
time or space. To determine the actual localization sites at a given time in
this work we propose sliding window object detection method together with
developed classifiers.

Such classifiers can be efficiently trained with different given labeled
datasets which are not only limited to numerical simulation data. Trained
classifiers then can be used to detect localization regions, e.g., in
numerical simulations, which is the first step towards fully automated tool
for quantitative data-based analysis of complex numerical experiments.
Importantly, our analysis and methodology, which follows, extends, in general,
to any one-dimensional crystal lattice models which support intrinsic
localized mode solutions.

The paper is organized in the following way. In Section 2 mathematical model
of one-dimensional crystal lattice is formulated. Spectral properties of
discrete lattice waves are discussed in Section 3. Dimensionality reduction
methods and classification algorithms are presented in Section 4. In Section 5
application of the classification algorithms is demonstrated for detecting
localized modes in numerical simulations. Discussion and conclusions are
provided in Section 6.

##  2 Mathematical model

In this work we are concerned with one-dimensional crystal lattice models of
classical molecular dynamics exhibiting ILM solutions. In particular, we
consider Hamiltonian dynamics of $N$ particles arising from the dimensionless
Hamiltonian:

$H=K+U+V=\sum_{n=1}^{N}\left(\frac{1}{2}\dot{q}_{n}^{2}+U(q_{n})+V(|q_{n+1}-q_{n}|)\right),$
(1)

where $K$ is the kinetic energy, $U$ is the on-site potential energy, $V$ is
the interaction potential of particles, $q_{n}\in\mathbb{R}$ is the position
of the $n^{th}$ particle, $\dot{q}_{n}$ is its time derivative and $|\cdot|$
is the Euclidean distance. In the dimensionless form (1) particle equilibrium
distance $\Delta=q_{n+1}-q_{n}$ as well as all masses are equal to one.

For example, such model (1) arises as one-dimensional model of the K-K layer
of layered silicate muscovite mica [9] or simplified model of two-dimensional
hexagonal crystal lattice [11, 14]. In such case, the on-site potential $U$
models forces from atoms above and below the K-K sheet. Thus, the on-site
potential is modeled as smoothed periodic function [9]:

$U(q_{n})=1-\cos(2\pi q_{n}).$ (2)

Note that more realistic potentials for muscovite mica can also be considered
[18].

There are multiple options for modeling the interaction potential $V$ [23]. In
this work we consider well known scaled Lennard-Jones potential in the
following form:

$V(r_{n})=\epsilon\left(\left(\frac{1}{r_{n}}\right)^{12}-2\left(\frac{1}{r_{n}}\right)^{6}\right),\quad
r_{n}=|q_{n+1}-q_{n}|,$ (3)

where the parameter $\epsilon>0$ describes the ratio of the atomic well depths
associated to the interaction and the on-site potential. That is, large values
of $\epsilon$ lead to a model describing one-dimensional Lennard-Jones fluid,
while small values of $\epsilon$ lead to system of decoupled nonlinear
oscillators. Further simplifications in (1) have been made to only consider
close neighbor interactions with periodic boundary conditions, i.e.,
${q}_{N+1}={q}_{1}$. Note that proposed methodology of Section 4 is not
limited to the particular choice of potentials (2)–(3) and only fixed neighbor
interactions.

From the Hamiltonian (1) obtained Hamiltonian equations:

$\displaystyle\dot{q}_{n}$ $\displaystyle=p_{n},$ (4)
$\displaystyle\dot{p}_{n}$
$\displaystyle=-U^{\prime}(q_{n})+V^{\prime}(|q_{n+1}-q_{n}|)-V^{\prime}(|q_{n}-q_{n-1}|),$
(5)

where $p_{n}$ is momentum, are solved numerically with the second order time
reversible symplectic Verlet method [24]. In the following, without loss of
generality, all analysis and numerical simulations are performed with time
step $\tau=0.01$ and $\epsilon=0.05$.

To excite DB solutions we consider particles in their dynamical equilibrium
states, i.e., $q^{0}_{n}=n-1$ and $p^{0}_{n}=0$, while exciting four
neighboring particle momenta with the pattern:

$\mbox{\boldmath${p}$\unboldmath}_{0}=\gamma\,(-1,2,-2,1)^{T},\quad\gamma>0.$
(6)

Such excitation, due to nonlinearity, generates in our model stationary
discrete breather. In the current study we only consider numerical data from
stationary DB solutions while the extension to differentiate between
stationary and moving breather solutions is left for future work.

In general, larger values of $\gamma$ generate localized waves with larger
particle displacements. While such excitation (6) does not generate exact wave
solutions, the amount of generated linear phonon waves has been observed to be
minimal. Not all values of $\gamma$ can be considered. Too small values will
generate phonons, while too large values will force particles to move out from
their potential wells or in addition to the stationary breather will also
generate moving breathers. Moreover, appropriate values of $\gamma$ depend on
$\epsilon$. With spectral analysis of Section 3 and $\epsilon=0.05$ we find
$\gamma$ lower bound $0.15$ for the initial momenta pattern (6) to generate
stationary DB. On the other hand, to generate arbitrary phonon waves we can
choose small random momenta and particle displacement values as initial
conditions, as confirmed by the spectral analysis in Section 3.

Both described initial conditions will be used to generate datasets in Section
4 to obtain classification and detection algorithms which are capable of
distinguishing between localized and phonon waves. To see discrepancies
between both types of waves in Figure 1 we illustrate particle displacement
function $u_{n}=q_{n}-q_{n}^{0}$ from two numerical simulations with $N=24$ at
final computational time $T_{end}=5$. In Figure 1 we show displacement
function of a localized discrete breather solution computed with $\gamma=0.3$,
while in Figure 1 we demonstrate particle displacements for phonon waves
computed with random initial momentum and particle displacement values. It is
easy to see larger displacement values for the breather solution and that it
is spatially localized compared to the approximately linear waves. In
addition, in Figure 1 we have added markers to eight particle displacement
values on which the breather is localized. Exactly eight data point values
will be used in the classification algorithms of Section 4.

Figure 1: Particle displacements from the numerical simulations at $T_{end}=5$
with $\epsilon=0.05$ and $N=24$. (a) stationary discrete breather solution
initiated with $\gamma=0.3$ and markers indicating eight particles on which
the breather is localized. (b) phonon wave solution which does not exhibit
localization.

##  3 Spectral properties of lattice waves

It is well known that linear phonon waves in crystal lattice models can be
characterized by their dispersion relation [3]. For the model (4)–(5)
linearization at the equilibrium $(q_{n}^{0},p_{n}^{0})=(n-1,0)$ yields the
linearized equations:

$\displaystyle\ddot{u}_{n}$
$\displaystyle=-U^{\prime\prime}(q_{n}^{0})u_{n}+V^{\prime\prime}(r_{n}^{0})u_{n+1}-V^{\prime\prime}(r_{n}^{0})u_{n}-V^{\prime\prime}(r_{n-1}^{0})u_{n}+V^{\prime\prime}(r_{n-1}^{0})u_{n-1}$
$\displaystyle=-U^{\prime\prime}(q_{n}^{0})u_{n}+V^{\prime\prime}(1)(u_{n+1}-2u_{n}+u_{n-1})$
$\displaystyle=-\omega_{0}^{2}u_{n}+c^{2}(u_{n+1}-2u_{n}+u_{n-1}),$ (7)

where $r_{n}^{0}=|q_{n+1}^{0}-q_{n}^{0}|=1$,
$r_{n-1}^{0}=|q_{n}^{0}-q_{n-1}^{0}|=1$,
$\omega_{0}^{2}:=U^{\prime\prime}(q_{n}^{0})=4\pi^{2}$ and
$c^{2}:=V^{\prime\prime}(1)=72\epsilon$. With ansatz of plane wave solutions
$u_{n}=\exp(\mathrm{i}(kn-\omega t))$ we derive the dispersion relation:

$\omega^{2}=\omega_{0}^{2}+4c^{2}\sin^{2}\left(\frac{k}{2}\right),$ (8)

where $\omega$ and $k$ are the phonon frequency and wave number, respectively,
$\omega_{0}$ is the oscillation frequency of the isolated oscillators when
$\epsilon=0$ and $c$ is the speed of sound when $\omega_{0}=0$ and $k\to 0$
which coincides with the maximal phase and group velocity values, i.e.,
$V_{ph}=c\frac{\sin(k/2)}{k/2}$ and $V_{g}=c\cos(k/2)$, respectively. In the
case with the on-site potential $(\omega_{0}\neq 0)$, according to [18], since
$V_{ph}\to\infty$ when $k\to 0$, the more appropriate definition of the sound
speed is the maximum of the group velocity $V_{g}=c^{2}\sin(k)/\omega$ which
is $V_{g}\approx 0.5285$ at $k\approx 0.4753\pi$ if $\epsilon=0.05$.

Figure 2: Contours of two-dimensional discrete Fourier transform in space and
time of numerical simulation with $N=128$ and $T_{end}=400$. Dispersion
relation (8) is indicated by the solid line. (a) stationary discrete breather
amplitude spectrum computed with the initial pattern (6) and $\gamma=0.3$. (b)
phonon wave amplitude spectrum computed with random initial condition and
coinciding with the dispersion curve (8).

The dispersion relation (8) is optic-like [3] since it is bounded from below
by $\omega_{0}=2\pi$ at $k=0$. In addition, it is also bounded from above at
$k=\pm\pi$ as can be seen in Figure 2. In Figure 2 we plot two-dimensional
amplitude spectrum of numerical simulation data in space and time, i.e.,
contours of 2D discrete Fourier transform, where the color axis indicate
Fourier mode amplitude values. We consider a lattice of $N=128$ particles and
integrate in time until $T_{end}=400$. The solid line indicates the dispersion
relation (8).

Excited stationary breather solutions by (6) are high-frequency optical
breathers, meaning that their frequencies are above the optic-like phonon
linear spectrum as can be seen in Figure 2, with adjacent particles moving out
of phase. Only small amount of energy is localized on the dispersion curve
indicating minimal presence of phonons in numerical simulation. It is easy to
see that the stationary DB is not only localized in space, see Figure 1, but
also in frequency, a so called fundamental frequency of an exact time periodic
localized in space wave solution [18], which in our model varies depending on
the value of $\gamma$.

In Figure 2 we considered random initial conditions for the momentum and
displacement variables $p_{n}$ and $u_{n}$, respectively, drawn from the
uniform distribution $\mathcal{U}(-0.01,0.01)$. Obtained momentum values are
then rescaled such that the total energy values (1) are the same in both
simulations. Such random initial condition leads to full spectrum of phonons
which coincides with the dispersion curve (8). Thus, this provides
representative data of phonons for data analysis.

Spectral analysis of this section not only provides insights into the lattice
wave properties but also allows to justify the methods for obtaining so called
labeled training and testing data for classification problem presented in the
following section.

##  4 Classification of lattice waves

In this section we describe classification algorithms of crystal lattice waves
based on locally sampled data. Machine learning techniques [25, 26] are
considered for data dimensionality reduction and building an efficient
classifier. All calculations are performed in Python with built-in functions
of the open source scikit-learn111https://scikit-learn.org/stable/ library.

Figure 3: Energy density $E_{n}+\epsilon$ in time of two numerical simulations
of Figure 2 with $N=128$ and $T_{end}=400$. (a) stationary discrete breather
simulation with the initial pattern (6) and $\gamma=0.3$. (b) phonon wave
simulation with random initial condition.

We differentiate localized discrete breather solutions, see Figure 1, from
nonlocalized linear phonons, see Figure 1. Training (including validation) and
testing data is obtained from total of $N_{sim}=1000$ numerical simulations
with $0.75$ localized and nonlocalized wave data proportion. Local data of
particle displacements, momenta and energies is sampled from eight neighboring
particles indicated by red dots in Figure 1. For presentation purposes we
consider three different datasets for classification: $X_{u}$, $X_{u,p}$ and
$X_{u,p,E}$, where the dataset $X_{u}$ only contains particle displacement
$u_{n}$ values, the dataset $X_{u,p}$ contains particle displacement $u_{n}$
and momenta $p_{n}$ values, while the dataset $X_{u,p,E}$ in addition to
displacement and momenta values contains also particle energy density data
defined by

$E_{n}:=\frac{1}{2}{p}_{n}^{2}+U(q_{n})+\frac{1}{2}\left(V(|q_{n}-q_{n-1}|)+V(|q_{n+1}-q_{n}|)\right).$
(9)

Considering available numerical data, $u_{n}$, $p_{n}$ and $E_{n}$, we have
explored all different combinations, i.e., in total seven different datasets.
Not all results are shown but the main differences are stated. Since we
collect data from eight neighboring particles, i.e., $N_{d}=8$, we have that

$X_{u}\in\mathbb{R}^{N_{sim}\times N_{d}},\quad
X_{u,p}\in\mathbb{R}^{N_{sim}\times 2N_{d}},\quad
X_{u,p,E}\in\mathbb{R}^{N_{sim}\times 3N_{d}}.$ (10)

For further reference we define $K:=jN_{d}$, where $j=1$ or $j=2$, or $j=3$.
Thus, any given dataset is a matrix $X\in\mathbb{R}^{N_{sim}\times K}$.

To appreciate the importance of the particle energy density values (9), which
also will be evident in the following sections, in Figure 3 we illustrate
particle energy density function of two numerical simulations of Figure 2. In
the illustration we have added constant $\epsilon$ to $E_{n}$ such that the
minimal value is equal to zero. It is easy to see that the localized wave in
Figure 3 carries higher value of energy in the localized region compared to
nonlocalized phonon waves in Figure 3. Thus, energy value is an important
quantity to consider in datasets for building a classifier.

In what follows, we excite discrete breather solutions with momenta pattern
(6) where parameter $\gamma$ values are randomly chosen from the uniform
distribution $\mathcal{U}(0.15,1.2)$. Linear phonon waves are generated with
random initial particle displacement and momenta values from the uniform
distribution $\mathcal{U}(-0.01,0.01)$. Data is collected after performing
numerical integration of (4)–(5) in time until $T_{end}=3$ and in lattice with
$N=64$ particles.

###  4.1 Dimensionality reduction

Figure 4: Pearson’s correlation coefficients. (a) correlation matrix of
discrete breather data of the dataset $X_{u,p,E}$. (b) correlation matrix of
phonon data of the dataset $X_{u,p,E}$.

To increase the efficiency of the classification algorithms we consider two
dimensionality reduction techniques, i.e., Principal Component Analysis (PCA)
and Locally Linear Embedding (LLE) [25]. The objective is to find a low-
dimensional representation of high-dimensional correlated datasets (10) and
identify clusters of localized and nonlocalized wave data. In Figure 4 we plot
Pearson’s correlation coefficients of the dataset $X_{u,p,E}$, i.e., the
covariance matrix (11) of all pairs of particle sampled data variables,
$u_{n}$, $p_{n}$ and $E_{n}$, divided by the product of their standard
deviations. Pearson’s correlation coefficient shows linear dependence (linear
regression) between datasets, i.e., positive value indicates that positive
increase in one variable yields also a positive increase in the other
variable, while the negative value indicates a negative increase in the other
variable. In addition, zero value indicates that the variables are not
correlated and the absolute value of coefficients indicates the correlation
strength.

In Figure 4 we consider only discrete breather data from $X_{u,p,E}$ and
observe complex correlations between different particle dynamics data. Notice
only positive correlations between particle energy values, see the right
bottom square of Figure 4, and two squares indicated by dashed lines
suggesting adjacent particle movement out of phase in Figure 1. On the other
hand, in Figure 4 we plot correlation matrix of phonon data of the dataset
$X_{u,p,E}$. Notice very small correlations between sampled particle dynamics
data of phonons.

PCA is a data projection algorithm onto so called principal components that
account for the largest amount of variance in sampled data. The algorithm is
based on the eigen-decomposition of the covariance matrix

$C:=\frac{1}{N_{sim}-1}Y^{T}Y,$ (11)

where $Y:=X-\bar{X}$, $X\in\mathbb{R}^{N_{sim}\times K}$ is given dataset,
such as $X_{u}$ or $X_{u,p}$, or $X_{u,p,E}$, and $\bar{X}$ is the matrix of
the row-wise mean values $\bar{x}\in\mathbb{R}^{1\times K}$, i.e.,

$\bar{x}_{j}:=\frac{1}{N_{sim}}\sum_{i=1}^{N_{sim}}X_{ij},\quad
j=1,\dots,K,\quad\bar{X}:=\begin{pmatrix}1\\\ \vdots\\\
1\end{pmatrix}\bar{x}.$

The covariance matrix (11) is symmetric and semi-positive definite. Thus,
eigen-decomposition yields principal component orthogonal coordinate system
with nonnegative eigenvalues which are squared singular values of the matrix
$Y/\sqrt{N_{sim}-1}$. Lower dimensional representation
$Z\in\mathbb{R}^{N_{sim}\times d}$, where $d\ll K$, is obtained by projecting
data matrix $Y$ onto orthogonal coordinates with associated largest singular
values. In our experiments we observe that with projection down to only two
principal components, i.e., $d=2$, we are able to preserve more than $95\%$ of
dataset’s total variance what is found by computing variance of projected data
onto each principal component.

On the other hand, LLE is a dimensionality reduction technique that computes
low dimensional, neighborhood preserving embeddings of high dimensional data.
In contrast to PCA it does not rely on projections but evaluates linear
relations between closest neighbors and constructs low-dimensional
representation of data by attempting to preserve these relationships. Consider
a dataset $X\in\mathbb{R}^{N_{sim}\times K}$. Mathematically [25] the
algorithm is formulated as two optimization problems. In the first constrained
optimization problem for each data instance
$x^{i}:=X_{i}^{T}\in\mathbb{R}^{K}$, $i=1,\dots,N_{sim}$, i.e., $x^{i}$ is the
transpose of the $i^{th}$ row of the matrix $X$, we search for its $k$ closest
neighbors $x^{j}:=X_{j}^{T}$ and weight values $w_{i,j}$ which give the best
linear representation for the data instance $x^{i}$ of $k$ closest neighbors
$x^{j}$ measured by the square distance error, i.e.,

$\min_{w}\sum_{i=1}^{N_{sim}}\bigg{|}x^{i}-\sum_{j=1}^{N_{sim}}w_{i,j}x^{j}\bigg{|}^{2},$

subject to the following constraints:

$\begin{cases}w_{i,j}=0,&\mbox{if $x^{j}$ is not close neighbor of $x^{i}$ or
$j=i$},\\\ \displaystyle\sum_{j=1}^{N_{sim}}w_{i,j}=1,&\mbox{for
all}\,i=1,\dots,N_{sim},\end{cases}$

where the second constraint imposes invariance under translation.

Obtained matrix $W\in\mathbb{R}^{N_{sim}\times N_{sim}}$ of all weight values
$w_{i,j}$ contains the local linear relationships between the data instances.
We assume that localized wave data in high-dimensional space will be closer
and further away from the linear nonlocalized data points. In our experiments
performing LLE we considered twenty neighbors $(k=20)$ for each data point. In
the second unconstrained optimization problem we minimize the square distance
error to find low-dimensional representation $z^{i}\in\mathbb{R}^{d}$ of
$x^{i}\in\mathbb{R}^{K}$, where $d=2\ll K$ in our results, while preserving
local linear relationships, i.e.,

$\min_{z}\sum_{i=1}^{N_{sim}}\bigg{|}z^{i}-\sum_{j=1}^{N_{sim}}w_{i,j}z^{j}\bigg{|}^{2}.$

Once the low-dimensional dataset $Z\in\mathbb{R}^{N_{sim}\times d}$, i.e.,
$Z_{i}:={z^{i}}^{T}$, is obtained, we subtract row-wise mean values and scale
columns to unit variance before building classification algorithms.

Figure 5: Dimensionality reduction with PCA and LLE of three datasets $X_{u}$,
$X_{u,p}$ and $X_{u,p,E}$ and decision boundaries (12) of SVC indicated by
solid gray curve which separates two prediction regions. Colorbars indicate
$\gamma$ values of the pattern (6) used to obtain DB solution data. (a) PCA of
$X_{u}$ with nonlinear SVC, $\lambda=20$. (b) LLE of $X_{u}$ with nonlinear
SVC, $\lambda=80$. (c) PCA of $X_{u,p}$ with nonlinear SVC, $\lambda=80$. (d)
LLE of $X_{u,p}$ with nonlinear SVC, $\lambda=0.01$. (e) PCA of $X_{u,p,E}$
with linear SVC, $\lambda=1$. (f) LLE of $X_{u,p,E}$ with linear SVC,
$\lambda=0.01$.

In Figure 5 we illustrate dimension reduction of all three datasets $X_{u}$,
$X_{u,p}$ and $X_{u,p,E}$ applying both dimensionality reduction techniques
PCA and LLE, respectively. Colorbars in Figure 5 indicate $\gamma$ values of
the pattern (6) used to obtain particular DB solution data. Apart from the PCA
projection of the datasets $X_{u,p}$, see Figure 5, and $X_{p}$, not shown,
the $\gamma$ values can be traced along the continuous curve in two
dimensions. Interestingly, in Figure 5 the curve is a straight line parallel
to $y$-axis, which was observed for all datsets containing energy values,
i.e., $X_{E}$, $X_{u,E}$, $X_{p,E}$ and $X_{u,p,E}$.

In Figures 5, 5 and 5 we projected datasets down onto two principal components
$(PC1,PC2)$, while in Figures 5, 5 and 5 we show local data linear embedding
in two-dimensional $(LLE1,LLE2)$ coordinates. Notice that linear wave data is
hardly separated from nonlinear localized wave data in Figures 5, 5 and 5. For
the dataset $X_{u,p}$ Figure 5 suggests that LLE is superior over PCA for the
classification problem. Locally linear embeddings of the dataset $X_{u,p,E}$
(shown) and datasets $X_{E}$, $X_{u,E}$ and $X_{p,E}$ (not shown) are able to
produce the largest separation between two types of wave data. Thus,
information of particle displacements and momenta alone may not be sufficient
to build an optimal classifier for detecting ILMs, e.g., in numerical
simulations, see Section 5.

###  4.2 Classification algorithms

Linear and nonlinear kernel Support Vector Machine [25, 26] classifiers (SVCs)
are trained for lattice wave data classification, see Figure 5. SVCs are
supervised learning algorithms that construct a decision boundary given by the
equation:

$w^{T}\phi(z)+b=0,\quad z\in\mathbb{R}^{d},$ (12)

with the largest distance, also called a margin, to the nearest training-data
point of each data class. The function $\phi:\mathbb{R}^{d}\to\mathbb{R}^{D}$
maps $z$ into higher $D$-dimensional feature space $(D\geq d)$.
$w\in\mathbb{R}^{D}$ and $b\in\mathbb{R}$ are unknown parameters to be
determined such that the prediction function
$\text{sgn}\left(w^{T}\phi(z)+b\right)$ for two class classification problem
is correct for most of the training-data instances $z=z^{i}:=Z_{i}^{T}$, where
$Z\in\mathbb{R}^{N_{sim}\times d}$. Special case includes linear classifier
with $\phi(z)=z$. Then, e.g., in two dimensions equation (12) defines a line
as decision boundary, see Figures 5 and 5.

SVC algorithm is formulated as a constrained optimization problem [26]:

$\min_{w,b,\zeta}\left(\frac{1}{2}w^{T}w+\lambda\sum_{i=1}^{N_{sim}}\zeta_{i}\right)$
(13)

with respect to a subject to the constraints

$\begin{cases}y_{i}\left(w^{T}\phi(z^{i})+b\right)\geq 1-\zeta_{i},\\\
\zeta_{i}\geq 0,\end{cases}$ (14)

for all $i=1,\dots,N_{sim}$, where $y_{i}\in\\{-1,1\\}$ is two class, i.e.,
localized and nonlocalized wave, label value for the data instance $z^{i}$.
Constant $\lambda\geq 0$ is the regularization hyperparameter which makes
possible more flexible soft margine classification through the introduction of
nonnegative variables

$\zeta_{i}:=\max(0,1-y_{i}\left(w^{T}\phi(z^{i})+b\right))$

which allows for some data instances to be misclassified, i.e.,

  * •

if $\zeta_{i}=0$, then the data point is correctly classified;

  * •

if $0<\zeta_{i}<1$, then the data point lies inside the margin on the correct
side of the decision boundary (12);

  * •

if $\zeta_{i}=1$, then the data point lies on the decision boundary (12);

  * •

if $\zeta_{i}>1$, then the data point lies on the incorrect side of the
decision boundary (12) and is misclassified.

In order to solve the optimization problem (13)–(14), especially when $D\gg 1$
or even infinite, we define Lagrange multipliers $a_{i}\geq 0$ and
$\mu_{i}\geq 0$ of the constraints (14), respectively, and derive a dual
quadratic programming problem [26, 27]

$\min_{a}\left(\frac{1}{2}\sum_{i=1}^{N_{sim}}\sum_{j=1}^{N_{sim}}a_{i}a_{j}y_{i}y_{j}\mathcal{K}(z^{i},z^{j})-\sum_{i=1}^{N_{sim}}a_{i}\right),$
(15)

subject to the following constraints:

$\begin{cases}\displaystyle\sum_{i=1}^{N_{sim}}y_{i}a_{i}=0,\\\ 0\leq
a_{i}\leq\lambda,&\mbox{for all }i=1,\dots,N_{sim},\end{cases}$ (16)

where $\mathcal{K}(z^{i},z^{j}):=\phi(z^{i})^{T}\phi(z^{j})$ is symmetric
Kernel function which can be defined without explicit knowledge of the
transformation $\phi$. Note that the Lagrange multipliers $\mu_{i}\geq 0$ do
not appear in the optimization problem (15)–(16) but imply that
$a_{i}\leq\lambda$. In addition, if $0<a_{i}<\lambda$, then $\zeta_{i}=0$, and
if $a_{i}=\lambda$, then $\zeta_{i}>0$, while the zero values $a_{i}=0$ do not
contribute to the prediction, since the optimal $w$ of the problem (13)–(14)
is given by

$w=\sum_{i=1}^{N_{sim}}y_{i}a_{i}\phi(z^{i})$

and the prediction function is

$\text{sgn}\left(w^{T}\phi(z)+b\right)=\text{sgn}\left(\sum_{i=1}^{N_{sim}}y_{i}a_{i}\mathcal{K}(z^{i},z)+b\right)$

where $b$ can be determined noting that for all $0<a_{i}<\lambda$ we have that
$\zeta_{i}=0$ which implies that $y_{i}\left(w^{T}\phi(z^{i})+b\right)=1$.
Thus,

$y_{i}\left(\sum_{\begin{subarray}{c}j=1\\\
0<a_{j}<\lambda\end{subarray}}^{N_{sim}}y_{j}a_{j}\mathcal{K}(z^{j},z^{i})+b\right)=1\quad\implies\quad
b=y_{i}-\sum_{\begin{subarray}{c}j=1\\\
0<a_{j}<\lambda\end{subarray}}^{N_{sim}}y_{j}a_{j}\mathcal{K}(z^{j},z^{i})$

where numerically stable solution is obtained by averaging over all indexes
$i$ for which $0<a_{i}<\lambda$.

In addition to linear kernel SVCs,
$\mathcal{K}(z^{i},z^{j}):={z^{i}}^{T}z^{j}$, the Gaussian Radial Basis
Function (RBF) kernel,
$\mathcal{K}(z^{i},z^{j}):=\exp(-\alpha|z^{i}-z^{j}|^{2})$ where $\alpha>0$,
is used in nonlinear SVCs, see Figures 5–5. It showed improved performance
over other nonlinear kernels, such as polynomial and hyperbolic tangent [25].
Soft margin SVCs (15)–(16) are considered where optimal hyperparameter values,
such as $\lambda$ and $\alpha$, can be found with grid search technique [25].
We set problem scaled $\alpha$ value, i.e., $\alpha=1/d=0.5$, while all
$\lambda$ values are found and indicated in Figure’s 5 caption.

To evaluate the performance of linear and nonlinear classifiers for each
dataset $\mathrm{K}$-fold cross-validation [25] was considered with $0.7$
training and testing data proportion, i.e., $70\%$ of random data samples of
both wave types from $N_{sim}$ simulations were considered for training and
cross-validation, while remaining $30\%$ of the data was used for testing. In
our experiments we considered $\mathrm{K}=5$. The precision and recall scores
[25] were recorded on cross-validation and testing datasets. The precision is
a metric of accuracy of the positive (localized wave) predictions while the
recall is a metric of positive instances that are correctly detected by the
classifier.

| Linear kernel | RBF kernel |
---|---|---|---
PCA | LLE | PCA | LLE |
Precision (validation set) | 0.9985 | 0.9623 | 1 | 0.9953 | $X_{u}$
Recall (validation set) | 0.9533 | 0.9478 | 0.9998 | 0.9354
Precision (testing set) | 0.9984 | 0.9658 | 1 | 0.9965
Recall (testing set) | 0.9516 | 0.9487 | 0.9992 | 0.9370
Precision (validation set) | 0.7947 | 0.9927 | 1 | 1 | $X_{u,p}$
Recall (validation set) | 0.9544 | 0.9993 | 0.9754 | 0.9990
Precision (testing set) | 0.7978 | 0.9928 | 1 | 0.9999
Recall (testing set) | 0.9519 | 0.9999 | 0.9758 | 0.9991
Precision (validation set) | 1 | 1 | 1 | 1 | $X_{u,p,E}$
Recall (validation set) | 1 | 1 | 1 | 1
Precision (testing set) | 1 | 1 | 1 | 1
Recall (testing set) | 1 | 1 | 1 | 1
Table 1: Precision and recall metrics averaged over $100$ random train-test-
splits of different classification algorithms, dimensionality reduction
techniques and datasets.

Averaged precision and recall values over $100$ random train-test-splits of
datasets $X_{u}$, $X_{u,p}$ and $X_{u,p,E}$ with $0.7$ data proportion are
listed in Table 1. For each train-test-split we compute precision and recall
scores for the testing data and validation dataset obtained from the
$\mathrm{K}$-fold of training data. Notice that in most cases recall values
are smaller than precision values which is attributed to the fact that more
localized wave data is classified as linear opposed to linear waves being
classified as localized. In addition to mean values we also computed (not
shown) variance values which in all cases, if not exactly zero, did not exceed
$10^{-2}$ value.

Examining Table 1 we observe exceptionally good results except for the PCA of
the dataset $X_{u,p}$ when linear SVC is used, what is also true for the
dataset $X_{p}$. That is also evident in Figure 5, where it is easy to see
that nonlinear waves can not be separated from linear waves by a straight
line. In that case the use of nonlinear SVC is required. Very good results are
demonstrated by PCA and LLE dimensionality reduction of the dataset
$X_{u,p,E}$, including all other datasets containing energy values, either
with linear or nonlinear SVC. In this case for simplicity we advocate the use
of linear SVC. Precision and recall values equal to one state that none of the
data samples were incorrectly classified by SVC. As follows we consider linear
SVC for the dataset $X_{u,p,E}$ and nonlinear classifier for the datasets
$X_{u}$ and $X_{u,p}$, which is illustrated in Figure 5 where the decision
boundaries (12) are indicated by solid gray curves clearly separating two
prediction regions where darker region indicates prediction region of
nonlocalized waves, while the lighter region indicates prediction region of
localized waves. Arguments for using the dataset $X_{u,p,E}$ instead of just
$X_{E}$ will be stated in the following section.

In summary, Table 1 shows that the classifier is not sensitive to random
train-test-split of the datasets and that for the dataset $X_{u,p,E}$ SVC is
essentially a perfect classifier for both dimensionality reduction techniques.
To further differentiate the classifiers we consider their classification
performance on imperfect data, see the following section.

###  4.3 Classification of imperfect data

So far we have considered datasets of perfectly sampled data of linear and
localized nonlinear waves. With application in mind, i.e., detecting ILMs in
numerical simulations, we consider additional datasets $X_{u,p,E}^{n_{off}}$
of localized discrete breather data from $N_{sim}=1000$ simulations with
randomly chosen $\gamma$ values from the uniform distribution
$\mathcal{U}(0.15,1.2)$. We define index $n_{off}$, where
$n_{off}=0,1,\dots,8$, which indicates the shift in particles from which the
data is sampled, i.e., the data is sampled from particles $q_{m+n_{off}}$,
where particle $q_{m}$ is the first particle in Figure 1 identified with red
dot. Thus, the sampled datasets $X_{u,p,E}^{n_{off}}$ consist of partially
linear and nonlinear wave data.

We classify samples from the datasets $X_{u,p,E}^{n_{off}}$ using linear SVC
and both dimensionality reduction methods PCA and LLE. In Figure 6 we plot the
percentage of wave data being classified as localized wave or linear wave
depending on the shift index $n_{off}$. When $n_{off}=0$, the data consists of
purely nonlinear wave data and we can see $100\%$ identification of it. On the
other hand, when $n_{off}=8$ the dataset does not contain information of
nonlinear waves and we see $100\%$ identification of linear wave data. Both
classifiers demonstrate gradual decrease in localized wave identification as
$n_{off}$ value increases. Notice that the curve of LLE drops faster to zero
compared to PCA. This has an impact on detection region size, i.e., SVC with
PCA produces wider detection regions of localization as indicated by numerical
simulations of Section 5.

In addition, in Figure 6 we compare two essentially perfect linear SVCs
trained on the datasets $X_{u,p,E}$ and $X_{E}$ with LLE dimensionality
reduction applied to the datasets $X_{u,p,E}^{n_{off}}$ and $X_{E}^{n_{off}}$.
For both methods we observe desirable trend but the classifier based on the
data consisting only of the energy data will produce much wider detection
regions of the localization. Thus, we postulate that energy data (9) is
important to be included in datasets for training ILM classifiers but may be
not sufficient.

Figure 6: Percentages of classified imperfect localized wave data. (a) linear
SVC with dimensionality reduction methods PCA and LLE of the datasets
$X_{u,p,E}^{n_{off}}$. (b) linear SVC with dimensionality reduction method LLE
of the datasets $X_{u,p,E}^{n_{off}}$ and $X_{E}^{n_{off}}$.

##  5 Applications

In this section we apply classification algorithms developed above to detect
localized DBs of the one-dimensional crystal lattice model (1). We consider
sliding window object detection method by sliding a window of eight particle
size at each time step over the computational data. To detect region of a
localized wave we define a discrete density function $\rho_{n}$,
$n=1,\dots,N$, which takes a positive value if data at particular eight
particles is classified as being of localized wave, otherwise the density
value is equal to zero. In the following Figures 7 and 8 we identify all
particles by red dots where the localization is detected.

For the first example we consider a lattice of $N=64$ particles and excite two
stationary breathers using momenta pattern (6) with $\gamma=0.55$ and
$\gamma=0.35$ at particles with indices $(16,17,18,19)$ and $(44,45,46,47)$,
respectively. We solve the Hamiltonian system (4)–(5) until time $T_{end}=1$
and detect the regions of localization, see Figure 7.

Figure 7: Detection of numerically simulated two stationary breathers excited
by momenta pattern (6) with $\gamma=0.55$ and $\gamma=0.35$ on the left and
right, respectively. Detected localization region particles are identified by
red dots. (a) linear SVC with PCA trained on the dataset $X_{u,p,E}$. (b)
linear SVC with LLE trained on the dataset $X_{u,p,E}$. (c) nonlinear SVC with
PCA trained on the dataset $X_{u}$. (d) nonlinear SVC with LLE trained on the
dataset $X_{u}$.

For this experiment we used four classifiers. In Figure 7 we demonstrate
results of linear SVC with PCA, while in Figure 7 we consider linear SVC with
LLE. Both classifiers are trained on the dataset $X_{u,p,E}$ and are able to
capture regions of localization very well. As previously suggested, see Figure
6, on average SVC classifier with PCA produces wider detection regions, as can
be seen from number of colored particles.

In comparison, nonlinear classifiers trained on the dataset $X_{u}$ containing
only the information of particle displacements, see Figures 7 and 7, are not
always able to capture small displacement ILMs. Similar issues (not shown)
were observed applying nonlinear classifiers trained on the dataset $X_{u,p}$.
Despite these shortcomings, the classifiers can still be applied if additional
information of the energy (9) is unknown.

To further demonstrate methods’ robustness and future research directions we
consider two additional examples. In the first example we consider a
stationary discrete breather in noisy background by initially randomly
exciting all particles and their momenta in addition to the pattern (6) with
$\gamma=0.5$. In this example we performed calculations until time
$T_{end}=10$ and used linear SVC with LLE trained on the dataset $X_{u,p,E}$,
see Figure 8. Figure 8 illustrates very well that the method was able to
detect the regions of localization despite the presence of large amount of
phonon waves. These results suggest that the classification algorithms can be
applied to thermostatic dynamics, i.e., during experiments at fixed
temperature.

Figure 8: Detection of localization regions in numerical simulations with
linear SVC and LLE trained on the dataset $X_{u,p,E}$. Particles of detected
regions are identified by red dots. (a) simulation of stationary discrete
breather ($\gamma=0.5$) in noisy background. (b) simulation of two mobile
breather collision initiated by momenta pattern (17) with $\gamma=0.45$ (left)
and $\gamma=-0.65$ (right).

In the second example we consider a simulation of two mobile breather
collision initiated by three neighboring particle momenta with pattern:

$\mbox{\boldmath${p}$\unboldmath}_{0}=\gamma\,(-1,2,-1)^{T}.$ (17)

The mobile breather propagating to right was excited with $\gamma=0.45$ while
the mobile breather propagating to left was excited with $\gamma=-0.65$.
Integration in time was performed until $T_{end}=100$. Around time $t=50$ both
breathers collide and pass through or reflect from each other and continue to
propagate, see Figure 8. Although the classifier was trained on data obtained
from stationary discrete breathers, we are still able to detect the
localization regions before, during and after the collision. As in Figure 8 we
used linear SVC with LLE trained on the dataset $X_{u,p,E}$. These
satisfactory results show that the methodology can be used beyond just
detection of stationary DBs.

##  6 Discussion and conclusions

The results of Section 5 demonstrate that classification algorithms, such as
SVC, either linear or nonlinear, can be used to detect localization regions in
numerical simulations of ILMs. The strength of the developed approach lies in
the fact that localization regions can be detected from locally sampled data
at a given time. In future research we will look for the method which makes
predictions with the most optimal width of the localization region and
incorporates spectral properties of the lattice waves. For example, in this
work we only considered data sampled from eight neighboring particles, other
options could also be investigated.

Numerical examples demonstrated that our classifiers are also able to detect
mobile DBs but lack the ability to differentiate between mobile and stationary
solutions. To address this we will explore multilabel and multioutput
classification algorithms in future work. In addition, we plan to extend our
work to higher dimensional crystal lattice models and include kink solutions.

For efficiency and visualization purposes with high success we applied two
dimensionality reduction methods before training SVC. The best results were
obtained with LLE on all datasets containing particle energy data, especially
on the dataset $X_{u,p,E}$, such that we could use linear SVC and find more
optimal localization regions in numerical experiments. In future work we plan
to explore, in particular, more other manifold learning techniques.

To demonstrate versatility of our approach we considered different datasets
for training. We found that the best results are obtained if in addition to
particle displacement values particle momenta and energy density values are
also added to the dataset. Nevertheless, classifiers can be built on all the
datasets and overall good performance can be achieved.

## Acknowledgements

J. Bajārs acknowledges support from PostDocLatvia grant
No.1.1.1.2/VIAA/4/20/617.

## References

  * [1] F.J. Montáns, F. Chinesta, R. Gómez-Bombarelli, and J.N. Kutz.  Data-driven modeling and learning in science and engineering.  Comptes Rendus Mécanique, 347:845–855, 2019.
  * [2] R.S. MacKay and S. Aubry.  Proof of existence of breathers for time-reversible or Hamiltonian networks of weakly coupled oscillators.  Nonlinearity, 7(6):1623–1643, 1994.
  * [3] S. Flach and A.V. Gorbach.  Discrete breathers — Advances in theory and applications.  Physics Reports, 467(1):1–116, 2008.
  * [4] J.F.R. Archilla et al., editors.  Quodons in mica: nonlinear localized travelling excitations in crystals, volume 221 of Springer Series in Materials Science.  Springer International Publishing, 2015.
  * [5] F.M. Russell, J.F.R. Archilla, and S. Medina-Carrasco.  Localized waves in silicates. What do we know from experiments?  In C.H. Skiadas and Y. Dimotikalis, editors, 13th Chaotic Modeling and Simulation International Conference, pages 721–734, Cham, 2021\. Springer International Publishing.
  * [6] F.M. Russell, A.W. Russell, and J.F.R. Archilla.  Hyperconductivity in fluorphlogopite at 300 K and 1.1 T.  EPL, 127:16001, 2019.
  * [7] A. Nicolaï, P. Delarue, and P. Senet.  Intrinsic localized modes in proteins.  Sci Rep, 5:18128, 2015.
  * [8] F. Piazza and Y.-H. Sanejouand.  Discrete breathers in protein structures.  Phys. Biol., 5:026001, 2008.
  * [9] Q. Dou, J. Cuevas, J.C. Eilbeck, and F.M. Russell.  Breathers and kinks in a simulated crystal experiment.  Discrete & Continuous Dynamical Systems - S, 4(5):1107–1118, 2011\.
  * [10] J.L. Marín, F.M. Russell, and J.C. Eilbeck.  Breathers in cuprate-like lattices.  Physics Letters A, 281(1):21–25, 2001.
  * [11] J. Bajars, J.C. Eilbeck, and B. Leimkuhler.  Nonlinear propagating localized modes in a 2D hexagonal crystal lattice.  Physica D: Nonlinear Phenomena, 301-302:8–20, 2015.
  * [12] I.A. Shepelev, S.V. Dmitriev, A.A. Kudreyko, M.G. Velarde, and E.A. Korznikova.  Supersonic voidions in 2D Morse lattice.  Chaos, Solitons & Fractals, 140:110217, 2020.
  * [13] A.P. Chetverikov, W. Ebeling, and M.G. Velarde.  Collisions of quasi-one-dimensional solitons in triangular Morse lattice.  Letters on Materials, 6(1):82–85, 2016.
  * [14] J. Bajārs, J.C. Eilbeck, and B. Leimkuhler.  Two-dimensional mobile breather scattering in a hexagonal crystal lattice.  Phys. Rev. E, 103:022212, 2021.
  * [15] S. Aubry.  Discrete breathers: Localization and transfer of energy in discrete Hamiltonian nonlinear systems.  Physica D: Nonlinear Phenomena, 216(1):1–30, 2006.
  * [16] S. Flach and K. Kladko.  Moving discrete breathers.  Physica D: Nonlinear Phenomena, 127:61–72, 1999.
  * [17] F.M. Russell and J.C. Eilbeck.  Evidence for moving breathers in a layered crystal insulator at 300 K.  Europhys. Lett., 78:10004, 2007.
  * [18] J.F.R. Archilla, Y. Doi, and M. Kimura.  Pterobreathers in a model for a layered crystal with realistic potentials: Exact moving breathers in a moving frame.  Phys. Rev. E, 100:022206, 2019.
  * [19] J. Bajārs and J.F.R. Archilla.  Frequency-momentum representation of moving breathers in a two dimensional hexagonal lattice.  arXiv preprint arXiv:2201.03207, 2022.
  * [20] K. Forinash and W.C. Lang.  Frequency analysis of discrete breather modes using a continuous wavelet transform.  Physica D: Nonlinear Phenomena, 123:437–447, 1998.
  * [21] A. Rivière, S. Lepri, D. Colognesi, and F. Piazza.  Wavelet imaging of transient energy localization in nonlinear systems at thermal equilibrium: The case study of NaI crystals at high temperature.  Physical Review B, 99:024307, 2019.
  * [22] K. Hori.  Wavelet analysis of anharmonic self-localized modes.  Journal of the Physical Society of Japan, 62(6):1819–1822, 1993\.
  * [23] M.P. Allen and D.J. Tildesley.  Computer Simulation of Liquids.  Clarendon Press, USA, 1989.
  * [24] E. Hairer, C. Lubich, and G. Wanner.  Geometric numerical integration: structure-preserving algorithms for ordinary differential equations.  Springer Science & Business Media, 2006.
  * [25] Aurélien Géron.  Hands-on Machine Learning with Scikit-Learn, Keras, and TensorFlow.  O’Reilly Media, Inc., 2019.
  * [26] Christopher M. Bishop.  Pattern Recognition and Machine Learning (Information Science and Statistics).  Springer-Verlag, Berlin, Heidelberg, 2006.
  * [27] Chih-Chung Chang and Chih-Jen Lin.  LIBSVM: A library for support vector machines.  ACM Transactions on Intelligent Systems and Technology, 2:27:1–27:27, 2011.  Software available at http://www.csie.ntu.edu.tw/~cjlin/libsvm.