# Interpretability of linear regression models of glassy dynamics

Anand Sharma, Chen Liu, Misaki Ozawa, Daniele Coslovich

► To cite this version:

Anand Sharma, Chen Liu, Misaki Ozawa, Daniele Coslovich. Interpretability of linear regression models of glassy dynamics. 2026. hal-05469895

HAL Id: hal-05469895

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

Preprint submitted on 21 Jan 2026

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 CC BY 4.0 - Attribution - International License
# Interpretability of linear regression models of glassy dynamics

Anand Sharma

Indian Institute of Science Education and Research, Dr. Homi Bhabha Road, Pashan, Pune 411008, India

Univ. Grenoble Alpes, CNRS, LIPhy, 38000 Grenoble, France and

Department of Physics and Materials Science, University of Luxembourg, L-1511 Luxembourg, Luxembourg

Chen Liu

Innovation and Research Division, Ge-Room Inc., 93160 Noisy le Grand, France

Misaki Ozawa

Univ. Grenoble Alpes, CNRS, LIPhy, 38000 Grenoble, France

Daniele Coslovich*

Dipartimento di Fisica, Università di Trieste, Strada Costiera 11, 34151, Trieste, Italy

(Dated: August 25, 2025)

Data-driven models can accurately describe and predict the dynamical properties of glass-forming liquids from structural data. Accurate predictions, however, do not guarantee an understanding of the underlying physical phenomena and the key factors that control them. In this paper, we illustrate the merits and limitations of linear regression models of glassy dynamics built on high-dimensional structural descriptors. By analyzing data for a two-dimensional glass model, we show that several descriptors commonly used in glass-transition studies display multicollinearity, which hinders the interpretability of linear models. Ridge regression suppresses some of the shortcomings of multicollinearity, but its solutions are not succinct enough to be physically interpretable. Only by using dimensional reduction techniques we eventually obtain linear models that strike a balance between prediction accuracy and interpretability. Our analysis points to a key role of local packing and composition fluctuations in the glass model under study.

## I. INTRODUCTION

As the temperature of glass-forming liquids decreases, the relaxation timescale increases dramatically, while structural changes remain relatively modest [1–3]. Moreover, the dynamics exhibits pronounced spatial fluctuations, characterized by regions of high and low mobility, referred to as dynamic heterogeneities [4–6]. Despite the vivid patterns of dynamic heterogeneities, however, static snapshots of the system appear homogeneous and lack distinct features, at least to the naked eye. This has been confirmed through direct observations in computer simulations and colloidal glass experiments [7, 8]. To identify subtle but significant structural changes linked to dynamics, various structural order parameters based on physical intuition have been investigated [9–18]. Similar efforts have been done to identify the structural origins of plastic events in amorphous solids under loading [19–22].

Recent advances in machine learning have demonstrated that glassy dynamics, including dynamic heterogeneities, can be accurately described and predicted from local structural information [23, 24]. A range of machine learning techniques has been applied to this problem, including support vector machines [23, 25], multi-layer perceptrons [26], and graph neural networks [27–30]. This research has primarily progressed in two directions. The first involves increasing the complexity of machine learning architectures, employing deep neural networks with a large number of parameters and taking advantage

of cutting-edge methods [27, 28, 30]. The second direction focuses on integrating domain knowledge from physics in data-driven models [26, 31]. These physics-informed approaches enable high predictive accuracy while maintaining relatively simple model architectures.

However, achieving accurate predictions alone does not guarantee an understanding of the underlying mechanisms driving the phenomenon under investigation [32, 33]. From the perspective of fundamental physics research, it is crucial that data-driven models provide physically interpretable results, which must be robust and expressed in a succinct form. A growing body of glass transition studies explores the issue of interpretability in deep-learning models or non-linear models, using a range of approaches [34–45]. While these works offer useful clues, the underlying deep networks remain highly complex and only partially transparent. Consequently, there is still a clear need for explicitly interpretable models whose solutions admit a direct physical reading.

Remarkably, recent studies have demonstrated that simple linear models, when combined with domain knowledge of glassy dynamics, such as coarse-graining techniques, can describe dynamic heterogeneities with remarkable accuracy [31, 46]. In some cases, their accuracy is comparable to that of more complex deep learning models. This outcome is particularly desirable, because it is often believed that simpler models not only provide accurate predictions at a minimal computational cost, but can also offer greater interpretability [47]. We caution, however, that employing a linear model *per se* is not sufficient to extract meaningful physical information: in high-dimensional settings, linear models are often plagued by numerical instabilities, due to the so-called multicollinearity [48], that can

<footnote>* Corresponding author: dcoslovich@units.it</footnote>
obscure interpretation. Moreover, their solution may still be too high-dimensional to convey a physical meaning. The aim of this study is to address these issues in linear models for glassy dynamics and introduce remedies that restore interpretability.

To tackle the issue of interpretability in a rigorous manner, we will study quantitatively the consequences of multicollinearity in linear regression models of glassy dynamics, using a two-dimensional glass-forming liquid model. First, we will demonstrate how multicollinearity leads to instability in weight estimation within linear regression, thereby hindering the interpretation of feature importance. These problems affect several structural descriptors used in recent glass transition studies. Next, we will explore strategies to mitigate the effects of multicollinearity and to identify a set of low-dimensional linear models that achieve good performance accuracy with a minimum of structural information. We will critically evaluate and discuss the advantages and limitations of each approach, providing a comprehensive assessment of linear models of the structure-dynamics relationship in a glass-forming liquid model.

The paper is organized as follows. Section II defines the problem under investigation. Section III describes the physical model and our structure-dynamics dataset. Sections IV and Section V introduce simple linear regression models of glassy dynamics and demonstrate the effects of multicollinearity on the instability of the estimated weights. Section VI examines feature selection and extraction approaches that cope with multicollinearity while also reducing the dimensionality of the problem. Sections VII and VIII offer a critical outlook on our results and summarize the key findings of the work.

## II. PROBLEM STATEMENT

In data-driven modeling of glassy dynamics, the problem is typically formulated as follows [24]. Several variables, called features, are computed in order to characterize the local structural environment around each particle. These structural features are collected into a descriptor that provides a high-dimensional representation of the local structure. Linear regression models estimate the dynamical observable $\mathbf{Y}$ (in this study, the dynamic propensity [49]) as a linear combination of $M$ input structural features, $\mathbf{X}^{(1)}, \mathbf{X}^{(2)}, \dots, \mathbf{X}^{(M)}$, using the following expression:

$$ \hat{\mathbf{Y}} = \hat{w}^{(1)}\mathbf{X}^{(1)} + \hat{w}^{(2)}\mathbf{X}^{(2)} + \dots + \hat{w}^{(M)}\mathbf{X}^{(M)}, \quad (1) $$

where $\hat{w}^{(1)}, \hat{w}^{(2)}, \dots, \hat{w}^{(M)}$ are the weights determined by minimizing a specific loss function.

How can we extract physical insights from the linear model in Eq. (1)? When $\mathbf{Y}$ and $\mathbf{X}^{(f)}$ ($f = 1, 2, \dots, M$) are properly normalized, for example, to have zero mean and unit variance, the sign and magnitude of the weights $\hat{w}^{(f)}$ provide information about the influence of the corresponding input feature $\mathbf{X}^{(f)}$ on the dynamical output $\mathbf{Y}$ (and its prediction $\hat{\mathbf{Y}}$). In other words, $\hat{w}^{(f)}$ serves as a measure of feature importance. Thus, linear regression models provide a strong case for physical interpretability, in which the connection between $\mathbf{X}$ and $\mathbf{Y}$ via the weights is direct and mechanistic.

In practice, however, the estimation of $\hat{w}^{(f)}$ can be unstable if $\mathbf{X}^{(1)}, \mathbf{X}^{(2)}, \dots, \mathbf{X}^{(M)}$ are strongly correlated with one another. As an extreme case, if $\mathbf{X}^{(1)}, \mathbf{X}^{(2)}, \dots, \mathbf{X}^{(M)}$ are indeed linearly dependent, the weights $\hat{w}^{(1)}, \hat{w}^{(2)}, \dots, \hat{w}^{(M)}$ are not uniquely determined. In real datasets, situations close to linear dependency frequently arise. Then, small perturbations in the dataset, such as numerical errors or limited statistics, can cause large variations in the weights. This pathological yet common phenomenon is known as multicollinearity [48].

Typically, multicollinearity is not a significant concern in machine learning studies, as the primary objective is to improve the prediction accuracy. However, it becomes problematic when interpreting regression models in terms of feature importance. In fact, if the weight estimation is unstable, interpretation becomes unreliable. The weight estimate must be stable against fluctuations in the dataset or small changes in the hyperparameters. Such robustness is, of course, a basic requirement for interpretability of data-driven models, see Ref. 50 for a recent review on this topic. In our view, interpretable models must also provide a succinct description of the relationships between the variables of interest. This is in line with the physicists' expectation that good models significantly compress the information involved in the problem at hand [51]. Phenomenological models in liquid state theory or statistical physics, for instance, are often based on a handful of independent variables – think of the hydrodynamic description of liquids [52], two-state models [53], etc. Hence, to be physically interpretable, data-driven models must be able to pinpoint a robust relationship between a small number of relevant variables, possibly carrying an intuitive physical meaning.

To close this section, a comment on terminology is in order. To evaluate the quality of a data-driven model, one computes the estimate $\hat{\mathbf{Y}}$ for a set of input data that were not used to train the model. In the following, we will do that by splitting the original data set into a training and test set, as is customary. The estimate $\hat{\mathbf{Y}}$ is then considered as prediction for the test set; standard goodness-of-fit metrics can be used to evaluate their quality. From a conventional physicist's viewpoint, this is a weak form of prediction: in the problem we will be dealing with, the data points will be sampled at precisely the same physical conditions both in the train and test set and no attempt will be made to extrapolate to different conditions, e.g., different temperatures or time scales. In this work, we will nonetheless stick to the word “prediction” for consistency with previous work [24] and with the widespread usage in the machine learning context. Should the readers feel uncomfortable with it, they can mentally replace “prediction” with “description” in the following.

## III. SIMULATION MODEL AND DATASET

### A. Simulation model

We use a three-component glass-forming liquid model composed of small (S), medium (M), and large (L) particles in two spatial dimensions with periodic boundary conditions [54]. The model is a variant of a well-studied binary mixture model [55–57], to which we add particles of type M with an intermediate
character between S and L [58]. The interaction between two particles is described by the Lennard-Jones potential

$$v_{\alpha\beta}(r) = 4\epsilon_{\alpha\beta} \left[ \left( \frac{\sigma_{\alpha\beta}}{r} \right)^{12} - \left( \frac{\sigma_{\alpha\beta}}{r} \right)^{6} \right],$$

where $\alpha, \beta = S, M, L$. The potential is modified to ensure that it is twice continuously differentiable at the cutoff, following Ref. [56].

The parameters $\sigma_{\alpha\beta}$ and $\epsilon_{\alpha\beta}$ are:

$$\begin{align*}
\sigma_{LL} &= 2 \sin\left(\frac{\pi}{5}\right) \approx 1.18, & \sigma_{SS} &= 2 \sin\left(\frac{\pi}{10}\right) \approx 0.62, & \sigma_{LS} &= 1, \\
\sigma_{LM} &= \frac{\sigma_{LL} + \sigma_{LS}}{2}, & \sigma_{MS} &= \frac{\sigma_{LS} + \sigma_{SS}}{2}, & \sigma_{MM} &= \frac{\sigma_{LL} + \sigma_{SS}}{2}, \\
\epsilon_{LL} &= \frac{1}{2}, & \epsilon_{SS} &= \frac{1}{2}, & \epsilon_{LS} &= 1, \\
\epsilon_{LM} &= \frac{\epsilon_{LL} + \epsilon_{LS}}{2}, & \epsilon_{MS} &= \frac{\epsilon_{LS} + \epsilon_{SS}}{2}, & \epsilon_{MM} &= \frac{\epsilon_{LL} + \epsilon_{SS}}{2}.
\end{align*}$$

The total 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 numbers of small, medium, and large particles, respectively. We use the NVT canonical ensemble, with a number density $\rho = N/L^2 = 1.024$, where $L$ is the linear length of the square simulation cell.

We perform Monte Carlo (MC) simulations using translational displacements [59]. The MC move consists in picking a particle at random and displacing it by a vector drawn randomly within a square box of linear size $\delta_{max} = 0.12$. The move is accepted on the basis of the Metropolis acceptance rule, which ensures the detailed balance condition. Although MC simulations do not possess a physical timescale, time $t$ can be measured in units of MC sweeps, each comprising $N$ attempts to perform the MC move. In the regime of slow glassy dynamics of interest in this work, the Monte Carlo dynamics behaves similarly to other types of physical dynamics, e.g., Newtonian and Brownian dynamics [60]. Therefore, we analyze the MC dynamics by following particle trajectories and calculating time-dependent observables as usual. The glassy dynamics of the model has been studied in Ref. [54] by computing the self intermediate scattering function. In this study, we focus on $T = 0.30$, which is the lowest temperature at which we can equilibrate the system within our computational timescale.

## B. Dynamic propensity

To investigate the heterogeneity of glassy dynamics in real space and assess its connection with the static structure, we use the iso-configurational ensemble [11, 49]. A set of $n = 10$ statistically uncorrelated configurations are obtained at equilibrium conditions at temperature $T = 0.30$. From each of these equilibrium configurations, we generate an ensemble of trajectories using MC dynamics at $T = 0.30$ using 30 different initial random seeds. We then compute the dynamic propensity

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

where $\langle \cdots \rangle_{\text{iso}}$ denotes an average over all the trajectories originating from the same initial configuration, and the cage-relative displacement, $\Delta \mathbf{r}_i^{\text{CR}}(t)$, is given by

$$\Delta \mathbf{r}_i^{\text{CR}}(t) = \Delta \mathbf{r}_i(t) - \frac{1}{n_i} \sum_{j \in N_i} \Delta \mathbf{r}_j(t),$$

