# Selecting relevant structural features for glassy dynamics by information imbalance

Anand Sharma, Chen Liu, Misaki Ozawa

► To cite this version:

Anand Sharma, Chen Liu, Misaki Ozawa. Selecting relevant structural features for glassy dynamics by information imbalance. The Journal of Chemical Physics, 2024, 161 (18), pp.184506. 10.1063/5.0235084. hal-04895913

HAL Id: hal-04895913

https://hal.science/hal-04895913v1

Submitted on 19 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Selecting Relevant Structural Features for Glassy Dynamics by Information Imbalance

Anand Sharma,<sup>1,2</sup> Chen Liu,<sup>3</sup> and Misaki Ozawa*<sup>2</sup>

<sup>1)</sup>Indian Institute of Science Education and Research, Dr. Homi Bhabha Road, Pashan, Pune 411008, India

<sup>2)</sup>Univ. Grenoble Alpes, CNRS, LIPhy, 38000 Grenoble, France

<sup>3)</sup>Innovation and Research Division, Ge-Room Inc., 93160 Noisy le Grand, France

(*Electronic mail: * misaki.ozawa@univ-grenoble-alpes.fr)

(Dated: 14 November 2024)

We investigate numerically the identification of relevant structural features that contribute to the dynamical heterogeneity in a model glass-forming liquid. By employing the recently proposed information imbalance technique, we select these features from a range of physically motivated descriptors. This selection process is performed in a supervised manner (using both dynamical and structural data) and an unsupervised manner (using only structural data). We then apply the selected features to predict future dynamics using a machine learning technique. One of the advantages of the information imbalance technique is that it does not assume any model *a priori*, i.e., it is a non-parametric method. Finally, we discuss the potential applications of this approach in identifying the dominant mechanisms governing the glassy slow dynamics.

## I. INTRODUCTION

Glass-forming liquids exhibit a significant slowing down in dynamics as the temperature decreases. Relaxation times and transport coefficients, like viscosity, increase by more than ten orders of magnitude within a small temperature range<sup>1-3</sup>. In addition to growing timescales, glass-forming liquids demonstrate heterogeneous spatial patterns of dynamics composed of mobile and immobile domains, called dynamical heterogeneity, which is another hallmark of glassy dynamics<sup>4,5</sup>. The magnitude of dynamical heterogeneity increases with decreasing temperature, implying the emergence of highly non-trivial spatial correlations at lower temperatures.

It has been reported that such heterogeneous dynamical patterns are statistically reproducible if one looks at the mobility or dynamic propensity of each particle obtained from an ensemble of many independent trajectories starting from a given static configuration, called the iso-configurational ensemble<sup>6-8</sup>. Such numerical observations suggest the existence of a relationship between the initial structure and future dynamics. Frustration-based theories attribute glassy dynamical behaviors to the existence of domains with locally preferred structures<sup>9</sup>. Indeed, extensive numerical and experimental works revealed that locally preferred structures are highly correlated with slow dynamics and dynamical heterogeneities<sup>10-14</sup> (see e.g., Ref. 15 and 16 for review).

Besides locally preferred structures, various physically-motivated descriptors characterizing local structural environments (e.g., local energy and density<sup>13,17</sup>, geometrical order parameters<sup>18,19</sup>, mesoscopic elastic stiffness<sup>20</sup>, local entropy<sup>21</sup>, etc.) have been proposed and tested to find a relation between structure and dynamics. Some of the proposed descriptors correlate well with future dynamics to some extent. While some descriptors have apparently different physical origins, some others share essentially the same information. Besides, recent studies have highlighted the role of coarse-graining on descriptors<sup>19,22</sup>. Coarse-graining with an optimal length scale enhances the correlation between dynamics and

averaged descriptors. These extensive investigations lead to a question: What are the relevant (or main) physical descriptors responsible for glassy dynamics among a pool of many proposed descriptors with various coarse-graining lengths? This question is also related to identifying the correct theory of glass transition among numerous proposals based on distinct physical mechanisms<sup>23</sup>. Typically, each theoretical scenario attributes a specific quantity to glassy dynamics, such as entropy<sup>24,25</sup> in the random first-order transition theory<sup>26,27</sup>, elasticity<sup>20</sup> in the solid-that-flow picture<sup>28</sup>, and so on. Thus, selecting the main physical descriptor(s) relevant to glassy dynamics would shed light on this fundamental theoretical issue.

Recently, machine learning techniques have been applied to reveal the structure-dynamics relationship from another angle. Machine learning techniques specific to this problem consist of two approaches, supervised and unsupervised learning. Supervised learning in this context employs both structure (configurations, structural descriptors) and dynamics (propensity, mobility) data in training and forecasts the future dynamics of each particle of an unseen configuration<sup>29-34</sup>. Moreover, very recent studies have identified important descriptors for aged<sup>35,36</sup> and active<sup>37</sup> glasses using post-hoc explainable machine learning methods. Unsupervised learning, instead, does not require dynamics data in training<sup>38-40</sup>. It consists of dimensional reduction of structural data and their clustering. A clustered map based solely on structural information is compared with a dynamical propensity map, and their correlations are evaluated<sup>40</sup>.

In the language of machine learning and data science, identifying relevant features (descriptors in our case) is referred to as feature selection. One can easily perceive that the feature selection process is reduced to comparing two features, say $X^{(1)}$ and $X^{(2)}$, and asking whether one is more informative than the other. By using the standard statistical measures such as the Pearson coefficient $\rho$, one can perform this comparison by computing $\rho$ between $X^{(1)}$ (or $X^{(2)}$) and $Y$, denoted as, $\rho_{X^{(1)}, Y}$ (or $\rho_{X^{(2)}, Y}$), where $Y$ is the desired output data (dynamic propensity in our case). Based on the val-
ues of $\rho_{X^{(1)}, Y}$ and $\rho_{X^{(2)}, Y}$, one can conclude which variable, $X^{(1)}$ or $X^{(2)}$, is more relevant to $Y$. However, when the output data $Y$ is not available (the case of unsupervised learning), it is hard to assess whether one is more informative than the other by using solely $X^{(1)}$ and $X^{(2)}$, because $\rho_{X^{(1)}, X^{(2)}}$ is symmetric with respect to exchanging $X^{(1)}$ and $X^{(2)}$ by construction, i.e., $\rho_{X^{(1)}, X^{(2)}} = \rho_{X^{(2)}, X^{(1)}}$. This is also the case for other statistical measures such as the mutual information $I$, i.e., $I(X^{(1)}, X^{(2)}) = I(X^{(2)}, X^{(1)})$.

Recently, Ref. 41 introduced a new information measure, called information imbalance, $\Delta(A \to B)$. It is defined in an asymmetric way to assess how much $A$ predicts (or determines) $B$. By comparing $\Delta(A \to B)$ and $\Delta(B \to A)$, one can ask if $A$ and $B$ are equivalent, independent, or if one is more informative than the other. In particular, this new method was applied to selection of relevant policies on the Covid-19 epidemic state relevant for the number of deaths in the future$^{41}$. Besides, it was applied to dimensional reduction of structural descriptors for atomic systems$^{41}$. This pioneering work is followed by Ref. 42, where they used the information balance techniques to detect causality in simulated dynamical systems and real-world data taken from the electrophysiological signals of human brains. In Refs. 43 and 44, the information imbalance analysis is used to assess the predictivity and interpretability of various descriptors for systems of water molecules. In addition to scientific applications, Ref. 45 conducted feature selection from clinical databases to obtain a maximally informative set of features relevant to the severity of the disease.

In this paper, we employ information imbalance techniques to select structural features relevant to glassy slow dynamics. In particular, we select several relevant features from a pool of physically motivated descriptors composed of local energy, local number density, local volume fraction, coordination number$^{46}$, bond-orientational order$^{47}$, and steric bond order$^{19}$, with various coarse-grained lengthscales. We perform feature selection in both supervised (using dynamic propensity) and unsupervised (without dynamical data) settings.

First, we analyze the correlation between dynamic propensity and various physical descriptors varying the coarse-grained length scale using the conventional Pearson coefficient. We then confirm that the information imbalance reproduces the general trend obtained by the Pearson coefficient. Next we select relevant features in a supervised way by minimizing the information imbalance between descriptors and dynamic propensity. We also perform feature selection in an unsupervised way using only structural data. Finally, we use the selected features (from supervised and unsupervised manners) as input for the prediction of dynamics by machine learning and compare performance among different feature selection methods.

This paper is organized as follows. Section II describes our simulation model, physically motivated descriptors, and a concise introduction to the information imbalance techniques. We then present in Sec. III the results for the comparison between the Pearson coefficient and the information imbalance, feature selections by both supervised and unsupervised

learning, and their application to machine learning dynamic propensity. We finally discuss future perspectives in Sec. IV.

## II. METHODS

### A. Simulation model

We study a three-component glass-forming liquid model composed of small (S), medium (M), and large (L) particles in two-dimensions. We add an additional species (M) to a binary model (composed of S and L species) studied widely$^{48-50}$.