where $\Delta \mathbf{r}_i(t) = \mathbf{r}_i(t) - \mathbf{r}_i(0)$ is the displacement vector of the $i$-th particle at position $\mathbf{r}_i$. Here, $n_i$ is the number of neighboring particles, and the set of neighbors ($N_i$) is defined as the particles located within a circular cutoff radius of $1.4\sigma_{\alpha\beta}$. The choice of cage-relative displacements is necessary to filter out the effect of the so-called Mermin-Wagner fluctuations in two-dimensional systems [61]. In this paper, we focus on the dynamic propensity computed at the structural relaxation timescale, $p_i(\tau_\alpha)$, where $\tau_\alpha$ is defined as the time at which the self intermediate scattering function becomes $1/e$ [54]. $\tau_\alpha$ at $T = 0.30$ is $\tau_\alpha \approx 4 \times 10^6$. We also considered a shorter time scale, $t = 5 \times 10^4$, which corresponds to $\beta$ relaxation timescale. The main findings of this work remained qualitatively unchanged.

## C. Behler-Parrinello descriptor

To characterize the local structure around each particle, we use the Behler-Parrinello (BP) descriptor [62], which has been widely used to study structure-property relationships, including the description of glassy dynamics in two dimensions [23, 63]. The descriptor comprises two subsets of features that characterize radial and angular correlations, respectively. Following Refs. [23, 63], we further distinguish features according to the species of the neighboring particles.

For each particle $i$, the radial feature $G_i^\alpha$ is defined by

$$G_i^\alpha = \sum_{j \in N_i} e^{-(r_{ij} - \mu)^2 / \delta^2} f_c(r_{ij}), \quad (2)$$

where $r_{ij} = |\mathbf{r}_i - \mathbf{r}_j|$ is the distance between particles $i$ and $j$, and $\mu$ and $\delta$ are parameters. The sum is carried out over the subset $N_\alpha$ of particles of species $\alpha$. $\sum'$ indicates that the particle $i$ is removed from the sum. The cut-off function $f_c(r)$ is defined by $f_c(r) = \frac{1}{2} [\cos(\pi r / R_c) + 1]$ for $r \le R_c$ and $f_c(r) = 0$ for $r > R_c$ [64]. $R_c$ is a cut-off radius and we set $R_c = 5.0\sigma_{LS}$. We vary $\mu$ between $0.3\sigma_{LS}$ and $5.0\sigma_{LS}$ in increments of $0.1\sigma_{LS}$ with $\delta = 0.1\sigma_{LS}$. Thus, for each species $\alpha$, the radial features $G^\alpha(k)$ are parametrized by an integer $k$ that selects values of the parameter $\mu = \mu_k = 0.3\sigma_{LS} + k \times 0.1\sigma_{LS}$ and $0 \le k \le 47$. Thus, for each particle, we have $144(= 3 \times 48)$ different radial features.

The angular descriptor $\Psi_i^{\alpha\beta}$ is defined by

$$\Psi_i^{\alpha\beta} = 2^{1-\zeta} \sum'_{\substack{j \in \alpha, k \in \beta \\ (j \neq k)}} e^{-(r_{ij}^2+r_{ik}^2+r_{jk}^2)/\xi^2} \times (1+\lambda \cos \theta_{ijk})^{\zeta} f_c(r_{ij}) f_c(r_{ik}) f_c(r_{jk}), \quad (3)$$

where $\theta_{ijk}$ is the angle at the corner $i$ of the triangle defined by particles $i, j$, and $k$, and $\xi, \lambda$, and $\zeta$ are parameters that
<table><caption>TABLE I. Parameters of the angular features $\Psi^{\alpha\beta}(k)$ of the BP descriptor.</caption><thead><tr><th></th><th>ξ</th><th>ζ</th><th>λ</th></tr></thead><tbody><tr><td>Ψ<sup>αβ</sup>(0)</td><td>14.633</td><td>1</td><td>-1</td></tr><tr><td>Ψ<sup>αβ</sup>(1)</td><td>14.633</td><td>1</td><td>1</td></tr><tr><td>Ψ<sup>αβ</sup>(2)</td><td>14.638</td><td>2</td><td>-1</td></tr><tr><td>Ψ<sup>αβ</sup>(3)</td><td>14.638</td><td>2</td><td>1</td></tr><tr><td>Ψ<sup>αβ</sup>(4)</td><td>2.554</td><td>1</td><td>-1</td></tr><tr><td>Ψ<sup>αβ</sup>(5)</td><td>2.554</td><td>1</td><td>1</td></tr><tr><td>Ψ<sup>αβ</sup>(6)</td><td>2.554</td><td>2</td><td>-1</td></tr><tr><td>Ψ<sup>αβ</sup>(7)</td><td>2.554</td><td>2</td><td>1</td></tr><tr><td>Ψ<sup>αβ</sup>(8)</td><td>1.648</td><td>1</td><td>1</td></tr><tr><td>Ψ<sup>αβ</sup>(9)</td><td>1.648</td><td>2</td><td>1</td></tr><tr><td>Ψ<sup>αβ</sup>(10)</td><td>1.204</td><td>1</td><td>1</td></tr><tr><td>Ψ<sup>αβ</sup>(11)</td><td>1.204</td><td>2</td><td>1</td></tr><tr><td>Ψ<sup>αβ</sup>(12)</td><td>1.204</td><td>4</td><td>1</td></tr><tr><td>Ψ<sup>αβ</sup>(13)</td><td>1.204</td><td>16</td><td>1</td></tr><tr><td>Ψ<sup>αβ</sup>(14)</td><td>0.933</td><td>1</td><td>1</td></tr><tr><td>Ψ<sup>αβ</sup>(15)</td><td>0.933</td><td>2</td><td>1</td></tr><tr><td>Ψ<sup>αβ</sup>(16)</td><td>0.933</td><td>4</td><td>1</td></tr><tr><td>Ψ<sup>αβ</sup>(17)</td><td>0.933</td><td>16</td><td>1</td></tr><tr><td>Ψ<sup>αβ</sup>(18)</td><td>0.695</td><td>1</td><td>1</td></tr><tr><td>Ψ<sup>αβ</sup>(19)</td><td>0.695</td><td>2</td><td>1</td></tr><tr><td>Ψ<sup>αβ</sup>(20)</td><td>0.695</td><td>4</td><td>1</td></tr><tr><td>Ψ<sup>αβ</sup>(21)</td><td>0.695</td><td>16</td><td>1</td></tr></tbody></table>

are varied systematically. For each pair of species, $(\alpha, \beta)$, we employ the same set of 22 parameters $(\xi, \lambda, \zeta)$ given in Ref. [23] in unit of $\sigma_{LS}$. The features $\Psi^{\alpha\beta}(k)$ are parametrized by an integer $k$ and the parameters $(\xi = \xi_k, \lambda = \lambda_k, \zeta = \zeta_k)$ with $0 \le k \le 21$ are shown in Table I. Thus, we have $132(=6 \times 22)$ angular features.

The full BP descriptor comprises a total of $M = 276(=144 + 132)$ features for each particle. Contrary to previous work [23, 63], we coarse-grain each feature over a length scale $\ell = 1.5$ using the procedure described in Sec. III D. Coarse-graining improves the prediction accuracy of the descriptor, without changing qualitatively the conclusions of this work. These variables constitute a feature vector, given by

$$ \mathbf{X}_i = (X_i^{(1)}, X_i^{(2)}, \dots, X_i^{(M)}). \quad (4) $$

We will sort the different kinds of features as follows:

$$ \mathbf{X}_i = (G_i^S, G_i^M, G_i^L, \Psi_i^{SS}, \Psi_i^{SM}, \Psi_i^{SL}, \Psi_i^{MM}, \Psi_i^{ML}, \Psi_i^{LL}). \quad (5) $$

### D. Physically motivated descriptors

To assess the generality of our findings, we also consider two additional structural descriptors that have been recently used to study structure-dynamics relationships in glass-forming liquids [26, 54]. Both of them are physically motivated: they are based on single-particle structural variables that characterize the environment around a particle in a physically intuitive way. All these single-particle variables are coarse-grained [31] over

multiple length scales, as described at the end of this section, to compose the full descriptor.

The local potential energy $u_i$ for particle $i$ is defined by

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

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

The coordination number $z_i$ for particle $i$ is defined as the number of neighboring particles within $r_{ij} < 1.5\sigma_{\alpha_i \beta_j}$, which corresponds well to the first minimum of each partial radial distribution function, $g_{\alpha\beta}(r)$.

The bond-orientational order parameter in two dimensions, $\Psi_{6,i}$, is defined by

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

where $\theta_{ij}$ is the angle between $\mathbf{r}_{ij} = \mathbf{r}_j - \mathbf{r}_i$ and the x-axis. The nearest neighbors are again defined as those within $r_{ij} < 1.5\sigma_{\alpha_i \beta_j}$. $\Psi_{6,i}$ quantifies hexagonal order, taking the value 1 for perfect hexagonal packings and smaller values for disordered packings [65, 66].

The steric bond order parameter $\Theta_i$ [67] is a measure of how well packed is the local environment around particle $i$. For each pair $\langle jk \rangle$ of neighboring particles, the angle $\theta_{jk}^{\text{ref}}$, calculated using the cosine formula. The steric order parameter is given by

$$ \Theta_i = \frac{1}{z_i} \sum_{\langle jk \rangle} |\theta_{jk} - \theta_{jk}^{\text{ref}}|, $$

where $\langle jk \rangle$ denotes the summation over all pairs of neighbors. Smaller values of $\Theta_i$ indicate sterically favored configurations, while larger values reflect disordered packings.

The local number density is given by

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

where $N_i$ includes particle $i$, $\ell$ is a coarse-graining length and all the other $N_i$ particles are included in the sum.

The local volume fraction is defined by

$$ \bar{\varphi}_i(\ell) = \sum_{j \in N_i} (\sigma_{\alpha_i \beta_j})^2 e^{-r_{ij}/\ell}. $$

Finally, we also consider the perimeter $\pi_i$ of the Voronoi cell surrounding particle $i$, as obtained from a radical Voronoi tessellation [68], using the nominal interaction parameters $\sigma_{\alpha\alpha}$ as particle radii.

With these structural features at hand, we can proceed to define two physically motivated descriptors. First, we define the SLO descriptor from the paper by Sharma, Liu, and Ozawa [54]. The SLO descriptor comprises $\bar{\rho}_i$, $\bar{\varphi}_i$, $u_i$, $z_i$, $\Psi_{6,i}$, and $\Theta_i$. Each
structural feature, $x_i = u_i, z_i, \Psi_{6,i}, \Theta_i$, is coarse-grained using the procedure

$$ \bar{x}_i(\ell) = \frac{1}{\bar{\rho}_i(\ell)} \sum_{j \in N_i} x_j e^{-r_{ij}/\ell}, \quad (6) $$

yielding $\bar{u}_i(\ell), \bar{z}_i(\ell), \bar{\Psi}_{6,i}(\ell)$, and $\bar{\Theta}_i(\ell)$. The coarse-graining length scale $\ell$ is varied from $0.5\sigma_{LS}$ to $5.0\sigma_{LS}$. Each kind of physically motivated feature $X(k)$ is parametrized by an integer $k$, yielding features with $\ell = \ell(k) = 0.5\sigma_{LS} + k \times 0.5\sigma_{LS}$ and $0 \le k \le 9$. The same lengths are used for the calculation of $\bar{\rho}_i(\ell)$ and $\bar{\varphi}_i(\ell)$. This yields a total of $M = 60$ features.

We also use the JBB descriptor introduced in the paper by Jung, Biroli, and Berthier [26]. This descriptor is based on $\bar{\rho}_i$, $u_i$, and $\pi_i$. In this work, we ignore the variance of the potential energy, which was included in Ref. [26]. The JBB descriptor incorporates particle-species information more explicitly, as shown below. A first set of coarse-grained features is obtained using the same procedure as for the SLO descriptor, considering the whole set of neighbors, irrespective of their species, with 10 different coarse-graining length scales $\ell$. In addition, the JBB descriptor includes coarse-grained features obtained by a procedure similar to Eq. (6) but taking species into account:

$$ \bar{x}_i^\alpha(\ell) = \frac{1}{\bar{\rho}_i^\alpha(\ell)} \sum_{j \in N_i^\alpha} x_j e^{-r_{ij}/\ell}, \quad (7) $$

where

$$ \bar{\rho}_i^\alpha(\ell) = \sum_{j \in N_i^\alpha} e^{-r_{ij}/\ell}, $$

and the sums are restricted to neighbors of species $\alpha$. Thus, in addition to the species-independent coarse-graining defined by Eq. (6), we consider coarse-graining based on Eq. (7) using the three types of particles ($\alpha = S, M$, and L). This yields a total of $M = 120$ features (3 descriptors $\times 4$ types $\times 10$ length scales).

## E. Dataset

The dataset we use in this work comprises both dynamic and structural information. The propensity is computed using all the $n$ available configurations. The structural descriptors are instead computed on the inherent structures [69] of the initial configurations used for the propensity calculations. The dataset is then composed of one of the structural descriptors defined in Secs. III C and III D and the dynamic propensity.

As commonly done in machine learning studies, the dataset is feature scaled. In particular, we normalize the propensity data,

$$ Y_i = \frac{p_i - \mathbb{E}[p]}{\sqrt{\mathrm{Var}[p]}}, $$

where $\mathbb{E}[\cdot]$ and $\mathrm{Var}[\cdot]$ denote the mean value and the variance, respectively. This normalization ensures that $Y_i$ has zero mean and a standard deviation of one.

To incorporate the dynamic data for $N_S$ particles, possibly taken from several different configurations, we will use a vector notation,

$$ \mathbf{Y} = [Y_1, Y_2, \dots, Y_{N_S}]^T, \quad (8) $$

where the superscript $T$ is the transpose operation.

Each structural feature of a given descriptor is also normalized to have zero mean and unit variance. After normalization, the features form a vector $\mathbf{X}_i$ for particle $i$:

$$ \mathbf{X}_i = [X_i^{(1)}, X_i^{(2)}, \dots, X_i^{(M)}]^T. \quad (9) $$

This vector serves as the structural input for the regression models.