Interaction between two particles is given by the Lennard-Jones potential, $v_{\alpha\beta}(r) = 4\epsilon_{\alpha\beta}[(\sigma_{\alpha\beta}/r)^{12} - (\sigma_{\alpha\beta}/r)^{6}]$, where $\alpha, \beta = S, M, L$. The potentials are modified to be twice continuously differentiable functions at a cutoff, following Ref. 49. The parameters, $\sigma_{\alpha\beta}$ and $\epsilon_{\alpha\beta}$, for the potentials are given as follows: $\sigma_{LL} = 2\sin(\pi/5) \approx 1.18$, $\sigma_{SS} = 2\sin(\pi/10) \approx 0.62$, $\sigma_{LS} = 1$, $\sigma_{LM} = (\sigma_{LL} + \sigma_{LS})/2$, $\sigma_{MS} = (\sigma_{LL} + \sigma_{SS})/2$, and $\sigma_{MM} = (\sigma_{LL} + \sigma_{SS})/2$. $\epsilon_{LL} = 1/2$, $\epsilon_{SS} = 1/2$, $\epsilon_{LS} = 1$, $\epsilon_{LM} = (\epsilon_{LL} + \epsilon_{LS})/2$, $\epsilon_{MS} = (\epsilon_{LS} + \epsilon_{SS})/2$, and $\epsilon_{MM} = (\epsilon_{LL} + \epsilon_{SS})/2$. All particles have the same mass, which is set to 1. The number of particles is $N = N_S + N_M + N_L = 4000$, where $N_S = 1760$, $N_M = 800$, and $N_L = 1440$ are the number of particles for small (S), medium (M), and large (L) particles, respectively. We consider NVT canonical ensemble. The number density is $\rho = N/L^2 = 1.024$, where $L$ is the linear box length.

We perform Monte-Carlo (MC) dynamics simulations composed of translational displacements$^{51}$. For the MC moves, a particle is randomly chosen and displaced by a vector randomly drawn within a square box of linear size $\delta_{max} = 0.12$. The move is accepted according to the Metropolis acceptance rule, satisfying the detailed balance condition. Although MC simulations do not have an actual physical timescale, one can consider a time $t$ in units of MC sweeps, comprising $N$ MC trials. Indeed, it has been reported that Monte Carlo dynamics behaves essentially the same as other type of physical dynamics (e.g., Newtonian and Brownian dynamics) in terms of glassy slow dynamics$^{52}$. Thus, we analyze Monte Carlo dynamics following trajectories of particles and computing various time-dependent observables (see below).

To study dynamics of Monte-Carlo simulations, we measure the self-intermediate scattering function$^{53}$. It is known that the dynamics of two-dimensional glass-forming liquids is strongly affected by long-wavelength fluctuations whose magnitude increases with the system size$^{54-57}$. To remove this effect and study the intrinsic features of glassy dynamics, by convention, we employ the cage-relative version of the self-scattering function, which is defined by

$$F_s^{\text{CR}}(k,t) = \frac{1}{N} \sum_{i=1}^{N} \left\langle e^{-i\mathbf{k} \cdot \Delta \mathbf{r}_i^{\text{CR}}(t)} \right\rangle_{\text{time}}, \quad (1)$$

where $\langle \cdots \rangle_{\text{time}}$ denotes the time average. $\Delta \mathbf{r}_i^{\text{CR}}(t)$ is the cage-relative displacement$^{54}$ given by

$$\Delta \mathbf{r}_i^{\text{CR}}(t) = \Delta \mathbf{r}_i(t) - \frac{1}{n_i} \sum_{j \in NN} \Delta \mathbf{r}_j(t), \quad (2)$$
<figure><img src="image_2.png" /><figcaption>FIG. 1. Self-intermediate scattering function obtained by the Monte-Carlo dynamics.</figcaption></figure>

where $\Delta\mathbf{r}_i(t) = \mathbf{r}_i(t) - \mathbf{r}_i(0)$ is the displacement vector of the position $\mathbf{r}_i$ of the $i$-th particle and $n_i$ is the number of neighbor particles. The neighbors $NN$ are defined by the particles located within a circular cutoff whose radius is $1.4\sigma_{\alpha\beta}$. We note that the center of mass drifts at a longer time in the Monte-Carlo dynamics, and we subtracted such effect in our analysis. We set $k = 6.3$, which is close to the first peak of the static structure factor.

We show $F_s^{CR}(k, t)$ in Fig. 1. $F_s^{CR}(k, t)$ relaxes slowly with decreasing temperature. In particular, it shows a two-step relaxation at lower temperatures, which is a hallmark of glassy slow dynamics observed in various other types of dynamics, such as Newtonian and Brownian dynamics$^{52}$. We then define a relaxation time $\tau_\alpha$ by $F_s^{CR}(k, \tau_\alpha) = 1/e$. $\tau_\alpha$ grows steeply, reflecting glassy slow dynamics, as expected. In this study, we focus on $T = 0.30$, which is the lowest temperature we can reach within our computational timescale. The structural relaxation time $\tau_\alpha$ of this temperature is $\tau_\alpha \simeq 4 \times 10^6$.

To see heterogeneous dynamics in real space and assess the influence of static structure, we use the iso-configurational ensemble$^{6,8}$. We prepare an equilibrium configuration sampled at temperature $T = 0.30$. We then perform the MC dynamics at $T = 0.30$ many times with different seeds for random numbers, starting from the same initial equilibrium configuration. In particular, we compute a propensity parameter, $p_i$, which is given by

$$p_i(t) = \langle |\Delta\mathbf{r}_i^{CR}(t)| \rangle_{\text{iso}}, \quad (3)$$

where $\langle \cdots \rangle_{\text{iso}}$ is average over 30 different trajectories starting from a given static configuration. In this paper, we mainly focus on the propensity of heterogeneous dynamics at the structural relaxation timescale, namely, $p_i(\tau_\alpha)$. In Fig. 2(a), we present $p_i(\tau_\alpha)$ showing dynamical heterogeneity. Regions with larger and smaller values of $p_i$ correspond to mobile and immobile regions, respectively.

We also compute an another propensity parameter, $q_i(t)$,

<figure><img src="image_3.png" /><figcaption>FIG. 2. (a, b): Map of dynamic propensities, $p_i(\tau_\alpha)$ measured by the displacement vectors (a) and $q_i(\tau_\alpha)$ measured by the bond-breaking (b). (c): Map of the steric bond order parameter, $\Theta_i$, with the coarse-grained length $\ell = 2.0$. (d): Map of the predicted propensity, $\hat{p}_i$, obtained by linear regression with $d = 15$ features selected by the supervised information imbalance method in Sec. III B.</figcaption></figure>

based on the bond-breaking approach$^{58,59}$. $q_i(t)$ is defined by

$$q_i(t) = \langle n_i(t) / n_i(0) \rangle_{\text{iso}}, \quad (4)$$

where $n_i(0)$ is the number of neighbor particles within a cutoff radius $1.4\sigma_{\alpha\beta}$ of the particle $i$ and $n_i(t)$ is the number of particles that were initially part of the neighbors (at $t = 0$) and are still inside a cutoff $1.8\sigma_{\alpha\beta}$ at time $t$. We present $q_i(\tau_\alpha)$ in Fig. 2(b). $p_i(\tau_\alpha)$ and $q_i(\tau_\alpha)$ show very similar spatial patterns, suggesting that they share essentially the same information.

## B. Physically motivated descriptors

We consider various physically motivated descriptors characterizing the local structural environment.

The local potential energy $u_i$ for the particle $i$ is given by

$$u_i = \frac{1}{2} \sum_{j \neq i} v_{\alpha_i \beta_j}(r_{ij}), \quad (5)$$

where $v_{\alpha_i \beta_j}(r_{ij})$ is the pair-wise Lennard-Jones potential which has a cutoff at $2.5\sigma_{\alpha_i \beta_j}$.

The coordination number $z_i$ for the $i$-th particle is defined by the number of neighbor particles within $r_{ij} < 1.5\sigma_{\alpha_i \beta_j}$ (corresponding to near the first minimum of the radial distribution function).

We compute the conventional bond orientational order parameter in two-dimensions, $\Psi_{6,i}$, for the $i$-th particle, defined by

$$\Psi_{6,i} = \frac{1}{z_i} \left| \sum_{j=1}^{z_i} e^{\sqrt{-1} 6 \theta_{ij}} \right|, \quad (6)$$

where $\theta_{ij}$ is the angle between $\mathbf{r}_{ij} = \mathbf{r}_j - \mathbf{r}_i$ and the x-axis. The condition for the nearest neighbors is again set to $r_{ij} < 0$.
1.5σ<sub>αᵢβⱼ</sub>. Ψ₆ᵢ is often used to characterize hexagonal order in glassy and jamming systems<sup>47,60</sup>. Ψ₆ᵢ takes 1 for perfectly hexagonal packings, where as Ψ₆ᵢ takes lower values for disordered packings.

We also compute a steric bond order parameter Θᵢ introduced by Tong and Tanaka¹⁹. Here we consider a central particle *i* and its neighbor particles. The definition of the neighbors is the same as Ψ₆ᵢ. For each pair ⟨jk⟩ of neighbor particles next to each other, we measure the angle between **r**ᵢⱼ and **r**ᵢₗ, denoted as θⱼₖ. The reference configuration with these three particles, *i*, *j*, and *k*, perfectly just in touch, produces the angle indicated as θᵢⱼ<sup>ref</sup>. Practically θᵢⱼ<sup>ref</sup> is computed by σ<sub>αᵢβⱼ</sub>, σ<sub>αᵢβₖ</sub>, and σ<sub>αᵢβₖ</sub>, using the cosine formula. Then we define the order parameter for the *i*th particle as

$$ \Theta_i = \frac{1}{z_i} \sum_{\langle jk \rangle} \left| \theta_{jk} - \theta_{jk}^{\text{ref}} \right|, \qquad (7) $$

where ⟨jk⟩ denotes the summation over all pairs of neighbors. If particles form sterically favored, well-packed configurations, Θᵢ produces smaller value, because θⱼₖ would be close to the reference, θᵢⱼ<sup>ref</sup>. Instead, disordered packings generally take larger Θᵢ, since θⱼₖ would strongly deviate from θᵢⱼ<sup>ref</sup>. Thus, Θᵢ characterizes amount of disorder, which would play the similar roles as 1 − Ψ₆ᵢ, yet Θᵢ is more sensitive order parameter for multi-components (or polydisperse) systems¹⁹.

We next consider the coarse-grained local number density,

$$ \bar{\rho}_i(\ell) = \sum_{j \in \mathcal{N}_i} e^{-r_{ij}/\ell}, \qquad (8) $$

where $\mathcal{N}_i$ is the set of particles of *i* including *i* itself and *ℓ* is the coarse-graining lengthscale. These particles are inside the circle with the radius $R_N$ centered at the position of *i*. We set $R_N = 20\sigma_{LS}$ in this study.

We also compute local volume (area) fraction,

$$ \bar{\phi}_i(\ell) = \sum_{j \in \mathcal{N}_i} (\sigma_{\alpha_i \beta_j})^2 e^{-r_{ij}/\ell}. \qquad (9) $$

We perform the above coarse-graining procedure for all the other descriptors, $x_i = u_i, z_i, \Psi_{6,i}$, and $\Theta_i$, as given by

$$ \bar{x}_i(\ell) = (\bar{\rho}_i)^{-1} \sum_{j \in \mathcal{N}_i} x_j e^{-r_{ij}/\ell}, \qquad (10) $$

which provides us with $\bar{u}_i(\ell)$, $\bar{z}_i(\ell)$, $\bar{\Psi}_{6,i}(\ell)$, and $\bar{\Theta}_i(\ell)$. We vary the coarse-graining lengthscale $\ell$ from $\ell = 0.5$ to $\ell = 5.0$ for each descriptor. Thus in total, we have $M = 60$ features (6 different physically motivated descriptors times 10 different coarse-graining lengthscales).

For each feature, we normalize it to have zero mean and unit variance. After normalization, all features constitute a feature vector $X_{\text{Full},i}$ for a particle *i*, given by

$$ X_{\text{Full},i} = (X_i^{(1)}, X_i^{(2)}, \dots, X_i^{(M)}). \qquad (11) $$

This is the structural input for the information imbalance analysis described below.

### C. Information imbalance

We recap the information imbalance introduced in Ref. 41. Suppose that we have two features (physically motivated descriptors in our case), A and B, characterizing a data point. We wish to ask if A and B are equivalent, independent, or if one is more informative than the other. The last question, i.e., A is more informative than B, is particularly challenging for standard statistical analysis where A and B are treated equally or symmetrically. Recall, for example, that the Pearson coefficient $\rho_{A, B}$ (see Sec. III A for the definition) and mutual information $I(A, B)$ are defined in a symmetric way, $\rho_{A, B} = \rho_{B, A}$ and $I(A, B) = I(B, A)$, respectively. Hence, although $\rho_{A, B}$ and $I(A, B)$ can quantify the correlation between A and B, they alone cannot tell, say, A is more informative than B.

The information imbalance analysis between A and B is constructed in an asymmetric way, which allows us to tell if one is more informative than the other. We first consider a feature space spanned by A and B as schematically shown in Fig. 3(a). Each of A and B can generally be vectors, but for simplicity, we consider both to be scalars, forming a two-dimensional plane. Here we have four data points labeled by *i*, *j*, *k*, and *l*, respectively. Let us focus on the data point *i* and make a ranking based on distance along the A-axis and the B-axis separately. We define the Euclidean distance by $d_{ij}^A = |A_i - A_j|$ and $d_{ij}^B = |B_i - B_j|$, where $A_i$ and $B_i$ are the A and B component of the data point *i*. Note that one can also use non-Euclidean distance in general. For example, in Fig. 3(a), data point *j* is the closest data point of *i* in terms of A, which is ranked first and denoted as $r_{ij}^A = 1$. In general, $r_{ij}^A$ is the rank of *j* relative to *i*, obtained by sorting the pairwise distances between *i* and rest of the points from smallest to largest. The data points *k* and *l* are the second and third closest data points in terms of A, which are written as $r_{ik}^A = 2$ and $r_{il}^A = 3$, respectively. Similarly, one can construct ranking based solely on B, which are written as $r_{ij}^B = 3$, $r_{ik}^B = 1$, and $r_{il}^B = 2$.

We now imagine a situation where A contains information about B (to some extent). In such a case, the organization of N data points in the feature space would be as follows. When two data points, say *i* and *j*, are close along A, they would also be close along B. This idea can be quantified by monitoring an average rank of $r_{ij}^B$ when $r_{ij}^A = 1$, i.e., $\frac{1}{N}\sum_{i=1}^N r_{ij}^B$ given $r_{ij}^A = 1$. If $r_{ij}^B$ is close to 1 for most of *i*, we would get $\frac{1}{N}\sum_{i=1}^N r_{ij}^B \simeq 1$ on average. Instead, when A contains no information about B, $r_{ij}^B$ (given $r_{ij}^A = 1$) would be completely random, and hence we get $\frac{1}{N}\sum_{i=1}^N r_{ij}^B \simeq N/2$. We will use $N/2$ for a normalization factor below.

Finally, based on the above consideration, we define the information imbalance $\Delta(A \to B)$ using the average rank of
<figure><img src="image_4.png" /><figcaption>FIG. 3. (a): Schematic illustration of data points in the feature space composed of feature A and B. (b): Information imbalance plane. Four distinct relationships between A and B are presented schematically (see the main text).</figcaption></figure>

$r_{ij}^B$ given $r_{ij}^A = 1$ with the normalization factor:

$$ \Delta(A \to B) = \left(\frac{N}{2}\right)^{-1} \frac{1}{N} \sum_{i=1}^{N} \sum_{j=1}^{N} r_{ij}^{B} \\ = \frac{2}{N^2} \sum_{i=1}^{N} \sum_{j=1}^{N} \delta_{1,r_{ij}^A} r_{ij}^B, \qquad (12) $$

where $\delta_{1,r_{ij}^A}$ is the Kronecker's delta selecting $r_{ij}^A = 1$. We also used $r_{ii}^A = 0$ and $r_{ii}^B = 0$. Because of the normalization factor, $N/2$, we get $\Delta(A \to B) \simeq 0$ when A predicts the ranking of B strongly, namely, A contains information about B, whereas $\Delta(A \to B) \simeq 1$ when A cannot determine the ranking of B at all, or A is not informative about B.

As one expect from the definition in Eq. (12), one can also define $\Delta(B \to A)$ by computing the average rank of $r_{ij}^B$ given $r_{ij}^A = 1$. By measuring both $\Delta(A \to B)$ and $\Delta(B \to A)$ together, we can finally ask if A and B are equivalent, independent, or if one is more informative than the other. The limiting cases relevant to our study is summarized as follows$^{41}$.

* A and B are equivalent, carrying identical information when $\Delta(A \to B) \simeq 0$ and $\Delta(B \to A) \simeq 0$.

* A and B are independent (or orthogonal), carrying completely independent information when $\Delta(A \to B) \simeq 1$ and $\Delta(B \to A) \simeq 1$.

* A is more informative than B when $\Delta(A \to B) \simeq 0$ and $\Delta(B \to A) \simeq 1$.

* A and B share information symmetrically, carrying both identical and independent information when $0 \ll \Delta(A \to B) = \Delta(B \to A) \ll 1$.

In practice, it is convenient to plot a two-dimensional plane of $\Delta(A \to B)$ and $\Delta(B \to A)$, called the information imbalance plane, as shown in Fig. 3(b).

Up to here, we considered information imbalance between A and B. In applications, one often studies information imbalance between a feature (vector) composed of A and B, denoted as $(A, B)$, and a (scalar) feature C. In such a case, one has to define the distances properly. For example, one can define distances for the features, $(A, B)$ and C, by $d_{ij}^{(A,B)} = \sqrt{(A_i - A_j)^2 + (B_i - B_j)^2}$ and $d_{ij}^C = |C_i - C_j|$, respectively, using the Euclidean distance. Then the ranking $r_{ij}^{(A,B)}$ and $r_{ij}^C$ can be obtained based on $d_{ij}^{(A,B)}$ and $d_{ij}^C$, respectively. The above setting allows us to compute $\Delta((A, B) \to C)$ and $\Delta(C \to (A, B))$ appropriately.

In this paper, we compute the information imbalance, $\Delta(A \to B)$ and $\Delta(B \to A)$, as well as the Pearson coefficient (see below), for each configuration composed of $N = 4000$ particles. We average them over 10 different configurations.