To manage the whole dataset, it is convenient to introduce the following ($N_S \times M$) matrix

$$ X = \begin{bmatrix} X_1^{(1)} & X_1^{(2)} & \dots & X_1^{(M)} \\ X_2^{(1)} & X_2^{(2)} & \dots & X_2^{(M)} \\ \vdots & \vdots & \ddots & \vdots \\ X_{N_S}^{(1)} & X_{N_S}^{(2)} & \dots & X_{N_S}^{(M)} \end{bmatrix}, \quad (10) $$

which is sometimes called the design matrix.

Using Eq. (9), $X$ can be written as

$$ X = [\mathbf{X}_1, \mathbf{X}_2, \dots, \mathbf{X}_{N_S}]^T. \quad (11) $$

Thus, $\mathbf{X}_i$ ($i = 1, 2, \dots, N_S$) are the row vectors of $X$. In this paper, we also use the column vectors $\mathbf{X}^{(f)}$ ($f = 1, 2, \dots, M$) of $X$:

$$ X = [\mathbf{X}^{(1)}, \mathbf{X}^{(2)}, \dots, \mathbf{X}^{(M)}]. \quad (12) $$

The row and column vectors can be distinguished by the subscript and superscript.

In supervised learning studies, it is customary to train the model on one portion of the dataset and then test its performance on a different portion by computing some measure of correlation or statistical error. While this aspect is less crucial for linear models than for complex deep neural networks, we follow the standard procedure of splitting the full dataset into training and test sets. This is done by selecting a random subset of particles, $S_{train}$, as training set, from the full dataset. The fraction of selected particles in the training set is $x_{train}$. The test set is defined by selecting a random subset of particles, $S_{test}$, from the particles not included in $S_{train}$. The fraction of selected particles in the test set is, in general, $x_{test} \le x_{train}$.

For the supervised learning methods studied in Sec. IV we found that the goodness-of-fit metrics converge when the number of datapoints per feature is above $\approx 30$, which for our dataset corresponds to about $x_{train} = x_{test} = 0.2$. In the following, we consider $x_{train} = x_{test} = 0.5$. We use from 10 to 100 independent realizations of these sets to perform averages and assess the statistical accuracy of our results. As an exception, the analysis of the principal component regression in Sec. VI B involves the whole dataset.
<figure><img src="image_2.png" /><figcaption>FIG. 1. Pearson coefficient, $R[X^{(f)}, Y]$, between the dynamic propensity $Y$ and each structural feature $X^{(f)}$ of the BP descriptor for $f = 1, \dots, M$.</figcaption></figure>

## F. Pearson correlation coefficients

To illustrate some of key features of our datasets, we start with a simple analysis of correlations. To quantify the linear dependence between two variables, $A$ and $B$, computed for a subset $S$ of particles, we use the Pearson coefficient,

$$R[A, B] = \frac{1}{N_S} \sum_{i \in S} \frac{(A_i - \mathbb{E}[A])(B_i - \mathbb{E}[B])}{\sqrt{\mathrm{Var}[A] \mathrm{Var}[B]}}, \quad (13)$$

computed for a set $S$ comprising $N_S$ particles. By construction, $R[A, B]$ takes values in the range $-1 \le R[A, B] \le 1$.

We first compute the Pearson coefficients $R[X^{(f)}, Y]$ between the dynamic propensity $Y$ and each structural feature $X^{(f)}$ of the BP descriptor. Because of the very small dependence of the dynamic propensity on the species of the particles, we include all the particles, irrespective of their species, in the analysis of correlations. The Pearson coefficients are conveniently assembled in vector form

$$R[X, Y] = [R[X^{(1)}, Y], R[X^{(2)}, Y], \dots, R[X^{(M)}, Y]]^T. \quad (14)$$

In Fig. 1, we show results for $R[X, Y]$ obtained using the BP descriptor. The ordering of the features follows the logic described in Sec. III C. We first separate the features into a radial and an angular sector. Then, within each of these two sectors, the features are grouped into blocks according to the species or species pair for which the feature is computed. We see that the correlations are generally modest ($|R| < 0.5$) and spread over the whole set of structural features. In particular, both radial and angular features can have similar correlations or anticorrelations with the dynamic propensity. The presence of minima and maxima in each block of the radial features $G^\alpha$ can be connected to the peaks of the partial radial distribution functions, although we could not identify a straightforward interpretation. The angular sector displays systematic effects due to the chemical composition of the particle environment. Namely, features involving SS and LL pairs of neighbors are more strongly correlated or anti-correlated to the propensity. By contrast, $R$ shows no clear trend within each of the blocks of the angular sector.

The results for the two physically motivated descriptors introduced in Sec. III D display qualitatively similar trends, see the Appendix A. Among all the individual structural features entering our datasets, the steric order parameter $\bar{\Theta}_i(\ell)$, coarse-grained at $\ell = 2.0$, exhibits the strongest correlation ($R \approx 0.65$) with the dynamic propensity, making it the best-performing single feature within this simple correlation analysis. While the local volume fraction $\bar{\varphi}_i(\ell)$ and the coordination number $\bar{z}_i(\ell)$ are negatively correlated with the dynamic propensity, as expected, the local number density $\bar{\rho}_i(\ell)$ is positively correlated to it. This counter-intuitive result is likely a non-trivial effect of composition fluctuations in our ternary glass model.

<figure><img src="image_3.png" /><figcaption>FIG. 2. Correlation matrix C for the BP descriptor. The matrix elements are given by the Pearson correlation coefficient $R[X^{(f)}, X^{(f')}]$.</figcaption></figure>

A common feature of all the above descriptors is the presence of significant cross-correlations between groups of structural features. To quantify this effect, we introduce a correlation matrix, $C$, whose elements are given by

$$C_{f,f'} = R[X^{(f)}, X^{(f')}]. \quad (15)$$
<figure><img src="image_4.png" /><figcaption>FIG. 3. Weights obtained from OLS and Ridge regression of the dynamic propensity using the BP descriptor: (a) $\hat{w}_{OLS}$ and (b) $\hat{w}_{ridge}$ for $\alpha = 10^{-1}$. The error bar corresponds to the standard deviation estimated over independent random training sets.</figcaption></figure>

In Fig. 2, we show the correlation matrix for the BP descriptor. It has a distinct block structure, displaying strong positive or negative correlations within subsets of features. Correlations within species-wise blocks are particularly pronounced for the angular features, see the bottom right region of Fig. 2. Moreover, there are non-trivial cross-correlations also between radial and angular features. A similar block structure is also apparent in the physically motivated descriptors, see the Appendix A. In those cases, however, block correlations are due to coarse-graining a given structural feature over a range of similar distances $l$. As we will demonstrate in Sec. IV, the redundancy of these structural descriptors has a significant negative impact on the linear modeling of the dynamic propensity.

## IV. LEAST SQUARE REGRESSION

In this section, we introduce the simplest linear regression model to describe the dynamic propensity $\mathbf{Y}$ using a structural descriptor $\mathbf{X}$. The model yields a prediction of the dynamic propensity of the $i$-th particle based on the feature vector $\mathbf{X}_i$,

$$ \hat{Y}_i = \hat{\mathbf{w}}^T \mathbf{X}_i = \sum_{f=1}^{M} \hat{w}^{(f)} X_i^{(f)}, \quad (16) $$

where $\hat{\mathbf{w}} = [\hat{w}^{(1)}, \hat{w}^{(2)}, \dots, \hat{w}^{(M)}]^T$ are the weights.

With the vector and matrix notations in Eqs. (8) and (10), Eq. (16) is rewritten as

$$ \hat{\mathbf{Y}} = X\hat{\mathbf{w}}. \quad (17) $$

We consider linear models obtained by minimizing a loss function of the form

$$ \mathcal{L}(\hat{\mathbf{w}}) = \mathcal{L}^{\text{MSE}}(\hat{\mathbf{w}}) + \mathcal{L}^{\text{reg}}(\hat{\mathbf{w}}), \quad (18) $$

where

$$ \mathcal{L}^{\text{MSE}}(\hat{\mathbf{w}}) = \frac{1}{2N_S} \sum_{i \in S} (\hat{Y}_i - Y_i)^2 = \frac{1}{2N_S} \|\hat{\mathbf{Y}} - \mathbf{Y}\|^2 \quad (19) $$

is the mean square error (MSE) and $\mathcal{L}^{\text{reg}}$ is a regularization term. In this section, we will show results for two standard linear models [70], namely ordinary least square regression (Sec. IV A) and ridge regression (Sec. IV B), and discuss their limitations.

### A. Ordinary least squares regression

#### 1. Definition

In ordinary least square (OLS) regression, the weight vector $\hat{\mathbf{w}}$ is determined by minimizing the MSE, i.e., $\mathcal{L} = \mathcal{L}^{\text{MSE}}$. Setting $\nabla_{\hat{\mathbf{w}}} \mathcal{L}^{\text{MSE}}(\hat{\mathbf{w}}) = \mathbf{0}$ and remembering that all the features are normalized, the solution for the OLS regression (see the Appendix B for the derivation) is

$$ \hat{\mathbf{w}}_{\text{OLS}} = C^{-1} \mathbf{R}[X, Y], \quad (20) $$

where $C$ is the correlation matrix,

$$ C = \frac{1}{N_S} X^T X \quad (21) $$

and

$$ \mathbf{R}[X, Y] = \frac{1}{N_S} X^T \mathbf{Y} \quad (22) $$

is the vector of the Pearson coefficients.
Clearly, if all features were orthogonal to each other, $C = I$ ($I$ is the identity matrix) and $\hat{w}_{\text{OLS}}^{(f)}$ would equal the Pearson coefficient $R[X^{(f)}, Y]$ between the dynamic propensity and feature $f$. Then, the importance of a feature in linear regression (namely, the weight $\hat{w}_{\text{OLS}}^{(f)}$) would directly match its correlation with the propensity, as in the simple analysis of Sec. III F. Because of the presence of correlations between features, however, $C^{-1}$ possesses non-zero off-diagonal elements, and hence $\hat{w}_{\text{OLS}}^{(f)}$ is given by a linear combination of $R[X^{(1)}, Y], \dots, R[X^{(M)}, Y]$.

## 2. Oscillatory behavior of the weights

Figure 3(a) shows the weights $\hat{w}_{\text{OLS}}$ obtained for the BP descriptor. By comparing these results with Fig. 1, we immediately notice that the OLS solution gives a strong weight to the angular features. Moreover, $\hat{w}_{\text{OLS}}^{(f)}$ oscillates significantly, especially in the angular sector of the descriptor: the sign often changes considerably between successive features within a block. This consistent oscillatory behavior<fnref n="1" /> hampers the extraction of any meaningful physical insight into the relationship between the dynamical output $Y$ and the static features. In fact, a positive (negative) weight indicates that an increase in this feature enhances (reduces) the dynamic propensity. It would be hard to accept that very similar features, aligned along the x-axis in Fig. 3, can have large opposite effects on the dynamical variable $Y$.

We anticipate from Eq. (20) that the oscillations of $\hat{w}_{\text{OLS}}^{(f)}$ are related to the singularity of the matrix $C$. Indeed, the invertibility of the correlation matrix $C$ and the linear dependence of features are tightly connected. As an extreme case, one can show that the columns of the matrix $X = [X^{(1)}, \dots, X^{(M)}]$ are linearly independent if and only if the correlation matrix $C$ is invertible. Conversely, when some features are linearly dependent, $C$ is not invertible, and hence $\hat{w}_{\text{OLS}}^{(f)}$ is not uniquely determined. The results shown in Fig. 3(a) are representative of the instability due to strong correlations among features, a phenomenon known as multicollinearity in statistical analysis [48]. We will quantify it in detail in Sec. V.

Interestingly, however, the Pearson coefficient between the ground truth $Y$ and the prediction $\hat{Y}$ of OLS regression is very good, $R[Y, \hat{Y}] \approx 0.89$, on either the train or the test datasets. These results represent a paradigmatic machine learning case in which the prediction accuracy is very good, but the interpretability is poor. We confirmed that the other descriptors, namely the SLO and JBB descriptors, also demonstrate this pathological behavior.

<figure><img src="image_5.png" /><figcaption>FIG. 4. Normalized weights as a function of $\alpha$ corresponding to (a) the radial features $G^S(k)$ for $0 \le k \le 47$ and (b) the angular features $\Psi^{SS}(k)$ for $0 \le k \le 21$. The width of the shaded areas corresponds to the standard deviation estimated over independent random training sets. The color code indicates the feature index $k$.</figcaption></figure>

## B. Ridge regression

### 1. Definition

We now turn our attention to the so-called Ridge regression method [70], a simple variant of least square regression commonly used to alleviate the shortcomings of OLS. The loss function for Ridge regression, $\mathcal{L}^{\text{Ridge}}(\hat{\mathbf{w}})$,

$$ \mathcal{L}^{\text{Ridge}}(\hat{\mathbf{w}}) = \mathcal{L}^{\text{MSE}}(\hat{\mathbf{w}}) + \frac{\alpha}{2} \sum_{f=1}^{M} (\hat{w}^{(f)})^2, \quad (23) $$

includes a regularization term that penalizes solutions with large weights. The parameter $\alpha$ controls the magnitude of the regularization. Setting $\nabla_{\hat{\mathbf{w}}} \mathcal{L}^{\text{Ridge}}(\hat{\mathbf{w}}) = \mathbf{0}$ yields the estimated weights (see the Appendix B for the derivation)

$$ \hat{\mathbf{w}}_{\text{Ridge}} = (C + \alpha I)^{-1} \mathbf{R}[X, Y]. \quad (24) $$

Of course, when $\alpha \to 0$, $\hat{\mathbf{w}}_{\text{Ridge}}$ reduces to $\hat{\mathbf{w}}_{\text{OLS}}$. On the other hand, when $\alpha \gtrsim \lambda_{\max}$, where $\lambda_{\max}$ is the largest eigenvalue of $C$, we have $\hat{\mathbf{w}}_{\text{Ridge}} \approx \alpha^{-1} \mathbf{R}[X, Y]$, because the $\alpha I$ term dominates. In this regime, $\hat{\mathbf{w}}_{\text{Ridge}}$ cannot account for non-trivial correlations in $C$ and vanishes when $\alpha \to \infty$. Therefore, $\alpha$ should be chosen smaller than an order of $\lambda_{\max}$.

<footnote n="1">¹ The oscillations are not due to limited statistics or dataset-to-dataset fluctuations, since the error bars are relatively small compared to the typical absolute value of the largest weights.</footnote>
## 2. Sensitivity to the regularization parameter

To show the impact of the regularization, we show in Fig. 3(b) $\hat{w}_{Ridge}$ obtained for $\alpha = 0.1$. Two main effects are observed: the amplitude of the weights is more balanced between the radial and angular sectors compared to the OLS case, and the large oscillations between successive features are suppressed. While the main shortcoming of OLS regression seems to be solved, we notice that several features contribute with a finite weight, i.e., the solution is not sparse enough to be clearly interpretable. Note that by increasing $\alpha$ further toward $\lambda_{max} \approx 100$, we find that $\hat{w}_{Ridge} \approx \alpha^{-1}R[X, Y]$, as anticipated above.

How sensitive are the results to the regularization parameter? Analysis of the traces of the weights as a function of $\alpha$ provides a simple and intuitive way to assess the stability of the solutions in Ridge regression [71]. The idea is that as $\alpha$ increases from zero to small but finite values the weights will first change substantially, but there may be a range of $\alpha$ where $\hat{w}_{Ridge}$ they do not depend on $\alpha$ anymore. This, of course, should occur before reaching the trivial regime where $\hat{w}_{Ridge} \approx \alpha^{-1}R[X, Y]$.

To illustrate the results of this analysis, we show in Fig. 4 the Ridge traces corresponding to the features $G^S(k)$ and $\Psi^{SS}(k)$ in the top and bottom panels, respectively. In each panel and for each value of $\alpha$, the weights are scaled by the largest absolute value among the weights of the corresponding subset of features. This normalization adsorbs the huge change of scale evident from Fig. 3 and provides a vivid image of the sensitivity of the Ridge solutions. We find that the weights are extremely sensitive to $\alpha$, especially in the angular sector of the descriptor: the traces change chaotically over a broad range of $\alpha$ spanning several orders of magnitude. Only for $\alpha \ge 0.1$ the results seem to stabilize, in the sense that the order in the amplitude of the weights change less dramatically. Note that the chaotic behavior of the traces is not due to limited statistics, since the estimated error bars visible in the figure are relatively small.

## 3. Prediction accuracy

Crucially, very different solutions of the regression problem, corresponding to different $\alpha$, yield predictions for the dynamic propensity with nearly identical accuracy. To evaluate the prediction accuracy of the Ridge regression models, we use two standard metrics of performance as a function of $\alpha$. In addition to the Pearson coefficient, we also compute the coefficient of determination $R_2$, which focuses on the normalized total squared deviations. If $A$ is the variable to predict, $R_2[A, \hat{A}]$ is given by

$$R_2[A, \hat{A}] = 1 - \frac{\sum_{i \in S} (A_i - \hat{A}_i)^2}{\sum_{i \in S} (A_i - \mathbb{E}[A])^2} \quad (25)$$

When the prediction is perfect, $R_2[A, \hat{A}] = 1$. Conversely, when the prediction always outputs the mean value, $R_2[A, \hat{A}] = 0$, which serves as a baseline.

Figure 5 displays the performance metrics for the BP descriptor. We show results obtained on both the test set or the train set,

<figure><img src="image_6.png" /><figcaption>FIG. 5. Prediction performance metrics for Ridge regression of the dynamic propensity using the BP descriptor: (a) Pearson coefficient $R[Y, \hat{Y}]$ as a function of $\alpha$ and (b) coefficient of determination $R_2[Y, \hat{Y}]$ as a function of $\alpha$. Full and dashed lines correspond to results obtained using the test set and the train set, respectively.</figcaption></figure>

to verify the lack of overfitting. Both metrics indicate very good performance accuracy and display a flat maximum stretching over several orders of magnitudes in $\alpha$. The performance starts to degrade appreciably only when $\alpha > 0.1$, as can be seen from the similar drops of $R$ and $R_2$. We found quantitatively similar results for the physically motivated descriptors introduced in Sec. III D (not shown). The comparison between the chaotic traces of Fig. 4 and the features performance metrics of Fig. 5 is striking. Since very different solutions (see Fig. 4) yield nearly identical prediction accuracies (see Fig. 5), it is clear that prediction alone cannot be taken as a criterion to choose $\alpha$. This delicate aspect was overlooked in previous studies [31, 46], where the regularization parameter $\alpha$ was chosen to maximize the correlation. These issues are a manifestation of multicollinearity in the dataset and leave us with the question how to define an optimal model.

# V. MULTICOLLINEARITY AND ITS RESOLUTION

In this section, we introduce a simple metric that quantifies multicollinearity in the context of linear regression, namely the condition number. We first illustrate its qualitative behavior with a schematic two-features model (Sec. V A) and then use it to revisit the linear regression models for the dynamic propensity, including Ridge regression (Sec. V B). The output of this analysis provides us with a range of optimal weights for Ridge regression, solving the issue of multicollinearity in a statistical sense [48].

## A. Condition number

### 1. Definition

We start by quantifying the degree of multicollinearity using the condition number $\kappa(C)$ of the correlation matrix $C$. A brief introduction to the condition number and its meaning is provided in the Appendix C. In a nutshell, $\kappa(C)$ provides an upper bound for the relative error in the solution of a linear problem when small statistical perturbations are present in
the data. When $\kappa(C)$ is small, C is well-conditioned and the solution is stable. When $\kappa(C)$ is large, the matrix C is ill-conditioned and the error in the solution may become significant.

Since C is a positive semi-definite matrix, $\kappa(C)$ can be defined as

$$ \kappa(C) = \frac{\lambda_{\max}}{\lambda_{\min}}, \qquad (26) $$

where $\lambda_{\max}$ and $\lambda_{\min}$ are the largest and smallest eigenvalues of C, respectively. The instability of the matrix C is then identified by the divergence of its condition number. The largest eigenvalue is bounded as $\lambda_{\max} \le M$, because Tr(C) = $\sum_{f=1}^{M} \lambda^{(f)} = M$, where $\lambda^{(f)}$ are the eigenvalues. Therefore, the divergent behavior of $\kappa(C)$ arises from the vanishing of the smallest eigenvalue, $\lambda_{\min} \to 0$.

## 2. Two-features model

To understand the qualitative behavior of the condition number, we consider a simple toy model composed of only two features, $X = [\mathbf{X}^{(1)}, \mathbf{X}^{(2)}]$. The correlation matrix C is given by

$$ C = \begin{bmatrix} 1 & r \\ r & 1 \end{bmatrix}, \qquad (27) $$

where $r = R[X^{(1)}, X^{(2)}]$ is the Pearson coefficient between the two features. When $r = 0$, the two features are orthogonal, and C reduces to the identity matrix. When $r \to 1$ ($r \to -1$), the two features are perfectly correlated (anti-correlated), and hence $\mathbf{X}^{(1)}$ and $\mathbf{X}^{(2)}$ are linearly dependent.

Since we are interested in situations near the instability, we restrict ourselves to $0 < r \le 1$. In this setting, the eigenvalues of C are given by $\lambda_{\max} = 1+r$ and $\lambda_{\min} = 1-r$. Therefore, the condition number is

$$ \kappa(C) = \frac{1+r}{1-r}. \qquad (28) $$

When $r \to 1$, $\lambda_{\min} \to 0$ and hence the condition number diverges. Simultaneously, C becomes non-invertible since the determinant vanishes. Thus, $\kappa(C)$ effectively signals the presence of multicollinearity.

## 3. Origin of the oscillatory behavior of weights

We now seek to explain the oscillatory behavior of the weights observed in Fig. 3. This behavior arises from the ill-conditioning of the correlation matrix, due to the vanishing of its smallest eigenvalue $\lambda_{\min}$.

The key observation is that two strongly correlated features tend to acquire weights of opposite signs with large magnitude, even though such features are expected to be physically similar. This hampers the physical interpretation of the weights, which should reflect the importance of features contributing to the dynamics. The essence of this problem can be reduced to

<figure><img src="image_7.png" /><figcaption>FIG. 6. Condition number $\kappa(C + \alpha I)$ as a function of the Ridge regularization parameter $\alpha$. The vertical arrow marks the largest eigenvalue of C. The horizontal arrow is drawn at a condition number equal to 1000. The dashed line indicates an inverse power law.</figcaption></figure>

a two-feature model. Suppose that $R[X, Y] = [R^{(1)}, R^{(2)}]^T$, where $R^{(1)}$ and $R^{(2)}$ are some values with $-1 \le R^{(1)}, R^{(2)} \le 1$. Using Eq. (20), the estimated weights $\hat{w}_{\text{OLS}} = [\hat{w}_{\text{OLS}}^{(1)}, \hat{w}_{\text{OLS}}^{(2)}]^T$ are given by

$$ \hat{w}_{\text{OLS}}^{(1)} = \frac{1}{\lambda_{\max} \lambda_{\min}} (R^{(1)} - rR^{(2)}), \qquad (29) $$

$$ \hat{w}_{\text{OLS}}^{(2)} = -\frac{1}{\lambda_{\max} \lambda_{\min}} (rR^{(1)} - R^{(2)}). \qquad (30) $$

When $r \approx 1$, the magnitude of both weights becomes very large due to $\lambda_{\min} \to 0$, and the signs become opposite, $\hat{w}_{\text{OLS}}^{(1)} \approx -\hat{w}_{\text{OLS}}^{(2)}$.

## B. Least square regression revisited

We now turn our attention to the structural dataset given by the BP descriptor and quantify its degree of multicollinearity. A direct calculation of the condition number yields $\kappa(C) \approx 1.4 \times 10^{18}$. As a rule of thumb [48], condition numbers smaller than 100 are not problematic, while values greater than 1000 indicate problems with multicollinearity. We found similar orders of magnitude for the structural datasets based on the physically motivated descriptors, i.e., $4.5 \times 10^{17}$ and $1.2 \times 10^{15}$ for the JBB and SLO descriptors, respectively. Also the original BP descriptor used in Refs. [23, 63], in which features were not coarse-grained, yields a condition number of order $10^{17}$. We conclude that several structural datasets recently used in several recent glassy materials studies [23, 26, 54, 63] are severely affected by multicollinearity, much more than typical datasets analyzed in statistics textbooks [48].

We can now build on the multicollinearity analysis to identify a range of optimal solutions within Ridge regression. To see how the regularization term in Eq. (23) mitigates the effect of collinearity, we compute the condition number in the Ridge setting. We get

$$ \kappa(C + \alpha I) = \frac{\lambda_{\max} + \alpha}{\lambda_{\min} + \alpha}, \qquad (31) $$
which avoids the divergence of $\kappa(C+\alpha I)$ when multicollinearity is severe ($\lambda_{min} \to 0$). In Fig. 6(a), we show $\kappa(C + \alpha I)$ as a function of $\alpha$. In the range $\lambda_{min} \ll \alpha \ll \lambda_{max}$, the condition numbers decreases like $\lambda_{max}\alpha^{-1}$. To reduce the effects of multicollinearity to an acceptable degree, one should lower the condition numbers down to at least $10^3$ or less [48]. This occurs around $\alpha = 0.1$. The regularization becomes meaningless for $\alpha \ge \lambda_{max}$, since it entirely suppresses correlations in C. The crossover of $\kappa$ towards 1 occurs when $\alpha$ becomes of the order of $\lambda_{max} \approx 100$, as confirmed by the Fig. 6(a). Since values of $\alpha$ in the range $0.1 \le \alpha \le 100$ are acceptable but the performance accuracy decreases systematically in this range, see Fig. 5, we tentatively choose $\alpha = 0.1$ as an optimal value. Finally, one can also see how the oscillatory behavior originating from multicollinearity is suppressed in Ridge regression. This point can already be grasped by comparing Fig. 3(a) and Fig. 3(b), but is best understood in the principal component basis, see the Appendix D 3.

Before closing this section, we note that, in addition to the condition number, we also computed the variance inflation factor (VIF) [48], another popular metric for evaluating the degree of multicollinearity. The VIF acts as a scaling factor for the variance of the estimated weights arising from dataset-to-dataset fluctuations, and it can diverge under severe multicollinearity even when the sample size is sufficiently large. We find that the VIF provides results consistent with the condition number, when applying the usual rules of thumb [48]. However, this quantity does not offer additional insight for our dataset, because the sample size is quite large (and thus the estimated error bars in Fig. 3 are small). For this reason, we do not report these results in this paper and leave their investigation to future work, particularly in situations with limited sample sizes.

# VI. TOWARDS INTERPRETABLE LINEAR REGRESSION MODELS

The analysis presented in Sec. V shows that Ridge regression is an effective method to cope with the effects of multicollinearity in structure-dynamics datasets. However, the solutions of this linear model are still too high-dimensional to be physically interpretable. To achieve a simple physical picture, we must substantially reduce the dimensionality of the problem at hand. In this section, we consider a generalization of Ridge regression, called elastic net [72], that provides means to select the most relevant features, while reducing multicollinearity (Sec. VI A). As an alternative approach, we use principal component regression, which builds upon a simple linear transformation of original features (Sec. VI B).

## A. Elastic net regression

### 1. Definition

The elastic net model is an extension of Ridge regression to perform feature selection [72]. In this method, the regularization term in Eq. (18) reads

<figure><img src="image_8.png" /><figcaption>FIG. 7. Elastic net regression of the dynamic propensity using the BP descriptor. Main panel: Pearson coefficient $R[Y, \hat{Y}]$ for the optimal models (large circles) as a function of the number of selected features $P$. The values of $\alpha$ and $\beta$ corresponding to optimal models with $P \le 10$ are indicated in the panel. The small circles indicate results for values $(\alpha, \beta)$ corresponding to sub-optimal models. Inset: number of selected features $P$ as a function of $\beta$ for selected values of $\alpha$.</figcaption></figure>

$$ \mathcal{L}^{\text{reg}}(\hat{\mathbf{w}}) = a \sum_{f=1}^{M} |\hat{w}^{(f)}| + \frac{b}{2} \sum_{f=1}^{M} (\hat{w}^{(f)})^2, \quad (32) $$

where $a$ and $b$ are two regularization parameters. In addition to the L2 norm term already included in Ridge regression, Eq. (32) now includes an L1 norm regularization term. Like in Lasso regression [70], the L1 term shrinks some of the weights to zero within numerical accuracy to perform feature selection. Following the elastic net implementation of scikit-learn [73], we define $\alpha = a + b$ and $\beta = a/(a+b)$, where now $0 \le \beta \le 1$. $\alpha$ quantifies the overall strength of regularization, incorporating both L1 and L2 terms, while $\beta$ controls the relative contribution of L1 regularization. $\beta = 0$ and $\beta = 1$ correspond to pure Ridge and Lasso regression, respectively. For any pair $(\alpha, \beta)$, the method finds an optimal number of features, $P$, that minimizes the loss function.

### 2. Selection of parsimonious models

To get a feeling of the role of $\beta$, we show in the inset of Fig. 7 the dependence of the number of selected features $P$ on $\beta$, for a few fixed values of $\alpha$. The analysis is performed as usual on the test set. As expected, the number of selected features decreases as the strength of the L1 regularization increases, reaching its minimum for $\beta = 1$. A similar effect is observed with increasing $\alpha$.

Of course, the prediction performance of the model will vary with $\alpha$ and $\beta$. To define a set of optimal models that provide the best prediction accuracy on the test set, while minimizing the number of selected features, we proceed as follows. For
each fixed value of $\alpha$, we perform elastic net regression starting from $\beta = 1$ (Lasso regression). If the algorithm converges to a solution, we have found the best model for that value of $\alpha$. Otherwise, we reduce $\beta$ until a solution is found, and extract the corresponding value of $P$. We then compute the Pearson coefficient $R$ with the dynamic propensity for each optimal model at a given $\alpha$.

The results of this procedure are shown in the main panel of Fig. 7 as large black symbols. We also include the results of a more extensive sampling of $(\alpha, \beta)$ pairs, shown as small dots. We see that the optimal models identified by our procedure always provide the best performance for a given $P$. The envelope of $R(P)$ defines a set of optimal models, in the sense that they provide the best performance accuracy for a given size $P$ of the feature space. The typical values of the Pearson coefficients for the most parsimonious models range from 0.54 for $P = 1$, to 0.62 for $P = 5$, to 0.74 for $P = 10$. These values can be compared with the maximum value of 0.88 achieved when retaining the full descriptor, see also Sec. VII. Within our strategy, we found that pure Lasso regression ($\beta = 1$) yields the optimal model, except in the case $M = 1$.

In Table II, we summarize the features selected by parsimonious models ($P \le 10$) along the envelope of Fig. 7. We include the corresponding values of the Pearson coefficients with the dynamic propensity, both for the regression models and for the individual features. We see that the low-dimensional models of elastic net regression (Lasso regression for $P \ge 2$) successfully pinpoint the features that are most correlated to the dynamic propensity, found in either angular or radial sectors of the BP descriptor. We notice, however, that these low-dimensional models are still somewhat redundant. For instance, the Pearson coefficient between the two features selected by the $P = 2$ model is appreciable, $R = 0.62$. Even worse, the model with $P = 5$ contains two nearly identical angular features (163 and 164), having mutual correlation $R = 0.99$. This can be explained by the fact that for these models $\beta = 1$, hence, in practice, there is no suppression of multicollinearity. It would be possible to address this issue by using a minimum-redundancy-maximum-relevance selection scheme [74], which includes iteratively features that are highly correlated with the target, while penalizing those strongly correlated with the ones already taken. We leave this approach for a future investigation.

Although our numerical results indicate that the method is not able to fully remove multicollinearity, elastic net (Lasso) regression does overall a good job at selecting parsimonious models, achieving a reasonable correlation with the dynamic propensity ($R \approx 0.6$) using just a couple of features. The best performing low-dimensional model ($R \approx 0.74$), mixing both correlated and anti-correlated features, is obtained with $P = 10$.

## B. Principal component regression

### 1. Definition

In this section, we use principal component analysis (PCA) as an alternative approach to reducing the dimensionality of the problem. This method allows one to find a new set of

relevant features, defined in an orthogonal basis formed by the eigenvectors of the correlation matrix $C$. The corresponding modes are called principal components and can be used to perform linear regression. When only a subset of relevant principal components is selected, this regression model is called principal component regression (PCR).

Let us briefly remind the key ingredients of PCA. Since $C$ is a symmetric matrix, it can be diagonalized using the orthogonal matrix $U = [\mathbf{u}^{(1)}, \mathbf{u}^{(2)}, \dots, \mathbf{u}^{(M)}]$, where $\mathbf{u}^{(f)}$ are the eigenvectors of $C$, and $\lambda^{(f)}$ are the corresponding eigenvalues. This gives the equation

$$C = U\Lambda U^T, \quad (33)$$

where $\Lambda$ is a diagonal matrix with eigenvalues $\lambda^{(1)}, \dots, \lambda^{(M)}$, sorted as $\lambda^{(1)} \ge \dots \ge \lambda^{(M)} > 0$.

We now obtain new features $X'$ through a linear transformation, $X' = XU$. The original feature $\mathbf{X}_i$ in Eq. (9) is expanded using the orthogonal basis $\mathbf{u}^{(1)}, \mathbf{u}^{(2)}, \dots, \mathbf{u}^{(M)}$, and the coefficients $X_i^{(f)}$ represent the new features on this basis. Specifically, we have

$$\mathbf{X}_i = \sum_{f=1}^{M} X_i^{(f)} \mathbf{u}^{(f)}, \quad (34)$$

$$X_i^{(f)} = (\mathbf{u}^{(f)})^T \mathbf{X}_i. \quad (35)$$

From Eq. (35), we see that the mean value of $X_i^{(f)}$ is zero. The covariance matrix for $X'$ is given by

$$\frac{1}{N_S} X'{}^T X' = U^T C U = \Lambda,$$

hence the variance of $X_i^{(f)}$ is $\lambda^{(f)}$. We show $\lambda^{(f)}$ in Fig. 8(a). The first two PCs have very large eigenvalues, while the remaining components form a tail that is difficult to appreciate on a linear scale.

Let us now normalize the features

$$\tilde{X}_i^{(f)} = \frac{X_i^{(f)}}{\sqrt{\lambda^{(f)}}}, \quad (36)$$

so that the new feature $\tilde{X}_i^{(f)}$ has zero mean and unit variance, and determine how $R[\tilde{X}, Y]$ is transformed in the PC basis. Then, by applying $U^T$, we obtain

$$U^T R[\tilde{X}, Y] = \frac{1}{N_S} X'{}^T Y = \Sigma R[\tilde{X}, Y], \quad (37)$$

where $\Sigma$ (not to be confused with summation) is a diagonal matrix whose elements are $\sqrt{\lambda^{(1)}}, \sqrt{\lambda^{(2)}}, \dots, \sqrt{\lambda^{(M)}}$, meaning that $\Sigma^2 = \Lambda$. We show $R[\tilde{X}, Y]$ in Fig. 8(b). Only the first few PCs have non-negligible correlations with the dynamic propensity, say $|R| > 0.3$. This is best appreciated from panels (c) and (d) of Fig. 8, where we show the absolute value of $R[\tilde{X}, Y]$ as a function of $\lambda^{(f)}$. Surprisingly, PC1 has a negligible correlation with the dynamic propensity, while it is PC2 that has the largest correlation. Thus, selecting the PCs on
interpretation for test data predictions as well.

Let us first focus on the OLS regression case. As shown in the Appendix D, the dynamic propensity predicted by OLS regression is

$$ \hat{\mathbf{Y}} = \sum_{f=1}^{M} R[\tilde{\mathbf{X}}^{(f)}, Y] \tilde{\mathbf{X}}^{(f)}. \quad (38) $$

The coefficient in front of each feature $\tilde{\mathbf{X}}^{(f)}$, which reflects its importance in the prediction, is simply the Pearson correlation coefficient shown in Fig. 8(b). This is, of course, because all the PC eigenvectors are orthogonal to each other. The Pearson coefficient between the ground truth propensity and the predicted one is then easily found:

$$ R[Y, \hat{Y}] = \sqrt{\sum_{f=1}^{M} (R[\tilde{\mathbf{X}}^{(f)}, Y])^2}. \quad (39) $$

This expression provides a clear geometric interpretation of how each $R[\tilde{\mathbf{X}}^{(f)}, Y]$ contributes to the performance $R[Y, \hat{Y}]$ in training set.

### 4. Dimensional reduction and comparison with elastic net

To perform dimensional reduction, one selects only $\tilde{\mathbf{X}}^{(1)}, \tilde{\mathbf{X}}^{(2)}, \dots, \tilde{\mathbf{X}}^{(P)}$ with $P \ll M$, yielding

$$ \hat{\mathbf{Y}} = \sum_{f=1}^{P} R[\tilde{\mathbf{X}}^{(f)}, Y] \tilde{\mathbf{X}}^{(f)}. \quad (40) $$

This operation corresponds to neglecting features associated with smaller eigenvalues $\lambda^{(P+1)}, \dots, \lambda^{(M)}$. In Fig. 9, we analyze the performance of PCR by showing $R[Y, \hat{Y}] = \sqrt{\sum_{f=1}^{P} (R[\tilde{\mathbf{X}}^{(f)}, Y])^2}$ as a function of the number of components $P$. Inspired by Eq. (39), we also include results obtained by sorting the features according to the absolute value of $R[\tilde{\mathbf{X}}^{(f)}, Y]$. This corresponds to a supervised feature extraction scheme. On the side of the figure, we include as ticks the Pearson correlation coefficients achieved with a given number of features, $P$. The visual impression is that the first few PCs contribute to the bulk of the correlation, reaching $R = 0.7$ with $P = 5$. The rest of the PCs constitute a dense, poorly interpretable background of orthogonal variables that are essentially uncorrelated with the propensity, see Fig. 8(b). Linear regression methods can nonetheless harvest this background to reconstruct the dynamic propensity field precisely.

It is interesting to compare these results to those obtained by the optimal elastic net models, which are included as dashed lines in Fig. 9. We see that the elastic net models converge to the maximum correlation, as a function of selected features, approximately like in supervised PCR. For small numbers of features, the two regression models provide about the same performance accuracy, although the precise details depend on the descriptor. One advantage of PCR is that the feature extraction scheme is straightforward and efficient, thanks to

<figure><img src="image_9.png" /><figcaption>FIG. 9. Pearson coefficient $R$ between the ground truth dynamic propensity and PCR predictions using the first $P$ PCs sorted according to their eigenvalue (squares) or their Pearson coefficient with the dynamic propensity (circles). The results of the optimal elastic net regression models are also included (triangles). Note the change of scale on the x-axis after $P = 5$ components. The ticks on the right side of the figure indicate the values of the Pearson coefficient reached for a given value of number of PCs included in the regression model.</figcaption></figure>

<table><caption>TABLE IV. Pearson coefficient between the dynamic propensity and selected linear regression models for all the investigated descriptors. The performance of the JBB descriptor without regularization ($\alpha = 0$) is very low, due to severe overfitting.</caption><thead><tr><th rowspan="2"></th><th rowspan="2">Features</th><th colspan="2">Ridge</th><th colspan="2">PCR</th></tr><tr><th>α = 0</th><th>α = 0.1</th><th>P = 5</th><th>P = 2</th></tr></thead><tbody><tr><td>BP descriptor</td><td>276</td><td>0.86</td><td>0.85</td><td>0.69</td><td>0.55</td></tr><tr><td>JBB descriptor</td><td>120</td><td>0.25</td><td>0.87</td><td>0.80</td><td>0.53</td></tr><tr><td>SLO descriptor</td><td>60</td><td>0.86</td><td>0.85</td><td>0.84</td><td>0.81</td></tr></tbody></table>

the orthogonality of the PCA basis. By contrast, although the features selected by elastic net are not orthogonal to each other, they are more easily interpretable than the PC eigenvectors. The latter, in fact, mix all the features of the original descriptor and do not always lean themselves to a straightforward interpretation. We will further discuss this issue in Sec. VII.

## VII. DISCUSSION

Given the range of linear models that reproduce the dynamic propensity with comparable performance accuracy, a few questions arise naturally: is it possible to identify an “optimal” model that strikes a balance between prediction accuracy and physical interpretation? Does the choice of the descriptor play an important role? In this section, we will provide elements to try to address these questions.

**Real space structure of the models' prediction**– A first piece of information comes from direct visual comparison of the ground truth dynamic propensities and the predictions of the linear regression models. Here, we consider two low-dimensional PCR models, namely the ones with $P = 2$ and
<figure><img src="image_10.png" /><figcaption>FIG. 10. Dynamic propensity of a representative configuration: (a) ground truth and estimates from (b) the OLS regression model, (c) the PCR model with $P = 5$ components, and (d) the PCR model with $P = 2$ components. Dark and light particles correspond to fast and slow particles, respectively. The mid panels, from (e) to (h) show the same as in the top panels coarse-grained over a length $\ell = 1.5$. Particles are colored in black and white according to whether the corresponding coarse-grained variable is above or below the average, respectively. (i) Probability density function of the ground truth and estimated dynamic propensities for various models, using all the available configurations. The remaining panels in the bottom row display scatter plots of the ground truth propensity against its estimated value for (j) the OLS regression model, (k) the PCR model with $P = 5$ components, and (l) the PCR model with $P = 2$ components. In all the panels, $\hat{Y}$ is shifted and scaled to match the mean and variance of the ground truth $Y$.</figcaption></figure>

$P = 5$ within supervised PCR, as well as with the full $P = M$ model (corresponding to the OLS regression model) obtained with the BP descriptor. In this section, we revert our initial feature normalization (see Sec. III E): we scale and shift $\hat{Y}_i$ obtained from the PCR models so that their mean and variance match those of the actual unscaled dynamic propensity. Thus, we will consider $\hat{Y}_i$ vs $Y_i$ in their original LJ units.

Figure 10 shows the distributions of $\hat{Y}_i$ vs $Y_i$, the scatter plots of $\hat{Y}_i$ vs $Y_i$ and the corresponding snapshots for a representative configuration in our dataset. We see that the full descriptor ($P = M$ or OLS regression model) provides an excellent description of the dynamic propensity field, which is accurately reconstructed in most of its details, except for a slight discrepancy in the shape of the distribution, see Fig. 10(i). Thus, $R$ values of the order of 0.9 correspond to a very satisfactory description of the dynamic propensity. The model obtained with the $P = 5$ most correlated PCs gives $R \approx 0.7$ and reproduces most of the patterns of the dynamic propensity field, despite some additional noise. The correspondence is good when considering coarse-grained fields (see the mid panels), indicating that such models are able to grasp the relevant dynamic fluctuations on large length scales [75]. The model

with the $P = 2$ most correlated PCs, instead, only captures some of the fast and slow regions of the propensity, but the overall large-scale structure of the propensity field is not accurately reproduced. Qualitatively, these results suggest that correlation coefficients of 0.55, 0.7, 0.9 correspond to modest, satisfactory, and excellent estimation of the dynamic propensity, respectively.

*Dependence on the descriptor*— How sensitive are the results to the choice of the descriptor? From Table IV we see that the prediction performance of the optimal Ridge regression model is rather insensitive to it. Note that this similarity is partly due to using the inherent structure configurations: when the calculations are performed on the instantaneous configurations, the prediction performance improves slightly with increasing the dimensionality of the descriptor (not shown).

Even when using inherent structures, however, a clear difference can be seen at the level of parsimonious models. Interestingly, the low-dimensional models obtained from PCR perform better with the SLO descriptor than with the BP and JBB descriptors. This might be explained by the inclusion from the outset of some physically relevant features like the $\Theta$ parameter. With the SLO descriptor, in fact, even a two-features model
<figure><img src="image_11.png" /><figcaption>FIG. 11. Same as Fig. 9 but for the SLO descriptor.</figcaption></figure>

achieves a high Pearson coefficient ($R \approx 0.81$) and very little is gained by including the full set of features ($R \approx 0.86$). By contrast, high-dimensional descriptors accurately reconstruct the observable of interest, exploiting by brute force a broad spectrum of features, but they may fail to provide an interpretable parsimonious models. This could also partly depend on the details of the descriptor: the atomic cluster expansion, for instance, successfully predicts plasticity in amorphous solids even upon strong dimensional reduction [76].

**Interpretation of the structural modes**– The above findings call for an inspection of the dominant modes in PCR using a physically motivated descriptor. In Fig. 11, we show the PCR results for the SLO descriptor. Again, the PCs that are most correlated with the dynamic propensity are not necessarily the ones with the largest eigenvalues. We found that the modes with the largest correlations are PC2 and PC5, whose Pearson coefficients with the dynamic propensity sum up geometrically to 0.81.

In Fig. 12, we show the eigenvectors corresponding to PC1, PC2, and PC5. Let us first focus on the PC that is most relevant for the dynamics. PC2 account for fluctuations of $\bar{\Theta}$ that are anti-correlated to the local packing $\bar{\varphi}$ (and correlated to $\bar{\rho}$, as expected). This makes sense, as small values of $\bar{\Theta}$ should capture sterically favored environments with high local packing. Note that the coarse-graining length associated with the larger contribution on this PC is around $\ell \approx 2.5$, which roughly corresponds to the second coordination shell. We thus identify PC2 with a structural mode associated with the fluctuations of local packing on an intermediate length scale.

Interestingly, PCA reveals that some of the fluctuations of $\bar{\Theta}$ are positively correlated to $\bar{\varphi}$. This unexpected behavior is described by PC1, which captures the largest fraction of the variance of the normalized structural dataset. Our supervised PCR scheme effectively removes this source of fluctuations, which is irrelevant for the dynamics ($R = 0.18$). Note that the fluctuations described by PC2 are more strongly correlated to the dynamic propensity than the bare $\bar{\Theta}$ parameter.

Turning our attention to PC5, we see that the largest contribution to this mode comes from the 6-fold orientational order parameter, $\bar{\Psi}_6$, averaged over a short range distance ($\ell \approx 1.0$).

<figure><img src="image_12.png" /><figcaption>FIG. 12. Eigenvectors $\mathbf{u}^{(1)}, \mathbf{u}^{(2)},$ and $\mathbf{u}^{(5)}$ obtained using the SLO descriptor.</figcaption></figure>

We point out, however, that the relevant bond-orientational order of our ternary model need not be 6-fold. Indeed, previous analysis of the low-energy states of the binary version of the model indicated the presence of quasi-crystalline order with 5- and 10-fold symmetries. We leave a more systematic analysis of the preferred bond-orientational order of our ternary model for a future study.

To sum up on this point, PCR of the SLO descriptor yields a simple two-state structural model. The two relevant variables account for steric effects on an intermediate length scale and for some of the fluctuations of bond-orientational order. Superficially, these results suggest a connection with the phenomenological two-state model developed by Tanaka in the 1990s [53], which included local density and bond-order fluctuations as structural order parameters. We emphasize, however, that our ternary model does not display strong 6-fold orientational order and $\bar{\Psi}_6$ itself is poorly correlated to the dynamic propensity. It is only when the fluctuations of $\bar{\Psi}_6$ are coupled to those of the other features as described by $\mathbf{u}^{(5)}$ that some correlation emerges. We found that the other descriptors consistently identify a structural mode related to local packing (e.g., PC2 in the BP descriptor, PC1 in the JBB descriptor) that is moderately correlated to the dynamics, but none that is specifically connected to bond-orientational order.

**Optimality criteria**– Are there quantitative criteria to decide whether to accept a low-dimensional model, which is easier to interpret, at the expense of performance? We considered several optimality criteria, developed in the context of statistical learning, that attempt to define an optimum model by balancing the performance and the model complexity. In particular, we
employed two well-known scores: the Akaike Information Criterion (AIC) [77] and the Bayesian Information Criterion (BIC) [78]. Both criteria address over-fitting by combining the maximized log-likelihood, which quantifies data fidelity, with a complexity penalty term that grows with the number of parameters. In general, AIC applies a lighter penalty and therefore tends to keep more features if they help prediction, while BIC is stricter and usually picks the smaller model that still describes the data well. We evaluated AIC and BIC for our glass-forming liquids dataset using Lasso regression with varying magnitudes of the L1 penalty, which controls the number of selected features. Both AIC and BIC indicated that the optimal model corresponds to OLS (i.e., including all features). Thus, in our setting, AIC and BIC do not provide any practical guidance.

Our current standpoint on this issue is that there is at present no robust quantitative criterion to select an optimal model that strikes a balance between prediction and model size. We think that the simple empirical criterion of selecting a small number of features, say up to 5, is in line with a well-established approach when building phenomenological models in statistical physics [51]. We suggest that looking at the best two- or five-feature model provides a basic “handful criterion” to identify parsimonious data-driven models of glassy dynamics.

## VIII. CONCLUSIONS

Recent work has shown that data-driven models, based on either deep neural networks or linear regression, can accurately describe and predict the dynamic propensity of glass-forming liquids on the structural relaxation time scale [24]. In our opinion, the current level of prediction accuracy of these models is high enough to motivate a shift of focus toward interpretation [33]. In fact, a common criticism to these machine learning studies is that they still provide little physical insight into the underlying mechanisms behind glassy dynamics. Identifying interpretable models that provide a robust and succinct relationship between physically relevant variables is crucial to address this issue.

The main goal of this paper was to assess and improve the interpretability of linear regression models of glassy dynamics using a simple model of two-dimensional glass. We showed that, contrary to previous expectations [24], even linear regression models can be hard to interpret. A major issue can be traced back to the presence of strong linear dependencies between structural features, known as multicollinearity, which hinders the interpretability of linear models. We found that several structural descriptors used in recent studies are severely affected by multicollinearity. Ridge regression can be used to suppress some of the detrimental effects of multicollinearity, but the resulting models are not succinct enough to be interpretable.

To identify low-dimensional linear models of glassy dynamics, we used two simple dimensional reduction techniques. First, we considered elastic net regression, which yields a set of parsimonious models that strike a good balance between accuracy and physical interpretability. Second, we performed principal component regression using a supervised selection

scheme of the principal components. This approach yields an accurate enough description of the dynamic propensity with a handful of collective structural features. Overall, our work establishes that linear regression models, once properly fine-tuned, can be useful tools to identify the relevant structural modes associated to dynamic heterogeneities in glass-forming liquids. It would be interesting to extend our analysis to more sophisticated high-dimensional descriptors and to other glassy systems beyond the simple two-dimensional model liquid considered here, including amorphous solids subjected to external loading.

## ACKNOWLEDGMENTS

We thank Gerhard Jung for insightful discussions. MO thanks the support by MIAI@Grenoble Alpes and the Agence Nationale de la Recherche under France 2030 with the reference ANR-23-IACL-0006.

## DATA AVAILABILITY

The data and workflow necessary to reproduce the findings of this study will be openly available after publication of the paper in the Zenodo data repository.

## Appendix A: Selected results for the SLO and JBB descriptors

In this section, we present additional results for the physically motivated descriptors introduced in Sec. III D. Figures 13 and 14 display the Pearson coefficients $R[X^{(f)}, Y]$ between the dynamic propensity and each structural feature $X^{(f)}$ of the SLO and JBB descriptors, respectively.

The SLO descriptor emphasizes the substantial role of packing efficiency in determining dynamic fluctuations. This can be appreciated by the large positive correlations between $\bar{\Theta}$ and the dynamic propensity, as well as by the negative correlations with $\bar{\phi}$. Note, however, that the coarse-graining length $\ell$ that gives the largest correlation is not the same for $\bar{\Theta}$ and $\bar{\phi}$. The bond-orientational order parameter, $\bar{\Psi}_6$, by contrast, is poorly correlated with the dynamics, see the discussion in Sec. VI B.

It is instead more difficult to extract straightforward physical information from the analysis of the JBB descriptor, due to its complex dependence on the chemical composition. It is nonetheless interesting to note that the perimeter $\bar{\pi}$ of the Voronoi cells, which is included in the JBB descriptor but not in the SLO one, has an overall negative correlation with the dynamic propensity. This again points to a coupling between steric effects, i.e., larger local free volume, and mobility.

In Figures 15 and 16, we show the correlation matrices C of the SLO and JBB descriptors, respectively. As for the BP descriptor, the correlation matrices display a pronounced block structure: there are strong correlations within subblocks of features, this time due to similar coarse-graining lengths, as well as non-trivial correlations between different features. In
<figure><img src="image_13.png" /><figcaption>FIG. 13. Pearson coefficient, $R[X^{(f)}, Y]$, between the dynamic propensity $Y$ and each structural feature $X^{(f)}$ of the SLO descriptor for $f = 1, \dots, M$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIG. 14. Pearson coefficient, $R[X^{(f)}, Y]$, between the dynamic propensity $Y$ and each structural feature $X^{(f)}$ of the JBB descriptor for $f = 1, \dots, M$.</figcaption></figure>

In particular, the SLO descriptor reveals a negative correlation between local number density $\bar{\rho}$ and local packing fraction $\bar{\phi}$. This counter-intuitive result is likely a result of compositional fluctuations in this mixture, see Sec. III F.

## Appendix B: Ordinary least squares regression and Ridge regression

We derive the estimation of the weights $\hat{w}$ by the ordinary least squares regression (OLS) and Ridge regression for the dataset, $(X, Y)$, defined in Eqs. (8) and (10).

We first consider a linear model,

$$ \hat{\mathbf{Y}} = X\hat{\mathbf{w}}, \qquad (B1) $$

where $\hat{\mathbf{Y}} = [\hat{Y}_1, \hat{Y}_2, \dots, \hat{Y}_{N_S}]^T$ represents the predictions and $\hat{\mathbf{w}} = [\hat{w}^{(1)}, \hat{w}^{(2)}, \dots, \hat{w}^{(M)}]^T$ denotes the weight parameters.

The loss function $\mathcal{L}^{\text{Ridge}}(\hat{\mathbf{w}})$ in Eq. (23) can be rewritten as

$$ \mathcal{L}^{\text{Ridge}}(\hat{\mathbf{w}}) = \frac{1}{2N_S} \left( \hat{\mathbf{w}}^T X^T X \hat{\mathbf{w}} - 2 (X^T Y)^T \hat{\mathbf{w}} + Y^T Y \right) + \frac{\alpha}{2} \hat{\mathbf{w}}^T \hat{\mathbf{w}}. \qquad (B2) $$

<figure><img src="image_15.png" /><figcaption>FIG. 15. Correlation matrix C for the SLO descriptor.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>FIG. 16. Correlation matrix C for the JBB descriptor.</figcaption></figure>

The derivative of $\mathcal{L}^{\text{Ridge}}(\hat{\mathbf{w}})$ with respect to $\hat{\mathbf{w}}$ is

$$ \nabla_{\hat{\mathbf{w}}} \mathcal{L}^{\text{Ridge}}(\hat{\mathbf{w}}) = \left( \frac{1}{N_S} X^T X + \alpha I \right) \hat{\mathbf{w}} - \frac{1}{N_S} X^T Y, \qquad (B3) $$

where $I$ is the $M \times M$ identity matrix.

Since $Y_i$ (elements of $Y$) and $X_i^{(f)}$ (elements of $X$) are normalized to have zero mean and unit variance, we can express the terms in Eq. (B3) using the correlation matrix C
defined in Eq. (15) and the Pearson coefficients $R$,

$$ \frac{1}{N_S} X^T X = C, \qquad (B4) $$

$$ \frac{1}{N_S} X^T Y = \begin{bmatrix} R[X^{(1)}, Y], \dots, R[X^{(M)}, Y] \end{bmatrix}^T \\ = \mathbf{R}[X, Y]. \qquad (B5) $$

Setting $\nabla_{\hat{w}} \mathcal{L}^{\text{Ridge}}(\hat{w}) = 0$, the solution for Ridge regression is obtained as

$$ \hat{w}_{\text{Ridge}} = (C + \alpha I)^{-1} \mathbf{R}[X, Y]. \qquad (B6) $$

## Appendix C: Condition number of a matrix

We briefly review the condition number of a matrix, without going into the details of a mathematically rigorous treatment. It is defined using the norm of a matrix, and we also review its meaning as a measure of a matrix's instability. One can see detailed discussions in, e.g., Refs. 79 and 80.

### 1. Matrix norm

A matrix norm is a generalization of the absolute value $| \cdot |$ for scalars and the vector norm $|| \cdot ||$ for vectors, applied to matrices. The matrix norm is also denoted as $|| \cdot ||$. Here, we summarize only the important properties relevant for our purposes. For a scalar $\alpha$, a vector $\mathbf{x}$, and matrices A and B, the following properties hold: i) $||A|| \ge 0$, ii) $||\alpha A|| = |\alpha|||A||$, iii) $||A+B|| \le ||A||+||B||$, iv) $||A\mathbf{x}|| \le ||A|| ||\mathbf{x}||$, v) $||AB|| \le ||A|| ||B||$, and vi) $||I|| = 1$, where I is the identity matrix.