### III. RESULTS

#### A. Pearson coefficient and information imbalance

We first analyze the relationship between dynamic propensity and physically motivated descriptors by using the conventional Pearson coefficient $\rho_{A, B}$, as a reference for the information imbalance analysis below. $\rho_{A, B}$ quantifies the correlation between two variables A and B:

$$ \rho_{A, B} = \frac{1}{N} \sum_{i=1}^{N} \frac{(A_i - \mu_A)(B_i - \mu_B)}{\sigma_A \sigma_B}, \qquad (13) $$

where $\mu_A$ and $\mu_B$ are the mean of $A_i$ and $B_i$, respectively, and $\sigma_A$ and $\sigma_B$ are the standard deviation of $A_i$ and $B_i$, respectively.
Figure 4(a) shows the Pearson coefficient, $ρ_{X, p}$, between the propensity $p_i(τ_α)$ and all descriptors we considered (here simply denoted as $X$) with varying the coarse-graining length-scale $ℓ$. Note that because the volume fraction and coordination number have anti-correlation with propensity, we present the absolute value of the Pearson coefficient. A larger value of $|ρ_{X, p}|$ corresponds to a higher correlation (or anti-correlation) between the propensity and descriptor $X$.

We find that $\bar{\Theta}$ is the best descriptor in the plot, showing a non-monotonic $ℓ$ dependence. $ρ_{\bar{\Theta}, p}$ has the maximum value (around 0.56) at an optimal lengthscale near $ℓ = 2$, consistent with previous literature¹⁹. We show in Fig. 2(c) a map for $\bar{\Theta}_i(ℓ = 2)$. Indeed, it shows a similar pattern as the dynamical heterogeneity presented in Fig. 2(a). The other descriptors demonstrate a smaller correlation with dynamics either monotonically ($\bar{\rho}$, $\bar{\Psi}$) or non-monotonically ($\bar{u}$, $\bar{\phi}$, $\bar{z}$) dependence on $ℓ$.

Next, we analyze the same dataset by the information imbalance and ask if we get the same trend as the Pearson coefficient. Recall that when $Δ(X → p)$ is small, a feature $X$ predicts the propensity $p$ well. Hence, to compare with the Pearson coefficient, we plot the inverse of $Δ(X → p)$ in Fig. 4(c). We find that $1/Δ(X → p)$ follows a quite similar trend as the Pearson coefficient in terms of the order of the descriptors as well as their $ℓ$ dependence.

We also measure $Δ(p → X)$ and plot its inverse in Figs. 4(d). Interestingly, $1/Δ(p → X)$ also shows essentially the same behavior with a similar value, suggesting $Δ(X → p) \simeq Δ(p → X)$. Thus, in Fig. 4(b), we present the $Δ(X → p)$ versus $Δ(p → X)$ plot (information imbalance plane) for all features studied. Indeed, we find that all features concentrate on the straight line having a symmetrical form, $Δ(X → p) \simeq Δ(p → X)$. In particular, $\bar{\Psi}_6$ (for all $ℓ$) takes $Δ(\bar{\Psi}_6 → p) \simeq Δ(p → \bar{\Psi}_6) \simeq 1$, meaning that $\bar{\Psi}_6$ and propensity $p$ are independent (or orthogonal), sharing no information, at least in the simulation model we use. Instead, $\bar{\Theta}$ shows $Δ(\bar{\Theta} → p) \simeq Δ(p → \bar{\Theta}) \simeq 0.85$, meaning that $\bar{\Theta}$ and propensity $p$ share information to some extent. These trends are consistent with findings in the Pearson coefficient.

To better understand the information imbalance plane, we also compute the information imbalance between $p$ and the propensity $q$ measured by the bond-breaking in Eq. (4). We obtain $Δ(q → p) \simeq 0.21$ and $Δ(p → q) \simeq 0.26$, meaning that $p$ and $q$ share very similar information, as expected from the visual comparison between Figs. 2(a) and (b).

To conclude, the information imbalance analysis can characterize the relationship between structural features and dynamic propensity, reproducing the trend obtained by the Pearson coefficient. We also found that the information imbalance between each feature and propensity is symmetric. This means that although some descriptors such as $\bar{\Theta}$ and $\bar{u}$ have a higher correlation with propensity, no single feature is more informative than the propensity. This is the first non-trivial statement achieved by the information imbalance analysis.

## B. Feature selection by supervised learning

We next perform feature selection among a pool of $M = 60$ physically motivated descriptors by information imbalance in a supervised learning manner, using propensity data. This method is inspired by feature selection of relevant policies on the Covid-19 epidemic state in terms of the number of deaths in the future⁴¹. Recall that the full set of $M$ features is denoted as the vector, $X_{\text{Full}}$, which is given by

$$X_{\text{Full}} = (X^{(1)}, X^{(2)}, \dots, X^{(M)}). \quad (14)$$

Here, we drop off the particle index $i$ for simplicity. The goal is to select $d \ll M$ features out of $M = 60$ features and obtain a vector of selected features,

$$X_d = (X^{\text{1st}}, X^{\text{2nd}}, \dots, X^{\text{dth}}). \quad (15)$$

This operation corresponds to a dimensional reduction from $M$ to $d$.

The algorithm to obtain $X_d$ using the information imbalance is as follows. First, among $M(= 60)$ descriptors, we search $X^{\text{1st}}$ at which the information imbalance between the feature and propensity is minimized at $Δ(X_1 → p)$, where $X_1 = (X^{\text{1st}})$. We next search $X^{\text{2nd}}$ among the rest of $M - 1$ features such that the information imbalance is minimized at $Δ(X_2 → p)$, where $X_2 = (X^{\text{1st}}, X^{\text{2nd}})$. We repeat the above procedure until we find $X^{\text{dth}}$ among the rest of $M - (d - 1)$ features such that the information imbalance is minimized at $Δ(X_d → p)$. For each minimization, we perform a greedy algorithm to find the optimal feature. Since $X_d$ is found to have minimum at $Δ(X_d → p)$, at least within our algorithm, $X_d$ should have a high predictivity in terms of the propensity $p$. Note that in the above iterative procedure, we keep the previously selected features and find a new one among the rest of the features. Although it is not guaranteed to reach the global minimum of $Δ(X_d → p)$, this procedure avoids searching all combinatorially possible sets of $d$ features, which requires a huge computational cost⁴¹.

Figure 5(a) shows $Δ(X_d → p)$ as a function of $d$. $Δ(X_d → p)$ decays quickly as $d$ is increased and it saturates, say, above $d = 10$. This means that $X_{d=10}$ is already as informative as the full set, $X_{d=M} = X_{\text{Full}}$, in terms of prediction of the propensity. The selected features up to $d = 10$ from the 1st to $d$th, for the configuration presented in Fig. 2 for instance, are as follows: $\bar{\Theta}_i(ℓ = 1.5), \bar{\rho}_i(ℓ = 1), \bar{\rho}_i(ℓ = 4), \bar{\Theta}_i(ℓ = 5), \bar{u}_i(ℓ = 5), \bar{\phi}_i(ℓ = 4.5), \bar{z}_i(ℓ = 1.5), \bar{\Psi}_{6,i}(ℓ = 4), \bar{\phi}_i(ℓ = 0.5),$ and $\bar{z}_i(ℓ = 5)$. This lineup fluctuates slightly depending on configurations, yet the overall trend does not change.

One might assume that adding features, even irrelevant ones, would decrease (or at least not increase) $Δ(A → B)$. However, this is not always the case, as demonstrated in the following exercise⁴¹. In $Δ(X_d → p)$, we add an independent feature generated from Gaussian noise, which is entirely irrelevant. As shown in Fig. 5(a), $Δ(X_d → p)$ actually increases with the addition of this irrelevant feature.

Instead, as seen in Fig. 5(a), $Δ(p → X_d)$ is essentially flat as a function of $d$, which means that there is no improvement of predictivity from $p$ to $X_d$, as expected. This huge asymmetry
<figure><img src="image_5.png" /><figcaption>FIG. 4. (a): The absolute value of the Pearson coefficient between the propensity $p_i$ and various physical descriptors as a function of the coarse-grained length $\ell$. (b): Information imbalance plane composed of $\Delta(X \to p)$ and $\Delta(p \to X)$. Each physically motivated descriptor has 10 data points associated with different coarse-graining lengths (shown in the same color). (c, d): The inverse of the information imbalance $\Delta(X \to p)$ (c) and $\Delta(p \to X)$ (d).</figcaption></figure>

can be seen in the information imbalance plane in Fig. 5(b), where we show the evolution of ($\Delta(X_d \to p)$, $\Delta(p \to X_d)$) with increasing $d$. The first data point at $d=1$ associated with a single feature $X_1 = (X^{1\text{st}})$ locates on the straight line, showing symmetry, $\Delta(X_1 \to p) \simeq \Delta(p \to X_1)$, as we have already seen in Sec. III A. Yet, with increasing the number of selected features, ($\Delta(X_d \to p)$, $\Delta(p \to X_d)$) deviates from the straight line with a strong asymmetry, $\Delta(X_d \to p) \ll \Delta(p \to X_d)$. This means that $X_d$ is getting more and more informative than propensity $p$ with increasing $d$, while $p$ alone is getting harder to predict $X_d$. From this observation, we conclude that *multiple* features are required to construct a set of selected features that is more informative than the dynamic propensity.

The above trend is analogous to the feature selection of relevant policies on the Covid-19 epidemic studied in Ref. 41. The authors of Ref. 41 selected relevant policies for the number of deaths in the future among various conducted policies (e.g., school closing, stay-home restrictions, contact tracing, etc.). They found that no *single* policy (feature) is sufficient on its own, yet the combination of *multiple* policies is getting

effective in terms of predicting the number of future deaths.

### C. Feature selection by unsupervised learning

Feature selection done in Sec. III B used the propensity data in a supervised learning manner. Here, we wish to perform feature selection without using propensity (or any dynamical data), which corresponds to an unsupervised learning method.

Again, the goal is to select features, constructing a feature vector, $X_d$ in Eq. (15), yet without using the propensity. Ideally, the selected features $X_d$ and the original set $X_{\text{Full}}$ should be nearly equivalent, carrying essentially the same information, while $d \ll M$. By using the language of information imbalance, this idea corresponds to $\Delta(X_{\text{Full}} \to X_d) \simeq \Delta(X_d \to X_{\text{Full}}) \simeq 0$. To achieve this condition we wish to minimize both $\Delta(X_{\text{Full}} \to X_d)$ and $\Delta(X_d \to X_{\text{Full}})$ simultaneously. Thus, we define the symmetrized information imbalance$^{41}$, $\bar{\Delta}$, be-
<figure><img src="image_6.png" /><figcaption>FIG. 5. Feature selection by the supervised approach. (a): Information imbalance, $\Delta(X_d \to p)$ and $\Delta(p \to X_d)$, as a function of the number of selected features, $d$. The errorbars correspond to the standard deviation from 10 different configurations. The last triangular data point at $d = 16$ represents the addition of one feature generated from Gaussian noise. (b): The same data are presented in the information imbalance plane. The number inside a point represents $d$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIG. 6. Feature selection by the unsupervised approach. (a): Decay of the symmetric information imbalance, $\bar{\Delta}(X_{Full}, X_d)$, as a function of $d$, together with $\Delta(X_{Full} \to X_d)$ and $\Delta(X_d \to X_{Full})$. The errorbars correspond to the standard deviation from 10 different configurations. (b): The same data for $\bar{\Delta}(X_{Full} \to X_d)$ and $\Delta(X_d \to X_{Full})$ are presented in the information imbalance plane. The number inside a point represents $d$.</figcaption></figure>

tween $X_{Full}$ and $X_d$ as

$$ \bar{\Delta}(X_{Full}, X_d) = \frac{\Delta(X_{Full} \to X_d) + \Delta(X_d \to X_{Full})}{\sqrt{2}}. \quad (16) $$

We then minimize $\bar{\Delta}(X_{Full}, X_d)$ using a greedy algorithm similar to the supervised learning case in Sec. III B.

First, among $M(= 60)$ features, we search $X^{1st}$ at which the symmetrized information imbalance between $X_{Full}$ and $X_1$ is minimized at $\bar{\Delta}(X_{Full}, X_1)$, where $X_1 = (X^{1st})$. We next search $X^{2nd}$ among the rest of $M - 1$ features such

that the symmetrized information imbalance is minimized at $\bar{\Delta}(X_{Full}, X_2)$, where $X_2 = (X^{1st}, X^{2nd})$. We repeat the above procedure until we find $X^{dth}$ among the rest of $M - (d - 1)$ features such that the symmetrized information imbalance is minimized at $\bar{\Delta}(X_{Full}, X_d)$.

In Fig. 6(a), we show the decay of $\bar{\Delta}(X_{Full}, X_d)$ as a function of $d$, together with $\Delta(X_{Full} \to X_d)$ and $\Delta(X_d \to X_{Full})$. We confirm that both of the terms, $\Delta(X_{Full} \to X_d)$ and $\Delta(X_d \to X_{Full})$, are minimized well (notice that the Y-axis is a log scale), yet decaying with a different speed, and
show a saturation. The corresponding information imbalance plane is shown in Fig. 6(b). First, we find a strong asymmetry at $d = 1$, $\Delta(X_{\text{Full}} \to X_1) \ll \Delta(X_1 \to X_{\text{Full}})$, meaning that $X_{\text{Full}}$ is much more informative than $X_1$, as expected. Yet, with increasing $d$, symmetry is getting restored, converging to the origin, (0,0). The selected features up to $d = 10$ from the 1st to $d$th, for the configuration presented in Fig. 2 are as follows: $\Theta_i(\ell = 4.5)$, $\bar{z}_i(\ell = 3.5)$, $\bar{\Psi}_{6,i}(\ell = 5)$, $\bar{\rho}_i(\ell = 4.5)$, $\bar{\phi}_i(\ell = 0.5)$, $\bar{u}_i(\ell = 4.5)$, $\bar{\phi}_i(\ell = 3.5)$, $\bar{u}_i(\ell = 1.5)$, $\Theta_i(\ell = 1)$, and $\bar{\rho}_i(\ell = 1.5)$. This lineup is quite similar to the one obtained by the supervised learning method in Sec. III B. Yet the order to achieve each lineup is very different. This difference will be highlighted in the machine learning prediction (see below).

### D. Application to machine learning glassy dynamics

Having identified a reduced set of relevant structural features, $X_d$, from the original full set, $X_{\text{Full}}$, using information imbalance in both supervised and unsupervised manners, we now apply them to machine learning glassy dynamics to predict dynamic propensity.

The study of machine learning glassy dynamics was initiated by a series of papers by the group of Andrea Liu using a support vector machine that classifies mobile and immobile particles by making a hyperplane in a high-dimensional space of structural descriptors$^{29,30}$. This line of research was followed by neural network techniques, including cutting-edge graph neural networks, in order to increase the accuracy of prediction$^{31-34}$. Typically, neural network methods outperform other methods, yet their necessity is still highly debated. In particular, linear regression with elaborated structural descriptors provides high accuracy as good as a deep learning method$^{61}$.

Thus, in this study, we employ simple linear regression techniques. To predict the dynamic propensity from the structural information, we construct a linear model,

$$ \hat{p}_i = w^T X_{d,i} + w^{(0)}, \quad (17) $$

where $w = (w^{(1)}, w^{(2)}, \dots, w^{(d)})$ is the weight vector and $w^{(0)}$ is the intercept. $w$ and $w^{(0)}$ are obtained by minimizing the mean squared error loss function, $\mathcal{L}^{\text{MSE}} = \frac{1}{2N} \sum_{i=1}^{N} (\hat{p}_i - p_i)^2$.

In Fig. 2(d), we show a map for predicted propensity $\hat{p}_i$, using the supervised scheme with $d = 15$, for an unseen test configuration, which is the same configuration as Fig. 2(a) where the ground truth $p_i$ are measured. Comparing with Fig. 2(a), we visually confirmed that the linear model in Eq. (17) can predict the shape of dynamical heterogeneity qualitatively well. We also confirmed that the unsupervised method with $d = 15$ performs prediction similarly (see below).

Besides, as a benchmark, we employ a supervised feature selection using a L1 norm penalty inspired by Lasso regression$^{62}$. In this method, we first minimize the loss function of Lasso regression, i.e., $\mathcal{L}^{\text{Lasso}} = \mathcal{L}^{\text{MSE}} + \frac{\alpha}{2} \sum_{f=1}^{M} |w^{(f)}|$, with a fixed strength of penalty $\alpha$. Because of the penalty term, some

weights $w^{(f)}$ become strictly zero$^{62}$. Thus, the remaining features associated with non-zero weights can be considered relevant selected features. The number of selected features decreases with increasing $\alpha$. We tune the value of $\alpha$ to select $d$ features. We call this method of feature selection “Lasso selection”. We then perform the standard linear regression (minimizing $\mathcal{L}^{\text{MSE}}$ alone) by using $d$ selected features from the Lasso selection.

We also perform a random feature selection serving as a benchmark for the unsupervised method. We randomly select $d$ features from the pool of $M$ features. For each configuration, we performed 1000 realizations and the subsequent linear regression. We then averaged the results over them.

Scatter plots comparing the ground truth $p_i$ and the prediction $\hat{p}_i$ are shown in Appendix A. We then quantify the accuracy of the prediction by using the Pearson coefficient, $\rho_{p,\hat{p}}$, between $p_i$ and $\hat{p}_i$ for all feature selection methods we performed. In Fig. 7, we show $\rho_{p,\hat{p}}$ as a function of the number of selected features, $d$. $\rho_{p,\hat{p}}$ is computed for all particles regardless of species. Species-wise values for $\rho_{p,\hat{p}}$ are reported in Appendix B.

Both of the supervised methods, the Lasso selection and information imbalance ones, quickly reach a high value, $\rho_{p,\hat{p}} \simeq 0.7$ after $d = 2$, and seem to saturate, say, above $d = 10$. The first two selected features obtained by the supervised information imbalance method is typically $\Theta_i(\ell)$ with $\ell \simeq 2-3$ and $\bar{\rho}_i(\ell)$ with $\ell \simeq 0.5-1.5$, which is consistent with high Pearson coefficient with propensity observed in Fig. 4(a). We wish to emphasize that the relative merit of the information imbalance method with respect to Lasso selection is that the former does not assume any model a priori, i.e., a non-parametric method, while the latter assumes a specific model (the linear model in Eq. (17)).