There are various ways to define the matrix norm that satisfy these properties. One convenient approach is to use the maximum singular value, $s_{\max}$, of a matrix. The norm $||A|| = s_{\max}$ is referred to as the spectral norm.

For a positive semi-definite symmetric matrix, the singular values are equal to the eigenvalues. Thus, we have $||A|| = \lambda_{\max}$, where $\lambda_{\max}$ is the largest eigenvalue. Moreover, if A is invertible, then $||A^{-1}|| = \lambda_{\min}^{-1}$, where $\lambda_{\min}$ is the smallest eigenvalue.

### 2. Definition of condition number

In general, the condition number $\kappa(A)$ is defined as

$$ \kappa(A) = \frac{||A||}{||A^{-1}||}. \qquad (C1) $$

In this paper, we consider the spectral norm, and assume that A is a positive semi-definite symmetric matrix. Thus, we obtain

$$ \kappa(A) = \frac{\lambda_{\max}}{\lambda_{\min}}. \qquad (C2) $$

### 3. Instability of matrix

The condition number $\kappa(A)$ corresponds to a degree of instability of a matrix A. This concept can be best understood through perturbation analysis of a linear system.

Suppose we wish to solve the linear system

$$ \mathbf{A}\mathbf{x} = \mathbf{b}, \qquad (C3) $$