Turning our attention to the unsupervised methods, both the information imbalance and random feature selection methods converge slower than the supervised methods. In particular, low values of $d$ result in a significant error in predicting the dynamics compared with supervised methods, and the unsupervised information imbalance method converges to similar performance as the supervised methods only beyond, say, $d = 10$. This is consistent with the fact that $d = 10$ selected features by the supervised and unsupervised methods are quite similar. In contrast, random feature selection reaches a similar value only at the last data point, $d = 15$, or slightly beyond, suggesting that the information imbalance more effectively selects useful features in the range of $d \simeq 10-15$ by comparing among features.

As we have seen, while the accuracy of the supervised information imbalance method jumps at $d = 2$ directly to an approximate level of the large $d$ limit, the unsupervised one only gradually climbs up to that limit. Thus, the order of selection of each feature is quite different in the two methods, as we have confirmed. This reflects that while the unsupervised information imbalance selects the $d$ features based on representing the best cloud of points of the full feature space, the dynamic propensity is best encoded in a different set of $d$ features that the supervised method manages to capture. Thus, the information imbalance techniques provide us with insight
<figure><img src="image_8.png" /><figcaption>FIG. 7. Pearson coefficient between the ground truth propensity $p_i$ at the structural relaxation time $\tau_\alpha \simeq 4 \times 10^6$ and the predicted one $\hat{p}_i$ by the linear regression using the selected $d$ features. Results from four different methods (see the main text in detail) are compared. The errorbars correspond to the standard deviation from 10 different (unseen) configurations.</figcaption></figure>

into the organization of structural and dynamical data points in high-dimensional feature space.

We focused on predicting glassy dynamics on the timescale of structural $\alpha$ relaxation, using the propensity, $p_i(t)$, measured at $t = \tau_\alpha \simeq 4 \times 10^6$. Additionally, we made predictions for a shorter timescale, $t = 5 \times 10^4$, which is close to $\beta$ relaxation timescales, as described in Appendix C. Although the Pearson coefficient values are generally lower, consistent with previous studies$^{63}$, the overall trend remains similar to that observed for the $\alpha$ relaxation timescale.

## IV. DISCUSSION AND CONCLUSION

We applied the information imbalance techniques proposed recently$^{41}$ to select relevant structural features in glassy dynamics among a pool of various physically motivated descriptors with different coarse-graining lengthscales. Unlike the conventional statistical measures characterizing the relationship between variables A and B, such as the Pearson coefficient and mutual information, the information imbalance is constructed in an asymmetric manner, which allows us to ask not only if A and B are equivalent or independent, but also if one is more informative than the other, without referring to other variables.

We confirmed that the information imbalance between dynamic propensity and structural features can capture the general trend obtained by the Pearson coefficient, including its

dependence on the coarse-graining lengthscale. We then performed feature selection from the full set of structural features in a supervised manner using the propensity data. We found that although some descriptors (e.g., the steric order parameter) have higher correlations with propensity, no single descriptor is more informative than propensity, at least within our pool of descriptors. Instead, a set of selected multiple features is more informative than propensity, showing a strong asymmetry in the information imbalance plane. We also performed feature selection in an unsupervised manner without using any dynamics data. Remarkably, although the order of selection of each feature is quite different at smaller $d < 10$, the selected features at $d = 10$ are similar to those obtained by the supervised method using propensity. We confirmed that selected structural features in both supervised and unsupervised manners are indeed highly relevant and useful for predicting future glassy dynamics, comparing with other feature selection techniques.

This work corresponds to proof of the concept for the newly introduced idea of information imbalance. Thus, so far, we have focused only on one simulation model, one temperature (lowest temperature), and two timescales (structural $\alpha$ relaxation timescale in the main text and a $\beta$ relaxation timescale in Appendix C). It is interesting to explore and compare different models (2D vs. 3D$^{64}$, fragile vs. strong$^{65,66}$ etc.), different temperatures, and different timescales (including much lower temperatures and longer timescales$^{67}$). It is interesting to see if relevant descriptors and optimal coarse-graining lengthscale change depending on systems, the state point, and the timescale.

Besides, we studied only simple, basic physical descriptors, such as local energy and the steric order parameter. A wide variety of physical descriptors characterizing structure have been proposed, such as locally favored structures$^{10,15}$, vibrational modes$^{8,68}$, softness$^{30,69}$, caging potential$^{70,71}$, quasicalized excitations$^{72}$, mesoscopic elastic stiffness$^{20}$, local yield stress$^{49,50}$, overlap field$^{73,74}$, local configurational entropy$^{21}$, topological defects$^{75,76}$, and so on (see also Ref. 77 for a review). It would be insightful to include these descriptors in the information imbalance analysis to identify and select the most relevant features from the wide variety of proposed descriptors.

The above exercise would also be related to identifying the correct theory of glass transition. Many theories have been proposed based on quite distinct physical mechanisms, such as entropy$^{26,27}$, elasticity$^{28,78-80}$, locally preferred structures$^{9,14}$, etc. Information imbalance may help to identify the relevant physical mechanism associated with the selected feature. Furthermore, it would be worthwhile to investigate whether some of the descriptors mentioned above (and undiscovered ones) are more informative than dynamic propensity. Another scenario would be that none of the above features has a dominant role in dynamics. Instead, a combination of them becomes informative in terms of glassy dynamics. If this is the case, one cannot attribute a single dominant physical mechanism for glassy dynamics. Instead, multiple mechanisms would be hybridized in reality. Yet, there is also a possibility that structural data alone cannot contain complete
information about dynamics$^{81}$. The information imbalance techniques can be useful in assessing these different scenarios quantitatively.

## ACKNOWLEDGMENTS

We thank Ali Hassanali for introducing us to the information imbalance analysis. We also thank Gerhard Jung, Takumi Nagasawa, Rob Jack, and Daniele Coslovich for their discussions and suggestions. This work was supported by the Agence Nationale de la Recherche (ANR-19-P3IA-0003) under the Plan France 2030.

## DATA AVAILABILITY

All the source codes and dataset used in this paper are openly available at https://github.com/Anand-Sharma22/Feature-selection-Information-Imbalance.

## Appendix A: Scatter plots for propensity

In Fig. 8, we present scatter plots comparing the ground truth propensity $p_i$ with the predicted values $\hat{p}_i$, obtained using linear regression with $d = 15$ features selected by the supervised information imbalance method described in Sec. III B. These plots correspond to the configuration shown in Fig. 2 as an example. Figure 8(a) displays the scatter plot for all particles, while Figs. 8(b, c, d) show the scatter plots for small (b), medium (c), and large (d) particles, respectively. We note that in this paper, training (minimizing the loss function) is performed simultaneously for all particles for simplicity. Training each species separately could improve prediction performance$^{63}$.

The plot for smaller particles shows a wider range of propensity, as expected, while the data for medium and large particles span a narrower range. Overall, these plots exhibit a standard oval shape without skewness, indicating that the Pearson coefficient can appropriately quantify the degree of correlation.

## Appendix B: Species-wise plots

As anticipated from Fig. 8, the Pearson coefficient values differ across species. In Fig. 9, we present the species-wise Pearson coefficient, calculated separately for small, medium, and large particles. The value is highest for small particles, while medium and large particles show lower values, with the coefficient decreasing as particle size increases. Nevertheless, the overall trend remains consistent with the Pearson coefficient computed for all particles combined across species in Fig. 7.

<figure><img src="image_9.png" /><figcaption>FIG. 8. Scatter plots comparing the ground truth propensity $p_i$ with the predicted values $\hat{p}_i$ for all particles (a), small particles (b), medium particles (c), and large particles (d), obtained using linear regression with $d = 15$ features selected by the supervised information imbalance method.</figcaption></figure>

## Appendix C: Shorter timescale

In the main text, we studied the machine learning prediction of the propensity measured at the structural $\alpha$ relaxation time for $T = 0.30$. We also study a short timescale, $t = 5 \times 10^4$, close to $\beta$ relaxation timescales of the same temperature (see Fig. 1).

Figure 10 presents the Pearson coefficient obtained using the four feature selection methods we employed. We found that the Pearson coefficient at the shorter timescale is lower than at the structural relaxation time, consistent with previous reports$^{63}$. Additionally, for the supervised methods, the value at $d = 1$ is already close to the asymptotic value at large $d$. The feature identified at $d = 1$ by the supervised information imbalance method is $\Theta_i(\ell)$ with $\ell \simeq 1$, corresponding to a slightly smaller coarse-graining length scale than at the $\alpha$ relaxation timescale, which aligns with previous findings$^{19}$.

We also observed that the unsupervised information imbalance method underperforms compared to the random selection method around $d = 5$, but outperforms it above $d = 10$. This suggests that the features selected by the unsupervised information imbalance method (which rely solely on static information) near $d = 5$ are not relevant for short-time dynamics.

<footnote>$^{1}$L. Berthier and G. Biroli, “Theoretical perspective on the glass transition and amorphous materials,” Reviews of modern physics **83**, 587 (2011).</footnote>

<footnote>$^{2}$K. Binder and W. Kob, Glassy materials and disordered solids: An introduction to their statistical mechanics (World scientific, 2011).</footnote>

<footnote>$^{3}$M. D. Ediger, C. A. Angell, and S. R. Nagel, “Supercooled liquids and glasses,” The journal of physical chemistry **100**, 13200–13212 (1996).</footnote>

<footnote>$^{4}$L. Berthier, G. Biroli, J.-P. Bouchaud, L. Cipelletti, and W. van Saarloos, Dynamical heterogeneities in glasses, colloids, and granular media, Vol. 150 (OUP Oxford, 2011).</footnote>
<figure><img src="image_10.png" /><figcaption>FIG. 9. Pearson coefficient between the ground truth propensity $p_i$ at the structural relaxation time $\tau_\alpha \simeq 4 \times 10^6$ and the predicted propensity $\hat{p}_i$, obtained from linear regression using the selected $d$ features. The figure format follows that of Fig. 7, but the Pearson coefficients are computed separately for each species: small (a), medium (b), and large (c) particles, respectively.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIG. 10. Pearson coefficient between the ground truth propensity $p_i$ at a shorter time $t = 5 \times 10^4$ (near $\beta$ relaxation timescales) and the predicted one $\hat{p}_i$ by the linear regression using the selected $d$ features. The figure format follows that of Fig. 7.</figcaption></figure>

<footnote n="5">⁵S. Karmakar, C. Dasgupta, and S. Sastry, “Growing length scales and their relation to timescales in glass-forming liquids,” Annu. Rev. Condens. Matter Phys. **5**, 255–284 (2014).</footnote>

<footnote n="6">⁶A. Widmer-Cooper and P. Harrowell, “Predicting the long-time dynamic heterogeneity in a supercooled liquid on the basis of short-time heterogeneities,” Physical review letters **96**, 185701 (2006).</footnote>

<footnote n="7">⁷L. Berthier and R. L. Jack, “Structure and dynamics of glass formers: Predictability at large length scales,” Physical Review E—Statistical, Nonlinear, and Soft Matter Physics **76**, 041509 (2007).</footnote>

<footnote n="8">⁸A. Widmer-Cooper, H. Perry, P. Harrowell, and D. R. Reichman, “Irreversible reorganization in a supercooled liquid originates from localized soft modes,” Nature Physics **4**, 711–715 (2008).</footnote>

<footnote n="9">⁹G. Tarjus, S. A. Kivelson, Z. Nussinov, and P. Viot, “The frustration-based approach of supercooled liquids and the glass transition: a review and critical assessment,” Journal of Physics: Condensed Matter **17**, R1143 (2005).</footnote>

<footnote n="10">¹⁰D. Coslovich and G. Pastore, “Understanding fragility in supercooled lennard-jones mixtures. i. locally preferred structures,” The Journal of chemical physics **127** (2007).</footnote>

<footnote n="11">¹¹C. Patrick Royall, S. R. Williams, T. Ohtsuka, and H. Tanaka, “Direct observation of a local structural mechanism for dynamic arrest,” Nature materials **7**, 556–561 (2008).</footnote>

<footnote n="12">¹²G. M. Hocky, D. Coslovich, A. Ikeda, and D. R. Reichman, “Correlation of local order with particle mobility in supercooled liquids is highly system dependent,” Physical review letters **113**, 157801 (2014).</footnote>

<footnote n="13">¹³R. L. Jack, A. J. Dunleavy, and C. P. Royall, “Information-theoretic measurements of coupling between structure and dynamics in glass formers,” Physical review letters **113**, 095703 (2014).</footnote>

<footnote n="14">¹⁴F. Turci, C. P. Royall, and T. Speck, “Nonequilibrium phase transition in an atomistic glassformer: The connection to thermodynamics,” Physical Review X **7**, 031028 (2017).</footnote>

<footnote n="15">¹⁵C. P. Royall and S. R. Williams, “The role of local structure in dynamical arrest,” Physics Reports **560**, 1–75 (2015).</footnote>

<footnote n="16">¹⁶H. Tanaka, H. Tong, R. Shi, and J. Russo, “Revealing key structural features hidden in liquids and glasses,” Nature Reviews Physics **1**, 333–348 (2019).</footnote>

<footnote n="17">¹⁷A. Widmer-Cooper and P. Harrowell, “On the relationship between structure and dynamics in a supercooled liquid,” Journal of Physics: Condensed Matter **17**, S4025 (2005).</footnote>

<footnote n="18">¹⁸H. Tanaka, T. Kawasaki, H. Shintani, and K. Watanabe, “Critical-like behaviour of glass-forming liquids,” Nature materials **9**, 324–331 (2010).</footnote>

<footnote n="19">¹⁹H. Tong and H. Tanaka, “Revealing hidden structural order controlling both fast and slow glassy dynamics in supercooled liquids,” Physical Review X **8**, 011041 (2018).</footnote>

<footnote n="20">²⁰G. Kapteijns, D. Richard, E. Bouchbinder, T. B. Schrøder, J. C. Dyre, and E. Lerner, “Does mesoscopic elasticity control viscous slowing down in glassforming liquids?” The Journal of chemical physics **155** (2021).</footnote>

<footnote n="21">²¹L. Berthier, “Self-induced heterogeneity in deeply supercooled liquids,” Physical Review Letters **127**, 088002 (2021).</footnote>

<footnote n="22">²²E. Boattini, F. Smallenburg, and L. Filion, “Averaging local structure to predict the dynamic propensity in supercooled liquids,” Physical Review Letters **127**, 088007 (2021).</footnote>

<footnote n="23">²³G. Tarjus, “An overview of the theories of the glass transition,” Dynamical heterogeneities in glasses, colloids, and granular media **150**, 39 (2011).</footnote>

<footnote n="24">²⁴L. Berthier, M. Ozawa, and C. Scalliet, “Configurational entropy of glass-forming liquids,” The Journal of chemical physics **150** (2019).</footnote>

<footnote n="25">²⁵S. Sengupta, S. Karmakar, C. Dasgupta, and S. Sastry, “Adam-gibbs relation for glass-forming liquids in two, three, and four dimensions,” Physical review letters **109**, 095705 (2012).</footnote>

<footnote n="26">²⁶V. Lubchenko and P. G. Wolynes, “Theory of structural glasses and super-</footnote>
cooled liquids,” Annu. Rev. Phys. Chem. **58**, 235–266 (2007).

27J.-P. Bouchaud and G. Biroli, “On the adam-gibbs-kirkpatrick-thirumalai-
wolynes scenario for the viscosity increase in glasses,” The Journal of
chemical physics **121**, 7347–7354 (2004).

28J. C. Dyre, “Solid-that-flows picture of glass-forming liquids,” The Journal of Physical Chemistry Letters **15**, 1603–1617 (2024).

29E. D. Cubuk, S. S. Schoenholz, J. M. Rieser, B. D. Malone, J. Rottler, D. J. Durian, E. Kaxiras, and A. J. Liu, “Identifying structural flow defects in disordered solids using machine-learning methods,” Physical review letters **114**, 108001 (2015).

30S. S. Schoenholz, E. D. Cubuk, D. M. Sussman, E. Kaxiras, and A. J. Liu, “A structural approach to relaxation in glassy liquids,” Nature Physics **12**, 469–471 (2016).

31V. Bapst, T. Keck, A. Grabska-Barwińska, C. Donner, E. D. Cubuk, S. S. Schoenholz, A. Obika, A. W. Nelson, T. Back, D. Hassabis, et al., “Unveiling the predictive power of static structure in glassy systems,” Nature Physics **16**, 448–454 (2020).

32G. Jung, G. Biroli, and L. Berthier, “Predicting dynamic heterogeneity in glass-forming liquids by physics-inspired machine learning,” Physical Review Letters **130**, 238202 (2023).

33H. Shiba, M. Hanai, T. Suzumura, and T. Shimokawabe, “Botan: Bond targe-
ting network for prediction of slow glassy dynamics by machine learning
relative motion,” The Journal of Chemical Physics **158** (2023).

34F. S. Pezzicoli, G. Charpiat, and F. P. Landes, “Se (3)-equivariant graph neural networks for learning glassy liquids representations,” arXiv preprint arXiv:2211.03226 (2022).

35S. Ciarella, D. Khomenko, L. Berthier, F. C. Mocanu, D. R. Reichman, C. Scalliet, and F. Zamponi, “Finding defects in glasses through machine learning,” Nature Communications **14**, 4229 (2023).

36G. Janzen, C. Smit, S. Visbeek, V. E. Debets, C. Luo, C. Storm, S. Ciarella, and L. M. Janssen, “Classifying the age of a glass based on structural properties: A machine learning approach,” Physical Review Materials **8**, 025602 (2024).

37G. Janzen, X. L. Smeets, V. E. Debets, C. Luo, C. Storm, L. M. Janssen, and S. Ciarella, “Dead or alive: Distinguishing active from passive particles using supervised learning (a),” Europhysics Letters **143**, 17004 (2023).

38E. Boattini, S. Marín-Aguilar, S. Mitra, G. Foffi, F. Smallenburg, and L. Fil-
ion, “Autonomously revealing hidden local structures in supercooled liq-
uids,” Nature communications **11**, 5479 (2020).