where A is a square matrix, and $\mathbf{x}$ and $\mathbf{b}$ are vectors. When A is invertible, the solution is given by

$$ \mathbf{x} = \mathbf{A}^{-1}\mathbf{b}. \qquad (C4) $$

In practice, measuring A and b involves errors $\delta A$ and $\delta b$ due to, for example, lack of statistics, numerical errors, etc. We now ask how the error in $\mathbf{x}$, denoted as $\delta \mathbf{x}$, is induced by errors in A and/or b.

1) When b has error $\delta b$, the linear equation becomes

$$ A(\mathbf{x} + \delta \mathbf{x}) = \mathbf{b} + \delta \mathbf{b}, \qquad (C5) $$

which, using Eq. (C3), leads to $A\delta \mathbf{x} = \delta \mathbf{b}$, so that

$$ \delta \mathbf{x} = \mathbf{A}^{-1} \delta \mathbf{b}. \qquad (C6) $$

Thus, the relative error $||\delta \mathbf{x}||/||\mathbf{x}||$ is bounded as follows:

$$ \frac{||\delta \mathbf{x}||}{||\mathbf{x}||} = \frac{||\mathbf{A}^{-1} \delta \mathbf{b}||}{||\mathbf{x}||} \le \frac{||\mathbf{A}^{-1}|| ||\delta \mathbf{b}||}{||\mathbf{x}||}. \qquad (C7) $$