39N. Oyama, S. Koyama, and T. Kawasaki, “What do deep neural networks find in disordered structures of glasses?” Frontiers in Physics **10**, 1320 (2023).

40D. Coslovich, R. L. Jack, and J. Paret, “Dimensionality reduction of local structure in glassy binary mixtures,” The Journal of Chemical Physics **157** (2022).

41A. Glielmo, C. Zeni, B. Cheng, G. Csányi, and A. Laio, “Ranking the information content of distance measures,” PNAS nexus **1**, pgac039 (2022).

42V. Del Tatto, G. Fortunato, D. Bueti, and A. Laio, “Robust inference of causality in high-dimensional dynamical processes from the information imbalance of distance ranks,” Proceedings of the National Academy of Sciences **121**, e2317256121 (2024).

43E. D. Donkor, A. Laio, and A. Hassanali, “Do machine-learning atomic de-
scriptors and order parameters tell the same story? the case of liquid water,”
Journal of Chemical Theory and Computation **19**, 4596–4605 (2023).

44E. D. Donkor, A. Offei-Danso, A. Rodriguez, F. Sciortino, and A. Has-
sanali, “Beyond local structures in critical supercooled water through un-
supervised learning,” The Journal of Physical Chemistry Letters **15**, 3996–
4005 (2024).

45R. Wild, E. Sozio, R. G. Margiotta, F. Dellai, A. Acquasanta, F. Del Ben, C. Tascini, F. Curcio, and A. Laio, “Maximally informative feature selection using information imbalance: Application to covid-19 severity prediction,” Scientific Reports **14**, 10744 (2024).

46H. Tong and H. Tanaka, “Emerging exotic compositional order on approaching low-temperature equilibrium glasses,” Nature Communications **14**, 4614 (2023).

47T. Kawasaki, T. Araki, and H. Tanaka, “Correlation between dynamic het-
erogeneity and medium-range order<? format?> in two-dimensional glass-
forming liquids,” Physical review letters **99**, 215701 (2007).

48M. L. Falk and J. S. Langer, “Dynamics of viscoplastic deformation in amorphous solids,” Physical Review E **57**, 7192 (1998).

49A. Barbot, M. Lerbinger, A. Hernandez-Garcia, R. García-García, M. L. Falk, D. Vandembroucq, and S. Patinet, “Local yield stress statistics in model amorphous solids,” Physical Review E **97**, 033001 (2018).

50M. Lerbinger, A. Barbot, D. Vandembroucq, and S. Patinet, “Relevance of shear transformations in the relaxation of supercooled liquids,” Physical Review Letters **129**, 195501 (2022).

51D. Frenkel and B. Smit, *Understanding molecular simulation: from algorithms to applications* (Elsevier, 2023).

52L. Berthier and W. Kob, “The monte carlo dynamics of a binary lennard-jones glass-forming mixture,” Journal of Physics: Condensed Matter **19**, 205130 (2007).

53W. Kob and H. C. Andersen, “Testing mode-coupling theory for a supercooled binary lennard-jones mixture. ii. intermediate scattering function and dynamic susceptibility,” Physical Review E **52**, 4134 (1995).

54H. Shiba, Y. Yamada, T. Kawasaki, and K. Kim, “Unveiling dimensionality dependence of glassy dynamics: 2d infinite fluctuation eclipses inherent structural relaxation,” Physical review letters **117**, 245701 (2016).

55B. Illing, S. Fritschi, H. Kaiser, C. L. Klix, G. Maret, and P. Keim, “Mermin–wagner fluctuations in 2d amorphous solids,” Proceedings of the National Academy of Sciences **114**, 1856–1861 (2017).

56S. Vivek, C. P. Kelleher, P. M. Chaikin, and E. R. Weeks, “Long-wavelength fluctuations and the glass transition in two dimensions and three dimensions,” Proceedings of the National Academy of Sciences **114**, 1850–1855 (2017).

57G. Tarjus, “Glass transitions may be similar in two and three dimensions, after all,” Proceedings of the National Academy of Sciences **114**, 2440–2442 (2017).

58R. Yamamoto and A. Onuki, “Dynamics of highly supercooled liquids: Heterogeneity, rheology, and diffusion,” Physical Review E **58**, 3515 (1998).

59H. Shiba, T. Kawasaki, and A. Onuki, “Relationship between bond-
breakage correlations and four-point correlations in heterogeneous glassy
dynamics: Configuration changes and vibration modes,” Physical Review
E **86**, 041504 (2012).

60C. F. Schreck, C. S. O’Hern, and L. E. Silbert, “Tuning jammed frictionless disk packings from isostatic to hyperstatic,” Physical Review E—Statistical, Nonlinear, and Soft Matter Physics **84**, 011305 (2011).

61R. M. Alkemade, E. Boattini, L. Filion, and F. Smallenburg, “Comparing machine learning techniques for predicting glassy dynamics,” The Journal of Chemical Physics **156** (2022).

62C. Bishop, *Pattern Recognition and Machine Learning* (Springer, 2006).

63G. Jung, R. M. Alkemade, V. Bapst, D. Coslovich, L. Filion, F. P. Landes, A. Liu, F. S. Pezzicoli, H. Shiba, G. Volpe, et al., “Roadmap on machine learning glassy liquids,” arXiv preprint arXiv:2311.14752 (2023).

64E. Flenner and G. Szamel, “Fundamental differences between glassy dynamics in two and three dimensions,” Nature communications **6**, 7392 (2015).

65A. Furukawa and H. Tanaka, “Significant difference in the dynamics between strong and fragile glass formers,” Physical Review E **94**, 052607 (2016).

66K. Kim and S. Saito, “Multiple length and time scales of dynamic heterogeneities in model glass-forming liquids: A systematic analysis of multi-point and multi-time correlations,” The Journal of chemical physics **138** (2013).

67C. Scalliet, B. Guiselin, and L. Berthier, “Thirty milliseconds in the life of a supercooled liquid,” Physical Review **X** **12**, 041028 (2022).

68M. L. Manning and A. J. Liu, “Vibrational modes identify soft spots in a sheared disordered packing,” Physical Review Letters **107**, 108302 (2011).

69S. A. Ridout and A. J. Liu, “The dynamics of machine-learned “softness” in supercooled liquids describe dynamical heterogeneity,” arXiv preprint arXiv:2406.05868 (2024).

70R. Sahu, M. Sharma, P. Schall, S. Maitra Bhattacharyya, and V. Chikkadi, “Structural origin of relaxation in dense colloidal suspensions,” Proceedings of the National Academy of Sciences **121**, e2405515121 (2024).

71E. Anwar, P. Patel, M. Sharma, and S. Maitra Bhattacharyya, “Exploring the soft pinning effect in the dynamics and the structure–dynamics correlation in multicomponent supercooled liquids,” The Journal of Chemical Physics **161** (2024).

72D. Richard, G. Kaptejns, and E. Lerner, “Detecting low-energy quasilocalized excitations in computer glasses,” Physical Review E **108**, 044124 (2023).
<footnote n="73">⁷³P. Charbonneau, E. Dyer, J. Lee, and S. Yaida, “Linking dynamical heterogeneity to static amorphous order,” Journal of Statistical Mechanics: Theory and Experiment **2016**, 074004 (2016).</footnote>

<footnote n="74">⁷⁴B. Guiselin, L. Berthier, and G. Tarjus, “Statistical mechanics of coupled supercooled liquids in finite dimensions,” SciPost Physics **12**, 091 (2022).</footnote>

<footnote n="75">⁷⁵M. Baggioli, I. Kriuchevskyi, T. W. Sirk, and A. Zaccone, “Plasticity in amorphous solids is mediated by topological defects in the displacement field,” Physical Review Letters **127**, 015501 (2021).</footnote>

<footnote n="76">⁷⁶Z. W. Wu, Y. Chen, W.-H. Wang, W. Kob, and L. Xu, “Topology of vibrational modes predicts plastic events in glasses,” Nature Communications **14**, 2955 (2023).</footnote>

<footnote n="77">⁷⁷D. Richard, M. Ozawa, S. Patinet, E. Stanifer, B. Shang, S. Ridout, B. Xu, G. Zhang, P. Morse, J.-L. Barrat, et al., “Predicting plasticity in disordered solids from structural indicators,” Physical Review Materials **4**, 113609 (2020).</footnote>

<footnote n="78">⁷⁸A. Lemaître, “Structural relaxation is a scale-free process,” Physical review letters **113**, 245702 (2014).</footnote>

<footnote n="79">⁷⁹M. Ozawa and G. Biroli, “Elasticity, facilitation, and dynamic heterogeneity in glass-forming liquids,” Physical Review Letters **130**, 138201 (2023).</footnote>

<footnote n="80">⁸⁰A. Tahaei, G. Biroli, M. Ozawa, M. Popović, and M. Wyart, “Scaling description of dynamical heterogeneity and avalanches of relaxation in glass-forming liquids,” Physical Review X **13**, 031034 (2023).</footnote>

<footnote n="81">⁸¹D. Chandler and J. P. Garrahan, “Dynamics on the way to forming glass: Bubbles in space-time,” Annual review of physical chemistry **61**, 191–217 (2010).</footnote>