Additionally, since $||\mathbf{b}|| = ||\mathbf{A}\mathbf{x}|| \le ||\mathbf{A}|| ||\mathbf{x}||$, we obtain the bound:

$$ \frac{1}{||\mathbf{x}||} \le \frac{||\mathbf{A}||}{||\mathbf{b}||}. \qquad (C8) $$

Thus, Eq. (C7) becomes

$$ \frac{||\delta \mathbf{x}||}{||\mathbf{x}||} \le \frac{||\mathbf{A}|| ||\mathbf{A}^{-1}|| ||\delta \mathbf{b}||}{||\mathbf{b}||}. \qquad (C9) $$

Using the definition of the condition number $\kappa(A)$ in Eq. (C1), we conclude that

$$ \frac{\frac{||\delta \mathbf{x}||}{||\mathbf{x}||}}{\frac{||\delta \mathbf{b}||}{||\mathbf{b}||}} \le \kappa(A). \qquad (C10) $$

The ratio between the relative error $||\delta \mathbf{x}||/||\mathbf{x}||$ and $||\delta \mathbf{b}||/||\mathbf{b}||$ quantifies how tiny perturbations or fluctuations in b influence the error in $\mathbf{x}$, which can be interpreted as a form of susceptibility. Equation (C10) shows that this ratio is bounded by the condition number $\kappa(A)$. When $\kappa(A)$ is small, the linear system in Eq. (C3) is stable. However, when $\kappa(A)$ is large, it becomes unstable.

2) When A has error $\delta A$, one can similarly derive the inequality for the relative error $||\delta \mathbf{x}||/||\mathbf{x}+\delta \mathbf{x}||$ when A has an error $\delta A$:

$$ \frac{\frac{||\delta \mathbf{x}||}{||\mathbf{x}+\delta \mathbf{x}||}}{\frac{||\delta A||}{||\mathbf{A}||}} \le \kappa(A). \qquad (C11) $$
3) When both **b** and A have errors, δ**b** and δA, one can also
derive the inequality when both **b** and A have errors, δ**b** and
δA:

$$
\frac{\frac{\|\delta\mathbf{x}\|}{\|\mathbf{x}\|}}{\frac{\|\delta\mathbf{A}\|}{\|\mathbf{A}\|} + \frac{\|\delta\mathbf{b}\|}{\|\mathbf{b}\|}} \le \frac{\kappa(\mathbf{A})}{1 - \|\mathbf{A}^{-1}\delta\mathbf{A}\|}. \quad (C12)
$$

We assumed $||A^{-1}\delta A|| < 1$.

In summary, the condition number provides an upper bound for the relative error in the linear system. When the condition number $\kappa(A)$ is small, tiny perturbations $\delta A$ and/or $\delta b$ lead to only small relative errors in $\mathbf{x}$, and the linear system is stable. When the condition number is large, however, even tiny perturbations may induce huge errors, making the linear system unstable, even if A is invertible.

**Appendix D: Ridge regression in the PCA basis**

1. Expression of the Pearson coefficient

Let us consider the Ridge regression solution, Eq. (24), in the
PCA setting. Using Eq. (37) and $(C+\alpha I)^{-1} = U(\Lambda+\alpha I)^{-1}U^T$,
we get

$$
\begin{align*}
\hat{\mathbf{Y}} &= X\hat{\mathbf{w}}_{\text{Ridge}} = X(C + \alpha I)^{-1}\mathbf{R}[\mathbf{X}, \mathbf{Y}] \\
&= X'(\Lambda + \alpha I)^{-1}\Sigma\mathbf{R}[\tilde{\mathbf{X}}, \mathbf{Y}] \\
&= \sum_{f=1}^{M} R[\tilde{\mathbf{X}}^{(f)}, \mathbf{Y}] \left( \frac{\lambda^{(f)}}{\lambda^{(f)} + \alpha} \right) \tilde{\mathbf{X}}^{(f)}. \tag{D1}
\end{align*}
$$

Using the PCA basis, we obtain the Pearson coefficient between the ground truth $Y_i$ and the prediction $\hat{Y}_i$

$$
\begin{align*}
R[Y, \hat{Y}] &= \frac{1}{N_S \sqrt{\mathrm{Var}[\hat{Y}]}} \sum_{i \in S} Y_i \hat{Y}_i = \frac{1}{N_S \sqrt{\mathrm{Var}[\hat{Y}]}} \hat{Y}^T Y \\
&= \frac{1}{\sqrt{\mathrm{Var}[\hat{Y}]}} \hat{Y}^T \mathbf{R}[\mathbf{X}, Y] \tag{D2}
\end{align*}
$$

and

$$
\mathrm{Var}[\hat{Y}] = \frac{1}{N_S} \hat{\mathbf{Y}}^T \hat{\mathbf{Y}} = \hat{\mathbf{w}}^T C \hat{\mathbf{w}}, \quad (\text{D3})
$$

where we used Eqs. (17), (21), and (22).

By plugging in the Ridge regression weights, $\hat{\mathbf{w}}_{\text{Ridge}}$, from
Eq. (24), we get

$$
R[Y, \hat{Y}] = \frac{\mathbf{R}[\mathbf{X}, Y]^T (C + \alpha I)^{-1} \mathbf{R}[\mathbf{X}, Y]}{\sqrt{\mathbf{R}[\mathbf{X}, Y]^T (C + \alpha I)^{-1} C (C + \alpha I)^{-1} \mathbf{R}[\mathbf{X}, Y]}} \quad (\text{D4})
$$

Using the PCA basis, namely Eqs. (33) and (37), we finally
arrive at Eq. (D5) in the main text.

2. Role of the regularization parameter

Looking at the final expression in Eq. (D1), we see that,
or a given $\alpha > 0$, when $\lambda^{(f)} \gg \alpha$, the factor involving the
eigenvalues becomes $\lambda^{(f)}/(\lambda^{(f)} + \alpha) \approx 1$. On the other
hand, when $\lambda^{(f)} \ll \alpha$, $\lambda^{(f)}/(\lambda^{(f)} + \alpha) \approx 0$, hence such
features are suppressed in the sum. Therefore, through the
factor $\lambda^{(f)}/(\lambda^{(f)} + \alpha)$, $\alpha$ acts as a soft threshold, determining
whether we include a feature or not.

It is then easy to find the Pearson coefficient, see the Appendix D,

$$
R[Y, \hat{Y}] = \frac{\sum_{f=1}^{M} (R[\tilde{X}^{(f)}, Y])^2 \left( \frac{\lambda^{(f)}}{\lambda^{(f)} + \alpha} \right)}{\sqrt{\sum_{f=1}^{M} (R[\tilde{X}^{(f)}, Y])^2 \left( \frac{\lambda^{(f)}}{\lambda^{(f)} + \alpha} \right)^2}}. \quad (D5)
$$

We observe that the same threshold factor appears again in
the expression, offering a clear interpretation of the role of
Ridge regularization: in the PCA basis, the regularization term
softly suppresses modes associated with eigenvalues smaller
than $\alpha$. Thus, smaller eigenvalues, potentially associated to
multicollinearity, are systematically and softly suppressed as $\alpha$
increases.

It turns out, however, that the Ridge regularization does not
provide a practical benefit within the PCR scheme. We found,
in fact, that the order of the features is unaltered for small *P*,
both when sorting against *λ* and against *R*. Hence, *α* does not
affect feature selection discussed in Sec. VIB. Moreover, the
maximum correlation achieved when all the *M* features are
included decreases with increasing *α*.

3. Suppression of the oscillatory behavior of the weights

By considering the Ridge regression problem in the PCA
basis, we can now provide a simple theoretical argument to
illustrate how regularization suppresses the oscillatory behavior
of the weights.

We start from Eq. (24) and use the eigenvalue decomposition
of $(C + \alpha I)^{-1}$, which allows us to express $\hat{\mathbf{w}}_{\text{Ridge}}$ as a sum of
projections of the Pearson coefficient vector $\mathbf{R}[\mathbf{X}, Y]$ onto each
eigenmode:

$$
\hat{\mathbf{w}}_{\text{Ridge}} = \sum_{f=1}^{M} (\lambda^{(f)} + \alpha)^{-1} \mathbf{u}^{(f)} \mathbf{u}^{(f)T} \mathbf{R}[\mathbf{X}, \mathbf{Y}]. \quad (\text{D6})
$$

As discussed in Sec. V A, the oscillatory behavior can be
highlighted in the two-feature model, namely,

$$
\hat{\mathbf{w}}_{\text{Ridge}} = \left[ \frac{\mathbf{u}^{(1)}\mathbf{u}^{(1)T}}{\lambda_{\text{max}} + \alpha} + \frac{\mathbf{u}^{(2)}\mathbf{u}^{(2)T}}{\lambda_{\text{min}} + \alpha} \right] \mathbf{R}[\mathbf{X}, \mathbf{Y}], \quad (\text{D7})
$$

where $\mathbf{R}[\mathbf{X}, \mathbf{Y}] = [R^{(1)}, R^{(2)}]^T$ and $\hat{\mathbf{w}}_{\text{Ridge}} = [\hat{w}_{\text{Ridge}}^{(1)}, \hat{w}_{\text{Ridge}}^{(2)}]^T$.
This is further rewritten as

$$
\hat{w}_{\text{Ridge}}^{(1)} = \frac{R^{(1)} + R^{(2)}}{2(\lambda_{\text{max}} + \alpha)} + \frac{R^{(1)} - R^{(2)}}{2(\lambda_{\text{min}} + \alpha)}, \quad (\text{D8})
$$

$$
\hat{w}_{\text{Ridge}}^{(2)} = \frac{R^{(1)} + R^{(2)}}{2(\lambda_{\text{max}} + \alpha)} - \frac{R^{(1)} - R^{(2)}}{2(\lambda_{\text{min}} + \alpha)}. \quad (\text{D9})
$$
When $\alpha \to 0$, $\hat{w}_{\text{Ridge}} \to \hat{w}_{\text{OLS}}$, and one observes the oscillatory behavior, $\hat{w}_{\text{OLS}}^{(1)} \approx -\hat{w}_{\text{OLS}}^{(2)}$, with large magnitude, arising from the mode associated with $\lambda_{\min} \to 0$. For sufficiently large

$\alpha$, the contribution from such mode is suppressed, thereby mitigating the oscillatory behavior.

[1] L. Berthier and G. Biroli, Theoretical perspective on the glass transition and amorphous materials, *Reviews of modern physics* **83**, 587 (2011).

[2] K. Binder and W. Kob, *Glassy materials and disordered solids: An introduction to their statistical mechanics* (World Scientific, 2011).

[3] M. D. Ediger, C. A. Angell, and S. R. Nagel, Supercooled liquids and glasses, *J. Phys. Chem.* **100**, 13200 (1996).

[4] M. D. Ediger, Spatially heterogeneous dynamics in supercooled liquids, *Annu. Rev. Phys. Chem.* **51**, 99 (2000).

[5] L. Berthier, G. Biroli, J.-P. Bouchaud, L. Cipelletti, and W. van Saarloos, *Dynamical heterogeneities in glasses, colloids, and granular media* (OUP Oxford, 2011).

[6] 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 (2014).

[7] A. Widmer-Cooper, P. Harrowell, and H. Fynewever, How reproducible are dynamic heterogeneities in a supercooled liquid?, *Phys. Rev. Lett.* **93**, 135701 (2004).

[8] E. R. Weeks, Introduction to the colloidal glass transition, *ACS Macro Lett.* **6**, 27 (2017).

[9] D. Coslovich and G. Pastore, Understanding fragility in supercooled Lennard-Jones mixtures. I. Locally preferred structures, *J. Chem. Phys.* **127**, 124504 (2007).

[10] C. Patrick Royall, S. R. Williams, T. Ohtsuka, and H. Tanaka, Direct observation of a local structural mechanism for dynamic arrest, *Nature Mater.* **7**, 556 (2008).

[11] A. Widmer-Cooper, H. Perry, P. Harrowell, and D. R. Reichman, Irreversible reorganization in a supercooled liquid originates from localized soft modes, *Nat. Phys.* **4**, 711 (2008).

[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, *Phys. Rev. Lett.* **113**, 157801 (2014).

[13] R. L. Jack, A. J. Dunleavy, and C. P. Royall, Information-theoretic measurements of coupling between structure and dynamics in glass formers, *Phys. Rev. Lett.* **113**, 095703 (2014).

[14] F. Turci, C. P. Royall, and T. Speck, Nonequilibrium phase transition in an atomistic glassformer: The connection to thermodynamics, *Phys. Rev. X* **7**, 031028 (2017).

[15] C. P. Royall and S. R. Williams, The role of local structure in dynamical arrest, *Phys. Rep.* **560**, 1 (2015).

[16] H. Tanaka, H. Tong, R. Shi, and J. Russo, Revealing key structural features hidden in liquids and glasses, *Nat. Rev. Phys.* **1**, 333 (2019).

[17] 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?, *J. Chem. Phys.* **155**, 074502 (2021).

[18] M. Sharma, M. K. Nandi, and S. M. Bhattacharyya, Identifying structural signature of dynamical heterogeneity via the local softness parameter, *Phys. Rev. E* **105**, 044604 (2022).

[19] M. L. Manning and A. J. Liu, Vibrational modes identify soft spots in a sheared disordered packing, *Phys. Rev. Lett.* **107**, 108302 (2011).

[20] 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, *Phys. Rev. Mater.* **4**, 113609 (2020).

[21] M. Baggioli, I. Kriuchevskyi, T. W. Sirk, and A. Zaccone, Plasticity in amorphous solids is mediated by topological defects in the displacement field, *Phys. Rev. Lett.* **127**, 015501 (2021).

[22] Z. W. Wu, Y. Chen, W.-H. Wang, W. Kob, and L. Xu, Topology of vibrational modes predicts plastic events in glasses, *Nat. Commun.* **14**, 2955 (2023).

[23] E. 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, *Phys. Rev. Lett.* **114**, 108001 (2015).

[24] G. Jung, R. M. Alkemade, V. Bapst, D. Coslovich, L. Filion, F. P. Landes, A. J. Liu, F. S. Pezzicoli, H. Shiba, G. Volpe, et al., Roadmap on machine learning glassy dynamics, *Nat. Rev. Phys.* **7**, 91 (2025).

[25] S. S. Schoenholz, E. D. Cubuk, D. M. Sussman, E. Kaxiras, and A. J. Liu, A structural approach to relaxation in glassy liquids, *Nat. Phys.* **12**, 469 (2016).

[26] G. Jung, G. Biroli, and L. Berthier, Predicting dynamic heterogeneity in glass-forming liquids by physics-inspired machine learning, *Phys. Rev. Lett.* **130**, 238202 (2023).

[27] V. 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, *Nat. Phys.* **16**, 448 (2020).

[28] H. Shiba, M. Hanai, T. Suzumura, and T. Shimokawabe, *BOTAN: BOND Targeting Network for prediction of slow glassy dynamics by machine learning relative motion*, *J. Chem. Phys.* **158**, 084503 (2023).

[29] X. Jiang, Z. Tian, K. Li, and W. Hu, A geometry-enhanced graph neural network for learning the smoothness of glassy dynamics from static structure, *J. Chem. Phys.* **159**, 144504 (2023).

[30] F. S. Pezzicoli, G. Charpiat, and F. P. Landes, Rotation-equivariant graph neural networks for learning glassy liquids representations, *SciPost Phys.* **16**, 136 (2024).

[31] E. Boattini, F. Smallenburg, and L. Filion, Averaging local structure to predict the dynamic propensity in supercooled liquids, *Phys. Rev. Lett.* **127**, 088007 (2021).

[32] D. Teney, M. Peyrard, and E. Abbasnejad, Predicting is not understanding: Recognizing and addressing underspecification in machine learning, in *European Conference on Computer Vision* (Springer, 2022) p. 458.

[33] A. Swain, S. A. Ridout, and I. Nemenman, Machine learning that predicts well may not learn the correct physical descriptions of glassy systems, *Phys. Rev. R* **6**, 033091 (2024).

[34] K. Swanson, S. Trivedi, J. Lequieu, K. Swanson, and R. Kondor, Deep learning for automated classification and characterization of amorphous materials, *Soft Matter* **16**, 435 (2020).

[35] F. Font-Clos, M. Zanchi, S. Hiemer, S. Bonfanti, R. Guerra, M. Zaiser, and S. Zapperi, Predicting the failure of two-dimensional silica glasses, *Nat. Commun.* **13**, 2820 (2022).

[36] N. Oyama, S. Koyama, and T. Kawasaki, What do deep neural networks find in disordered structures of glasses?, *Frontiers in Physics* **10**, 1007861 (2023).
[37] S. Ciarella, D. Khomenko, L. Berthier, F. C. Mocanu, D. R. Reichman, C. Scalliet, and F. Zamponi, Finding defects in glasses through machine learning, Nat. Commun. **14**, 4229 (2023).

[38] G. 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, Europhys. Lett. **143**, 17004 (2023).

[39] G. 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, Phys. Rev. Mater. **8**, 025602 (2024).

[40] M. Liu, N. Oyama, T. Kawasaki, and H. Mizuno, Classification of solid and liquid structures using a deep neural network unveils origin of dynamical heterogeneities in supercooled liquids, J. Appl. Phys. **136**, 144702 (2024).

[41] G. Jung, G. Biroli, and L. Berthier, Dynamic heterogeneity at the experimental glass transition predicted by transferable machine learning, Phys. Rev. B **109**, 064205 (2024).

[42] Q. Wang, L.-F. Zhang, Z.-Y. Zhou, and H.-B. Yu, Predicting the pathways of string-like motions in metallic glasses via path-featurizing graph neural networks, Sci. Adv. **10**, eadk2799 (2024).

[43] K. Yoshikawa, K. Yano, S. Goto, K. Kim, and N. Matubayasi, Graph neural network-based structural classification of glass-forming liquids and its interpretation via self-attention mechanism, J. Chem. Phys. **163**, 024508 (2025).

[44] Q. Wang, J. Ding, L. Zhang, E. Podryabinkin, A. Shapeev, and E. Ma, Predicting the propensity for thermally activated β events in metallic glasses via interpretable machine learning, npj Computational Materials **6**, 194 (2020).

[45] C. Liu, Y. Wang, Y. Wang, M. Islam, J. Hwang, Y. Wang, and Y. Fan, Concurrent prediction of metallic glasses' global energy and internal structural heterogeneity by interpretable machine learning, Acta Materialia **259**, 119281 (2023).

[46] R. M. Alkemade, E. Boattini, L. Filion, and F. Smallenburg, Comparing machine learning techniques for predicting glassy dynamics, J. Chem. Phys. **156**, 204503 (2022).

[47] C. Rudin, Stop explaining black box machine learning models for high stakes decisions and use interpretable models instead, Nat. Mach. Intell. **1**, 206 (2019).

[48] D. C. Montgomery, E. A. Peck, and G. G. Vining, Introduction to linear regression analysis (John Wiley & Sons, 2021).

[49] A. Widmer-Cooper and P. Harrowell, Predicting the long-time dynamic heterogeneity in a supercooled liquid on the basis of short-time heterogeneities, Phys. Rev. Lett. **96**, 185701 (2006).

[50] S. J. Wetzel, S. Ha, R. Iten, M. Klopotek, and Z. Liu, Interpretable machine learning in physics: A review, arXiv preprint arXiv:2503.23616 10.48550/arXiv:2503.23616 (2025).

[51] S. Kivelson and S. Kivelson, Understanding complexity, Nat. Phys. **14**, 426 (2018).

[52] P. M. Chaikin, T. C. Lubensky, and T. A. Witten, Principles of condensed matter physics (Cambridge University Press, 1995).

[53] H. Tanaka, A simple physical model of liquid-glass transition: Intrinsic fluctuating interactions and random fields hidden in glass-forming liquids, J. Phys.: Condens. Matt. **10**, L207 (1998).

[54] A. Sharma, C. Liu, and M. Ozawa, Selecting relevant structural features for glassy dynamics by information imbalance, J. Chem. Phys. **161**, 184506 (2024).

[55] M. L. Falk and J. S. Langer, Dynamics of viscoplastic deformation in amorphous solids, Phys. Rev. E **57**, 7192 (1998).

[56] A. 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, Phys. Rev. E **97**, 033001 (2018).

[57] M. Lerbinger, A. Barbot, D. Vandembroucq, and S. Patinet, Rel-

evance of shear transformations in the relaxation of supercooled liquids, Phys. Rev. Lett. **129**, 195501 (2022).

[58] A. D. Parmar, M. Ozawa, and L. Berthier, Ultrastable metallic glasses in silico, Phys. Rev. Lett. **125**, 085505 (2020).

[59] D. Frenkel and B. Smit, Understanding molecular simulation: from algorithms to applications (Elsevier, 2023).

[60] L. Berthier and W. Kob, The Monte Carlo dynamics of a binary Lennard-Jones glass-forming mixture, J. Phys.: Condens. Matt. **19**, 205130 (2007).

[61] H. Shiba, Y. Yamada, T. Kawasaki, and K. Kim, Unveiling dimensionality dependence of glassy dynamics: 2d infinite fluctuation eclipses inherent structural relaxation, Phys. Rev. Lett. **117**, 245701 (2016).

[62] J. Behler and M. Parrinello, Generalized neural-network representation of high-dimensional potential-energy surfaces, Phys. Rev. Lett. **98**, 146401 (2007).

[63] J. W. Rocks, S. A. Ridout, and A. J. Liu, Learning-based approach to plasticity in athermal sheared amorphous packings: Improving softness, APL Materials **9**, 021107 (2021).

[64] J. Behler, Constructing high-dimensional neural network potentials: a tutorial review, Int. J. Quantum Chem. **115**, 1032 (2015).

[65] T. Kawasaki, T. Araki, and H. Tanaka, Correlation between dynamic heterogeneity and medium-range order in two-dimensional glass-forming liquids, Phys. Rev. Lett. **99**, 215701 (2007).

[66] C. F. Schreck, C. S. O'Hern, and L. E. Silbert, Tuning jammed frictionless disk packings from isostatic to hyperstatic, Phys. Rev. E **84**, 011305 (2011).

[67] H. Tong and H. Tanaka, Revealing hidden structural order controlling both fast and slow glassy dynamics in supercooled liquids, Phys. Rev. X **8**, 011041 (2018).

[68] C. Rycroft, VORO++: A three-dimensional Voronoi cell library in C++, Chaos **19**, 041111 (2009).

[69] F. H. Stillinger and T. A. Weber, Packing structures and transitions in liquids and solids, Science **225**, 983 (1984).

[70] G. James, D. Witten, T. Hastie, and R. Tibshirani, An introduction to statistical learning: with applications in R (Springer, 2013).

[71] D. W. Marquardt and R. D. Snee, Ridge regression in practice, Am. Stat. **29**, 3 (1975).

[72] H. Zou and T. Hastie, Regularization and variable selection via the elastic net, Journal of the Royal Statistical Society Series B: Statistical Methodology **67**, 301 (2005).

[73] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, et al., Scikit-learn: Machine learning in Python, J. Mach. Learn. Res. **12**, 2825 (2011).

[74] H. Peng, F. Long, and C. Ding, Feature selection based on mutual information criteria of max-dependency, max-relevance, and min-redundancy, IEEE Transactions on pattern analysis and machine intelligence **27**, 1226 (2005).

[75] L. Berthier and R. L. Jack, Structure and dynamics of glass formers: Predictability at large length scales, Phys. Rev. E **76**, 041509 (2007).

[76] J. Rottler and C. Ortner, Analysis of local structure of mechanical and thermal rearrangements in glasses with the atomic cluster expansion, J. Phys. Chem. B **128**, 11492 (2024).

[77] H. Akaike, A new look at the statistical model identification, IEEE transactions on automatic control **19**, 716 (2003).

[78] G. Schwarz, Estimating the dimension of a model, Ann. Statist. **6**, 461 (1978).

[79] G. H. Golub and C. F. Van Loan, Matrix computations (JHU Press, 2013).

[80] L. N. Trefethen and D. Bau, Numerical linear algebra (SIAM, 2022).