Machine Learning for the identification of phase-transitions
in interacting agent-based systems: a Desai-Zwanzig example

Abstract

Deriving closed-form, analytical expressions for reduced-order models, and judiciously choosing the closures leading to them, has long been the strategy of choice for studying phase- and noise-induced transitions for agent-based models (ABMs). In this paper, we propose a data-driven framework that pinpoints phase transitions for an ABM- the Desai-Zwanzig model - in its mean-field limit, using a smaller number of variables than traditional closed-form models. To this end, we use the manifold learning algorithm Diffusion Maps to identify a parsimonious set of data-driven latent variables, and show that they are in one-to-one correspondence with the expected theoretical order parameter of the ABM. We then utilize a deep learning framework to obtain a conformal reparametrization of the data-driven coordinates that facilitates, in our example, the identification of a single parameter-dependent ODE in these coordinates. We identify this ODE through a residual neural network inspired by a numerical integration scheme (forward Euler). We then use the identified ODE – enabled through an odd symmetry transformation– to construct the bifurcation diagram exhibiting the phase transition.

I Introduction:

Complex dynamic phenomena are ubiquitous in natural sciences, social sciences, and engineering [1, 2, 3, 4]. In many cases, their study has been performed using agent-based models (ABM) also known as interacting particle systems (IPS). The dynamics of those models in the thermodynamic (mean-field) limit undergo phase transitions – bifurcations in non-linear dynamics terminology  [5, 6]. The exploration of the dynamics of these systems typically involves extensive numerical simulation scenarios for a very large number, $N$, of interacting agents that can be truly challenging and impedes the widespread utilization of these models. Therefore, coarse-graining methodologies become necessary in order to reduce this complexity. Reducing the dimensionality of an ABM can be achieved by defining collective variables that are capable of accurately describing the full dynamics of such large systems in terms of a relatively small number of observables [7, 8, 9]; or discovering those by data mining [10, 11, 12].

In our previous work [9], we developed and successfully tested a model reduction approach based on the cumulants of the single-agent probability distribution for such large ABM systems. The proposed coarse-graining framework was based on an analytical closure methodology of the infinite hierarchy of equations for the moments or, equivalently, cumulants of the probability distribution of the infinite-dimensional system. The basic steps for building a reduced-order Desai-Zwanzig (DZ) model are: (i) Consider the mean-field ansatz by writing the $N$-particle distribution function, the solution of the $N-$particle Fokker-Planck equation, as the product of one-particle distribution functions $\rho(x,t)$ [13]. (ii) Represent $\rho(x,t)$ either in terms of its moments (DZ [14]) or in terms of its Fourier coefficients (Smoluchowski/noisy Kuramoto model [15]), thus obtaining an infinite system of ordinary differential equations (ODEs) that is exactly equivalent to the McKean-Vlasov PDE. Truncating the equations for the moments (or the Fourier coefficients), while choosing an appropriate closure scheme gives the reduced-order model described and analyzed in [9]. However, selecting the correct closure scheme for the truncated moments’ system is the trickiest part of the coarse-graining procedure.

Figure 1: Schematic of the overall workflow. (I) Sample data from the Agent-Based Model (ABM) across multiple initial conditions and parameter values, (II) Compute histograms for each snapshot of the ABM. (III) Apply the Diffusion Maps algorithm on the computed histograms to discover a reduced latent embedding. (IV) Use a conformal autoencoder (AE) to find a conformal reparametrization of the latent space. (V) Identify a data-driven ODE in terms of the latent coordinate $\nu_{2}$ of the AE. (VI) Construct the bifurcation diagram (enabled via a symmetry transformation) in terms of the latent coordinate $\nu_{2}$.

This analytical approximation requires choosing the “correct” observable(s), the right level of observation/analytical closure, and –importantly – the level at which we attempt the closure. To overcome these constraints, this paper proposes a data-driven framework for studying phase transitions of ABMs by (i) discovering the coarse observables in a data-driven fashion, (ii) determining the level of closure, and (iii) identifying the reduced dynamics directly from data, (iv) utilizing the data-driven reduced model and applying an odd symmetry transformation to it for the construction of the bifurcation diagram. To this end, we employ Diffusion Maps [16], a manifold learning algorithm, conformal autoencoders [17] and residual neural networks inspired by numerical integrators of ordinary differential equations [18, 19, 20, 21, 12]. This allows us to circumvent the difficulties arising from the choice of an analytical model and the selection of its closure. The main steps of our proposed data-driven workflow, are illustrated in Figure 1.

As an example, we apply our data-driven framework to the generalization of the DZ model with multiplicative noise considered in [9]; this very well-studied ABM model features a second-order (continuous) phase transition, which enables us to compare the performance of our data-driven approach with well-known analytical and computational results. We emphasize, however, that our method is quite general, since it does not depend on the detailed features of the interaction, and can be applied to many different ABMs that exhibit phase transitions in the thermodynamic limit.

II The Agent-Based Model:

We consider a model describing an ensemble of $N$ identical interacting agents subject to multiplicative noise. For this model, the agents are coupled via a mean reverting force and a second-order phase transition exists at a critical temperature that can be calculated analytically [9][Eqn. A12]. This transition appears as a symmetric pitchfork bifurcation of the mean-field dynamics. The dynamics of each agent, $x_{i}$, are described by a stochastic differential equation (SDE) of the form

$$\mathrm{d}x_{i}=\bigg{[}-x_{i}^{3}+(\alpha+\nu\sigma_{m}^{2})x_{i}-\theta(x_{i}-\overline{x})\bigg{]}dt+\sqrt{\sigma^{2}+\sigma_{m}^{2}x_{i}^{2}}\mathrm{d}W_{i};$$ (1)

where $\sigma$ is the bifurcation parameter, $\theta$ denotes the interaction strength, $\alpha$ is a parameter characterizing the amplitude of the multiplicative noise, $\nu$ corresponds to different mathematical prescriptions of the SDEs (Itô, Stratonovich, etc.) and, $\sigma_{m}$ is a rectifying parameter: when $\sigma_{m}\neq 0$ the phase-transition is pushed to higher values of $\sigma.$ The agents are coupled through $\bar{x}$, which denotes the center of mass of the system (equal to the first moment $M_{1}$). Furthermore, $dW_{i},i=1\dots N$ denotes independent Brownian motions. The values of the parameters were set to $\alpha=1$, $\theta=4$ , $\sigma_{m}=0.8$, $\nu=\frac{1}{2}$, and $N=12,000$ as in [9].

For this model, it was demonstrated in [9] that, away from the phase transition at $\sigma\cong 1.890$, the first moment, i.e., the order parameter of the system, is sufficient for an accurate description of the mean-field dynamics. As we get closer to the phase transition/the bifurcation, more moments need to be taken into account in order to accurately locate the bifurcation. More specifically, it was numerically demonstrated that at least a four-moment truncation is necessary in order to accurately predict the bifurcation/phase transition. In this work, by observing the eigenvalues of the Jacobian at the steady state (based on the moments equation proposed in [9]), we confirm that a clear separation of time scales prevails in the neighborhood of the bifurcation, and that the long-term dynamics live on a one-dimensional manifold. Our goal here, discussed in the next section, is to find a data-driven parametrization of this one-dimensional manifold, and identify the dynamics on it.

III Numerical Results:

III.1 Latent data-driven coordinates

Figure 2: (a) The first two moments ($M_{1}$, $M_{2}$) estimated from the ABM simulations colored with the parameter value $\sigma$ at which they were obtained. (b-c) Diffusion Maps on histograms for a single value of the parameter $\sigma$: (b) One-to-one relation between leading histogram moment $M_{1}$ and leading Diffusion Maps coordinate $\psi_{1}$; (c) shows the residual ($r_{k}$) based on the local linear regression algorithm indicating that $\psi_{1}$ might be enough to parametrize the data (larger value of $r_{k}$). (d-f) Diffusion Maps on collected histograms from the ABM simulation across multiple values of the parameter $\sigma$: (d-e) show the non-harmonic diffusion map coordinates $(\psi_{1},\psi_{2})$ colored with $\sigma$ and $M_{1}$ respectively; (f) shows the residual ($r_{k}$) based on the local linear regression algorithm indicating that $\psi_{1}$ and $\psi_{2}$ might be enough to parametrize the data (larger values of $r_{k}$).

We start by sampling data on an equidistant grid of $18$ distinct values of the parameter $\sigma$ in the interval $\sigma\in[0.5,2.2]$. A more detailed description of the sampling strategy is discussed in Section A.1.

Given the sampled data from the ABM we compute the first four moments $M_{1},M_{2},M_{3},M_{4}$ as $M_{k}=\frac{1}{N}\sum_{i=1}^{N}x_{i}^{k}(t)$, where $x_{i}^{k}$ denotes the $i$-th agent at a fixed time raised in the power $k$, $M_{k}$ is the $k$-th moment of the population from the sampled ABM data. The computation of the moments serves as a comparison between our data-driven framework and the one proposed in [9] where an analytical system of ODEs based on the four moments is proposed. The choice of four moment is also supported by studying the time-scale separation of the analytical reduced model proposed in [9] between the four-moment truncation and higher-order moment truncations of the dynamics, see A.5 in the Appendix. In Figure 2 we plot the first two moments $M_{1}$ and $M_{2}$ colored with the parameter $\sigma$. Based on this, we argue that the data are at least two-dimensional, and that the parameter $\sigma$ appears to correlate with the moments $M_{1},M_{2}$.

We then proceed by computing Diffusion Maps (see Section A.2 for a description of the algorithm) using data from the sampled ABM simulation. The Diffusion Maps algorithm was computed using: (a) histograms obtained for a single value of the parameter $\sigma$, (b) histograms obtained across our grid of 18 distinct $\sigma$ values, (c) on the four computed moments from data sampled across our 18 values of $\sigma$; those computations are reported in Section A.2.3 of the Appendix. We would like to reiterate that the Diffusion Maps computations on the moments are not necessary for the overall framework they serve as a comparison of our approach and the one by [9] et al.

Computing Diffusion Maps (based on histograms) for single values of $\sigma$ (far from the phase transition) gave, a one-dimensional manifold parametrized by $\psi_{1}$. In Figure 2 we plot the leading Diffusion Maps coordinate $\psi_{1}$ against the first moment $M_{1}$ and demonstrate that they are one-to-one. This suggests that Diffusion Maps is capable of discovering a coordinate that is one-to-one with the known order parameter $M_{1}$ for this model [14, 22, 9]. The claim that the manifold is one-dimensional in this case is also supported from the results of the local-linear algorithm shown in Figure 2; the residual $r_{k}$ of the first eigenvector $\psi_{1}$ is larger than the one of the remaining. In Section A.2.3 of the Appendix we illustrate the relation of the Diffusion Maps coordinate with the subsequent moments $M_{2},M_{3},M_{4}$.

We proceed with the Diffusion Maps computation on data sampled over multiple parameter values of $\sigma$. Details for the Diffusion Maps algorithm applied on histograms are provided in Section A.2.2. Here in Figures 2, 2 we show the two, non-harmonic eigenvectors that parameterize different eigendirections $\psi_{1}$ and $\psi_{2}$ colored with the parameter $\sigma$ and the first moment $M_{1}$. The selection of the non-harmonic eigenvectors was made by using the local linear regression algorithm proposed in [23] and implemented by the datafold Python package [24]. The residual $r_{k}$ of the first two eigenvectors, shown in Figure 2, is greater than the residual of the remaining eigenvectors. This suggests that the first two eigenvectors ($\psi_{1},\psi_{2}$) are non-harmonic, while the remaining ones are harmonic (functions) of the first two.

The Diffusion Maps results for the case of using collected data across multiple parameter values, shown in Figures 2, 2, and 2, suggest that the mean-field dynamics of the ABM live on a two-dimensional manifold. However, the parameter $\sigma$ appears to correlate with the behavior of the moments observed at that $\sigma$ value. This suggests that even though the data can be parametrized by two coordinates $\psi_{1}$, $\psi_{2}$, one might be able to disentangle this coupled two-dimensional description into a factorized (a single state)$\times$(a parameter) description. In the next section, we illustrate how a data-driven reparametrization of the Diffusion Maps coordinates can be achieved that will allow us to disentangle the effect of the parameter $\sigma$ from that of the latent variable.

III.2 Y-shaped conformal autoencoder

Figure 3: (a) A schematic of the Y-shaped conformal autoencoder. The inputs to the network ($\psi_{1}$ and $\psi_{2}$) are shown as green nodes. The outputs of the autoencoder ($\hat{\psi}_{1}$ and $\hat{\psi}_{2}$) and the estimator are shown as red nodes. The latent variables ($\nu_{1}$ and $\nu_{2}$) are shown as light blue nodes. (b-c) The obtained latent coordinates $\nu_{1},\nu_{2}$ colored with $\sigma$ and $M_{1}$ respectively. We can see that the $\sigma$ appears to vary only across $\nu_{1}$, compared to Fig. 2 (c) in which $\sigma$ varied across both $\psi_{1},\psi_{2}$. Similarly, $M_{1}$ appears to vary only across $\nu_{2}$. (d) The parameter $\sigma$ is plotted against the latent coordinate $\nu_{1}$ indicating a strong dependence. (e) Contour lines representing level sets of $\nu_{1}$ and $\nu_{2}$ are plotted in the Diffusion Maps ($\psi_{1}$,$\psi_{2}$) providing a visual illustration of the obtained conformality.

In this section, we illustrate the use of the Y-shaped conformal autoencoder, originally introduced in our previous work [17] as a means of disentangling the effect of different parameter combinations on model outputs in the context of parameter nonidentifiability. The Y-shaped conformal autoencoder seeks a reparametrization of the Diffusion Maps coordinates in which it disentangles the effect of the parameter $\sigma$ from that of the latent coordinates of the autoencoder. A more detailed description of the overall framework, the different loss function components of the network, and the training procedure we followed are discussed in Section A.3. A schematic of the Y-shaped conformal autoencoder is shown in Figure 3. The coordinates $\psi_{1}$ and $\psi_{2}$, obtained by computing the Diffusion Maps algorithm on estimated histograms of the ABM, were the ones used as input to the Y-shaped conformal autoencoder.

The bottleneck latent variables $\nu_{1},\nu_{2}$ are shown in Figures 3 and 3 colored with the parameter $\sigma$ and the first moment $M_{1}$ (that is the known coarse variable for this model [9]). We see that the parameter $\sigma$ appears to be varying only along $\nu_{1}$, compared to Figure 2 in which $\sigma$ varied across both $\psi_{1},\psi_{2}$. Similarly, $M_{1}$ varies only along $\nu_{2}$. The latent variable $\nu_{1}$ shows a strong correlation with the parameter $\sigma$, as depicted by Figure 3, indicating that they are effectively one-to-one. The conformality efficiency of the autoencoder becomes visually pronounced in Figure 3 where level sets of $\nu_{1}$ and $\nu_{2}$ are plotted in the Diffusion Maps coordinates. Furthermore, the average value of the cosine of the angle, $\cos\theta=\frac{\nabla\nu_{1}\cdot\nabla\nu_{2}}{||\nabla\nu_{1}||||\nabla\nu_{2}||}$, calculated over the test set, we found it to be $5\times 10^{-4}$.

This reparametrization of the latent coordinates is crucial since, as we show in the next section, it allows us to identify a single-state dimension ODE depending on $\sigma$ (rather than a needlessly two-dimensional ODE).

In Section A.3 we provide additional results obtained for the Y-shaped autoencoder illustrating the reconstruction of its input data ($\psi_{1},\psi_{2}$). Remarkably, in addition, the Y-shaped autoencoder allows us to estimate $\sigma$ from previously unseen histogram observation through $\nu_{1}$.

III.3 Identifying parameter dependent ODEs

Using the latent variable $\nu_{2}$ obtained from the Y-shaped conformal autoencoder as the state variable, we now identify a $\sigma$-dependent ODE; remember that $\nu_{2}$ is conformal to $\nu_{1}$ (and therefore to $\sigma$). The identification of the parameter-dependent ODE was achieved through a residual neural network inspired by the forward Euler numerical integration scheme. A schematic of the neural network is shown in Figure 4. The inputs to this network are the state variable at time $t$, $\nu_{2}(t)$, and the parameter $\sigma$, and the output is the state variable evolved to time $t+h$, $\nu_{2}(t+h)$. We provide a more detailed description of the forward Euler network and the constructed loss function scheme in Section A.4.

We test the ability of the ODE to produce accurate paths (trajectories), in $\nu_{2}$ space, for seven unseen values of the parameter $\sigma=\{0.57,0.85,1.11,1.75,1.9,2.06,2.25\}$, see Section A.1 for more details. In Figure 4 we illustrate the predicted values $\hat{\nu}_{2}(t+h)$ from the Euler neural network against the ground-truth (projected trajectories in $\nu_{2}$ generated by the ABM) for all seven unseen values of the parameter $\sigma$. The estimated MSE in this case was MSE = $6\times 10^{-4}$. In addition, in Figure 4 we provide a visual comparison between simulated by the Euler Neural Network trajectories and the ABM model (projected in $\nu_{2}$) for the test parameter values $\sigma=\{1.11,1.75,1.9\}$. The paths generated by the ABM and those generated by the neural network-identified ODE show excellent visual agreement across different values of $\sigma$ and different initial conditions. This suggests that the identified right-hand-side provides a successful approximation of the ground-truth dynamics.

Figure 4: (a) A schematic of the forward Euler residual neural network. The state variable $\nu_{2}(t)$, the parameter $\sigma$, and the time step $h$ are inputs to the neural network that estimates the right-hand side $f_{\theta}$ of the ODE. The right-hand-side is then used to estimate the state variable $\hat{\nu_{2}}(t+h)$ by using a forward Euler step. (b) The predicted value for $\hat{\nu}_{2}(t+h)$ estimated from the Euler neural network is plotted against the true value (projected ABM trajectories in $\nu_{2}$ space) for values of the parameter $\sigma=\{0.57,0.85,1.11,1.75,1.9,2.06,2.25\}$ not included in the training set. (c) For three values, of the parameter not included in the training set we contrast generated paths from the ODE (solid lines) with paths simulated by the ABM (dashed lines) embedded in $\nu_{2}$. The colors black, blue (gray) and green (light-gray) corresponds to $\sigma=1.11$, $\sigma=1.75$ and $\sigma=2.25$ respectively.

III.4 Bifurcation Diagram: phase transition

Figure 5: (a) A representative bifurcation diagram constructed by the identified right-hand-side of the Euler neural network suggests a “perturbed” pitchfork. (b) A representative bifurcation diagram after applying equation (2) to the identified right-hand-side shows a symmetric pitchfork.

In this section, we now proceed to test whether the parameter-dependent ODE identified in a data-driven way through our neural network can qualitatively and quantitatively capture the phase transition. To assess the robustness of our approach, we also estimate the error associated with the phase transition by training 5,000 neural networks in parallel, using different splits of the data for training, validation, and testing. All hyperparameters were kept the same across models during this computation. For each of the trained networks given the identified right-hand-side, we compute the steady-states across a range of values of the parameter $\sigma$ and construct the bifurcation diagram. For values of $\sigma\in[0.5,2.2]$ a representative constructed bifurcation diagram is shown in Figure 5. The bifurcation diagram shows that the model identifies the existence of three steady states, two stable and one unstable, for $\sigma<1.84$. For $\sigma\geq 1.84$ a unique stable steady state exists. However, the computed bifurcation diagram clearly possesses the flip symmetry of the generic pitchfork bifurcation. We have identified a perturbed pitchfork, in which the upper branch remains permanently stable and the lower branch exhibits a turning point where a stable and an unstable component collide. The consistency of identifying a perturbed pitchfork was maintained across the different models we trained.

To make sure we recover a generically symmetric pitchfork bifurcation diagram, a simple transformation suffices: given the identified right-hand-side $f_{\theta}(\nu_{2}(t);\sigma)$ of the neural network we compute

$$g(\nu_{2};\sigma)=\frac{f_{\theta}(\nu_{2}(t);\sigma)-f_{\theta}(-\nu_{2}(t);\sigma)}{2}.$$ (2)

The vector field $g(\nu_{2};\sigma)$ is then used to construct the bifurcation diagram, shown in Figure 5 correctly capturing the symmetry. Note that the computation of Equation (2) does not require a second neural network but only evaluations of the trained Euler neural network discussed in the previous section. The average critical parameter of the identified ODE (based on the $5000$ trained neural networks) was estimated at $\sigma^{*}=1.837$ and the standard deviation was $0.026$; the true critical transition has been computed to be at $\sigma^{*}=1.890$ [9]–see also[14][Eqn. 3.52] for the case of additive noise. This suggests that the identified ODE provides a qualitatively correct and arguably quantitatively accurate approximation of the phase transition. This small discrepancy can likely be attributed to the limited range of $\sigma$ values used in our training data. With only $18$ distinct values for $\sigma$ (details on sampling in Section A.1), the closest ones to the true bifurcation point are $\sigma=1.80$ and $\sigma=1.90$. Using a finer grid of parameter values we believe it would have led to even smaller discrepancy.

IV Summary and Conclusions

We presented a data-driven framework for identifying qualitatively and approximating quantitatively phase transitions of interacting agent systems through an interplay between direct simulations and machine-learning-assisted, coarse-grained system identification and bifurcation analysis. We demonstrated that our framework is capable of identifying the phase transition of the DZ system with a single, interpretable data-driven variable, in contrast to the derived closed-form model, proposed in [9] where a system of four approximate ODEs were required for the task. To discover the effective coarse (collective) variables that parameterize collected data from the ABM, the Diffusion Maps manifold learning algorithm [16] was used. We showed that, for a constant parameter value ($\sigma=1$), away from the phase transition, the Diffusion Maps algorithm discovers a one-dimensional manifold parametrized by the data-driven coordinate $\psi_{1}$; and that $\psi_{1}$ is one-to-one with the theoretical order parameter $M_{1}$. We then illustrated that, while data collected over a range of parameter values (bracketing the critical value) lie on a two-dimensional manifold, this manifold can be “disentangled” (factored) into a one-dimensional state-variability manifold “cross” a one-dimensional parameter variability manifold. To disentangle the variability due to the parameter $\sigma$ from the state variability, in terms of the latent Diffusion Maps coordinates ($\psi_{1},\psi_{2}$), we introduced a Y-shaped conformal autoencoder, initially proposed in [17] for addressing parameter non-identifiability. We illustrated that the autoencoder’s latent variables ($\nu_{1},\nu_{2}$) form a conformal set of coordinates that ultimately allow us to learn a single effective nonlinear ODE, in terms of a single state variable $\nu_{2}$ and depending on a single parameter $\sigma$ (the latter being one-to-one with $\nu_{1}$). This effective nonlinear ODE was identified via a residual neural network, templated on the forward Euler numerical integrator. We compared generated paths of the identified ODE for test values of the parameter $\sigma$ to paths generated with the full ABM embedded in the latent coordinate $\nu_{2}$ and demonstrated good agreement between the two. To construct the bifurcation diagram from the identified right-hand-side of the ODE, we used an odd symmetry transformation of the right-hand-side. This transformation provides a symmetric vector field and captures the pitchfork bifurcation that denotes the phase transition. Imposing known physical symmetries in the data-driven models in order to enhance accuracy and improve generalizability is currently an active area of research [25, 26, 27, 28, 29, 30, 31, 32, 33].

Our proposed framework can be extended to a broad class of complex, multiscale dynamical systems for which a fine scale, atomistic/stochastic mode exists, but for which accurate closed form macroscopic equations at the coarse-grained level are not explicitly known. We are particularly interested in models from the social sciences, where no physics-informed natural choice for the order parameter(s) exists. As an example, we mention the noisy Hegselmann-Krause model for opinion dynamics [34], for which it has been rigorously proved that a discontinuous disorder (no consensus)/order (consensus) phase transition exists [6][Prop. 6.2]. It would be interesting to identify the phenomenological order parameter introduced in [34] using our approach. Other examples for which our approach is expected to be applicable are the Keller-Segel model for chemotaxis, see the work in [35, 36, 37, 38], the Fitzhugh–Nagumo from Lattice Boltzmann in [39, 40], the mean-field Fitzhugh–Nagumo model for Neurons  [41][Sec. 5.6], and the Dean-Kawasaki stochastic PDE [3, 42], that takes into account finite-particle effects and fluctuations.

Codes

The code used to generate the results for this paper will become available as a public repository at https://gitlab.com/nicolasevangelou/ml_phase_transition upon acceptance of the paper.

Appendix A

A.1 Data Collection

Given the ABM model described in Section II in the main text, we generated trajectories at 18 equidistant values of $\sigma\in[0.5,2.2]$ [9]. For each value of $\sigma$ we sampled 100 trajectories for different initial conditions (agents’ distribution) drawn from the Pearson distribution (implemented in Matlab) with prescribed mean, standard deviation, skewness, and kurtosis. The values for the mean, standard deviation, skewness, and kurtosis were chosen randomly from a prescribed grid of equidistant points: The mean was chosen from the range $[-2.0,2.0]$, with increments of $0.2$. The standard deviation was selected from the range $[0.0,2.0]$, with increments of $0.1$. The skewness was chosen from the range $[-2.0,2.0]$, with increments of $0.2$ and the kurtosis was selected from the range $[0.0,15]$, with increments of $0.72$. This scheme ensured a dense sampling both, in parameter and state space.

The integration time of the ABM, containing $N=12,000$ agents, was set to $t_{f}=10$ with a time-step $dt=0.005$. Data were collected every five snapshots which led to a total of $400$ snapshots per trajectory. Therefore, for a single value of the parameter, the total number of initial data is $n=40,000$ ($400\times 100$) while for multiple values of the parameter $n=720,000$ ($400\times 100\times 18$). Note that a number of trajectories explodes and thus, are omitted from any further computations.

An additional preprocessing step was applied to the data before the Diffusion Maps computation that includes discarding a short transient $t_{cut}$ for each trajectory. This ensures that the fast transients have decayed and the collected data contains only the long-term dynamics (that live on the slow manifold). When dealing with multiple parameter values we chose $t_{cut}=1$ while for a single value of the parameter $\sigma=1$ we choose $t_{cut}=0.5$ to make sure that enough transients are close to the unstable steady state, otherwise, the manifold would appear as two clusters. As a test set, we sampled nine trajectories for each of the seven values of $\sigma=\{0.57,0.85,1.11,1.75,1.9,2.06,2.25\}$ not included in the training set. The selection of these test values was made to validate the predictions of our Euler neural network for dynamics before and after the bifurcation.

A.2 Diffusion Maps

The Diffusion Maps algorithm, introduced by Coifman and Lafon [16] can be used to discover a low-dimensional parameterization of high-dimensional data $\mathbf{X}=\{\bm{x}_{i}\}_{i=1}^{n}$ with each $\bm{x}_{i}\in\mathbb{R}^{m}$. Diffusion Maps constructs a weighted graph $\mathbf{K}\in\mathbb{R}^{n\times n}$ between the sampled data points by using a kernel function. A common choice, also used in our case, is the Gaussian kernel

$$K(\bm{x}_{i},\bm{x}_{j})=\exp\bigg{(}\frac{-||\bm{x}_{i}-\bm{x}_{j}||_{2}^{2}}{2\varepsilon}\bigg{)},$$ (3)

where $\varepsilon$ is a positive hyperparameter that controls the rate of the kernel’s decay. The metric $||\cdot||$ in our case was chosen as the $\ell^{2}$ norm but different metrics are also possible.

To discover a low-dimensional manifold, regardless of the sampling density, the following normalization is required

$$\mathbf{\tilde{K}}=\mathbf{P}^{-1}\mathbf{K}\mathbf{P}^{-1}$$ (4)

where $P_{ii}=\sum_{j=1}^{n}K_{jj}$. A second normalization of $\mathbf{\tilde{K}}$ recovers a row-stochastic, Markovian matrix

$$\mathbf{M}=\mathbf{D}^{-1}\mathbf{\tilde{K}}$$ (5)

where $\mathbf{D}$ is a diagonal matrix defined as $D_{ii}=\sum_{j=1}^{n}\tilde{K}_{ij}$.

The entries of matrix $\mathbf{M}$ can be seen as probabilities of jumping from one point to the other. The eigendecomposition of $\mathbf{M}$,

$$\mathbf{M}\bm{\phi}_{i}=\lambda_{i}\bm{\phi}_{i},$$ (6)

provides a set of eigenvectors $\bm{\phi}_{i}$ and corresponding eigenvalues $\lambda_{i}$. To obtain a more parsimonious representation of the original data set $\mathbf{X}$ proper selection of the eigenvectors is needed. If the intrinsic dimension of the data is small, this selection can be achieved by visual inspection of the non-harmonic eigenvectors (eigenvectors that span independent directions) [43]. Alternatively, the local-linear regression algorithm proposed by Dsilva et al. [23] can be used for selecting the non-harmonic eigenvectors. If, the number of independent non-harmonic eigenvectors is smaller than the dimension $m$ of the data $\mathbf{X}$, then dimensionality reduction has been achieved. In our work, the Python library datafold [24] was used for the Diffusion Maps and the local-linear regression algorithms.

A.2.1 Nystöm Extension

The Nyström Extension formula provides a numerical approximation of eigenfunctions of the form [44]

$$\int_{a}^{b}M(\bm{x}_{j},\bm{x}_{i})\phi_{i}(\bm{x}_{i})d\bm{x}_{i}=\lambda\phi({\bm{x}_{j}}).$$ (7)

In our work, Nyström extension is utilized when new out-of-sample data points are given, e.g., $\bm{x}_{new}\notin\mathbf{X}$. To generate the Diffusion Maps coordinates $\bm{\phi}_{new}$ Nyström extension uses an interpolation scheme based on the kernel computations and normalizations applied during the dimensionality reduction step, discussed in the previous Section. The Nyström extension formula reads

$$\phi_{i}(\bm{x}_{new})=\frac{1}{\lambda_{i}}\sum_{j=1}^{n}\tilde{\mathbf{M}}(\bm{x}_{new},\bm{x}_{j})\phi_{i}(\bm{x}_{j}),$$ (8)

where $\phi_{i}(x_{new})$ denotes the estimated $i$-th eigenvector for the data point $\bm{\phi}_{new}$, $\lambda_{i}$ denotes the corresponding eigenvalue, $\phi_{i}(\bm{x}_{j})$ denotes the $j$-th component of the $i$-th eigenvector and $\tilde{\mathbf{M}}(\cdot,\cdot)$ denotes the kernel function used to determine the similarity of $\bm{x}_{new}$ to all the points in $\mathbf{X}$.

A.2.2 Diffusion Maps on ABM data

In this section, we provide details on how the Diffusion Maps algorithm was computed on histograms and moments from the ABM. In the first case, for each snapshot of the ABM we constructed a histogram as an approximation of the agents’ density. Each histogram contains $40$ equidistant bins defined in the range $[-4,4]$. This range ensures that all agents in the collected (training) data lie in between. Note, that the method is insensitive to the selected number of bins. To reduce the computational cost of the Diffusion Maps, we subsampled the training data uniformly [24] which resulted in about $3,500$ data points when Diffusion Maps applied for a single parameter and about $19,000$ data points for multiple values of the parameter. The hyperparameter $\varepsilon$ was selected as the square of the median of the pairwise distances multiplied by a constant $c$, where $c=0.03$ for a single value of the parameter $\sigma$ and $c=20$ for multiple values of the parameter $\sigma$. In the second case, Diffusion Maps was computed on the sampled moments, $M_{1},M_{2},M_{3},M_{4}$. Again, to alleviate the computational cost we subsampled the data [24], which results in $N\sim 11,000$. The hyperparameter $\varepsilon$ was also selected for this case by computing the median of the pairwise distances multiplied with $c=10$.

A.2.3 Diffusion Maps: Additional Results

In this section, we provide additional results for the Diffusion Maps computation for $\sigma=1$ on the histograms and the Diffusion Maps computations on the moments $M_{1},M_{2},M_{3},M_{4}$. We showed in the main text that the Diffusion Maps coordinate $\psi_{1}$ is one-to-one with the first moment $M_{1}$. Here we show that $M_{3}$ is also one-to-one with $\psi_{1}$ (Figure 6) and that the even moments $M_{2}$ and $M_{4}$ can be seen as functions of $\psi_{1}$ (Figures 6,6).

To strengthen, the argument that the manifold in this case is one-dimensional we show in Figure 6 the eigenvectors $\psi_{2}-\psi_{9}$ plotted against the first non-trivial eigenvector $\psi_{1}$. From Figure 6 it appears that the eigenvectors $\psi_{2}-\psi_{9}$ are harmonics of $\psi_{1}$. This suggests that the manifold for a fixed value of $\sigma$ is one-dimensional.

Figure 6: (a-c) The Diffusion Maps coordinate $\psi_{1}$ is plotted against the computed moments $M_{2}$, $M_{3}$ and $M_{4}$ respectively. (d) The eigenvector $\psi_{1}$ is plotted against the eigenvectors $\psi_{2}-\psi_{9}$. This supports our argument that the eigenvectors $\psi_{2}-\psi_{9}$ are harmonics of $\psi_{1}$.

Figure 7: (a-b) The non-harmonic coordinates $\phi_{1}$ and $\phi_{2}$ colored with $\sigma$ and $M_{1}$ respectively. (c) The residual $r_{k}$ indicates $\phi_{1}$ and $\phi_{2}$ are the two non-harmonic coordinates.

The Diffusion Maps algorithm applied using the computed first four moments reveals a two-dimensional manifold embedded in four dimensions. The first two eigenvectors in this case $\phi_{1}$ and $\phi_{2}$, shown in Figures 7 and 7 are the non-harmonic, eigenvectors that parametrize the data. This is corroborated by the larger residuals $r_{k}$ of the first two eigenvectors shown in Figure 7. Also in this case, the parameter $\sigma$ and the first moment $M_{1}$ appear visually as functions of the Diffusion Maps coordinates (Figures 7 and 7).This pair of coordinates, obtained by Diffusion Maps on moments, could have been used for the CAE’s training instead of the Diffusion Maps coordinates ($\psi_{1},\psi_{2}$)obtained from the computations on the histograms; but we omit this for brevity.

A.3 Y-shaped conformal autoencoder

The Y-shaped conformal autoencoder was initially presented in [17]. In this work, the Y-shaped conformal autoencoder consists of the three connected (sub)networks

$$\text{Encoder : }(\psi_{1},\psi_{2})\mapsto(\nu_{1},\nu_{2})$$ (9)

$$\text{Decoder : }(\nu_{1},\nu_{2})\mapsto(\hat{\psi}_{1},\hat{\psi}_{2})$$ (10)

$$\text{Estimator : }\nu_{1}\mapsto\hat{\sigma}.$$ (11)

The Encoder receives as inputs the two Diffusion Maps coordinates $\psi_{1}$, $\psi_{2}$ and maps them to the latent variables $\nu_{1},\nu_{2}$. The Decoder aims to reconstruct the Diffusion Maps coordinates from the latent variables $\nu_{1}$, $\nu_{2}$. The Estimator has as input the latent coordinate $\nu_{1}$ and aims to learn a map from $\nu_{1}$ to the parameter $\sigma$.

The loss function used to train the Y-shaped conformal autoencoder consists of three parts: (a) The loss function of the Encoder-Decoder (autoencoder) $\mathcal{L}_{\text{ae}}$ that aims to reconstruct the input itself (b) the loss function of the Estimator, $\mathcal{L}_{\text{est.}}$, that aims to reproduce the parameter $\sigma$ given $\nu_{1}$ and (c) the loss function for imposing the conformality constrain, $\mathcal{L}_{\text{con.}}$, between $\nu_{1}$ and $\nu_{2}$,

$$\langle\nabla\nu_{1},\nabla\nu_{2}\rangle=0$$ (12)

where the gradient $\nabla$ is in terms of the Diffusion Maps coordinates ($\psi_{1}$,$\psi_{2}$) of the input and $\langle\cdot,\cdot\rangle$ denotes the inner product between the two vectors. The gradients were computed by using the automatic differentiation of Pytorch [45]. In practice, instead of using equation 12, one can minimize the angle between the vectors

$$\text{cos}\theta=\frac{\nabla\nu_{1}\cdot\nabla\nu_{2}}{||\nabla\nu_{1}||||\nabla\nu_{2}||}$$ (13)

which stabilizes the training of the network. We describe details for training the network, specifics about the architecture used, and the choice of hyperparameters in the next section.

A.3.1 Hyperparameter selection and training procedure

The implementation of the Y-shaped conformal autoencoder was done with the Pytorch Python library [45].

Each (sub)network (Encoder, Decoder, Estimator) in the architecture of the Y-shaped conformal autoencoder consists of five fully connected layers. The first four hidden layers have 20 neurons and tanh(t) activation functions and the fifth has no activation function (linear activation) and its size depends on the size of the desired output. The ADAM optimizer was chosen for training the overall network. We chose mini-batches of size $32$ to train the network. The learning rate was selected as $\eta=0.001$ and $500$ epochs. Training the network for a larger number of epochs leads to overfitting.

To train the network and test its generalization capability we used about $19,000$ data points. We split the data into train$|$test$|$validation as $80$:$10$:$10$. We then rescaled the training data by using the MinMaxScaler Python preprocessing scheme from sklearn. We applied the same transformation for the validation and test sets. During the training of the network, we used only the training set to perform backpropagation and the validation set to get insight into the model’s performance. The network did not see the test set during training.

The optimization process we performed was heuristic: for a fixed mini-batch two updates (backpropagation steps) were performed. The first step updates the weights of the Encoder-Decoder and the second step the weights of the Estimator-Encoder. Altering this training protocol is possible. In Algorithm 1 below we provide a more detailed description of the network’s training.

Upon training of the neural network, the estimated MSE for the autoencoder’s reconstruction, was $\mathcal{L}_{\text{ae}}=1.36e-08$ on the train set and $\mathcal{L}_{\text{ae}}=1.37e-08$ on the test set. The MSE for the Estimator was $\mathcal{L}_{\text{est.}}=1.2e-03$ for the train set and $\mathcal{L}_{\text{est.}}=1.2e-03$ for the test. The average value of the $\text{cos}\theta$, on the train set was $\mathcal{L}_{\text{con.}}=2.72e-05$ and on the test set $\mathcal{L}_{\text{con.}}=2.63e-05$.

Input: Diffusion Maps coordinates $\psi_{1}$,$\psi_{2}$ and parameter values $\sigma$.

Output: The weights of: (i) Encoder ($\theta_{\text{encoder}}$),

(ii) Decoder ($\theta_{\text{decoder}}$), (iii) Estimator ($\theta_{\text{estimator}}$).

For $i=1,2,...,T$

1. Predict:

   $$(\nu_{1},\nu_{2}) =\text{Encoder}(\psi_{1},\psi_{2})$$

   $$(\hat{\psi}_{1},\hat{\psi}_{2}) =\text{Decoder}(\nu_{1},\nu_{2})$$

2. Compute Autoencoder (Encoder-Decoder) and Conformality Losses:

   $$\mathcal{L}_{1} =\mathcal{L}_{\text{ae}}+\mathcal{L}_{\text{con.}}$$

   $$=\text{MSE}(\hat{\bm{\psi}},\bm{\psi})+\alpha\text{MSE}(\text{cos}\theta,0)$$

3. Backpropagation step - update weights (illustration with gradient descent):

   $$\theta_{\text{encoder}}-=\eta{\theta_{\text{encoder}}}\mathcal{L}_{1}$$

   $$\theta_{\text{decoder}}-=\eta{\theta_{\text{decoder}}}\mathcal{L}_{1}$$

4. Predict:

   $$(\nu_{1},\nu_{2})=\text{Encoder}(\psi_{1},\psi_{2})$$

   $$\hat{\sigma}=\text{Estimator}(\nu_{1})$$

5. Compute Estimator Loss:

   $$\mathcal{L}_{\text{est.}}=\text{MSE}(\hat{\sigma},\sigma)$$

6. Backpropagation step - update weights (illustration with gradient descent)

   $$\theta_{\text{estimator}}-=\eta{\theta_{\text{estimator}}}\mathcal{L}_{\text{est.}}$$

   $$\theta_{\text{encoder}}-=\eta{\theta_{\text{encoder}}}\mathcal{L}_{\text{est.}}$$

Algorithm 1 The algorithm illustrates a full iteration during training of the Y-shaped conformal autoencoder. We set the scale parameter to $\alpha=10$. The learning rate is denoted as $\eta$.

A.3.2 Y-shaped conformal autoencoder: Additional Results

In this section, we provide additional results regarding the Y-shaped conformal autoencoder described in the main text. The ability of the Estimator to predict the parameter $\sigma$ from $\nu_{1}$ is shown in Figure 8 for train (black dots) and test (red dots) points. We also illustrate, in Figures 8, 8, the reconstruction of the autoencoder for train and test points.

Figure 8: (a) The true values of the parameter $\sigma$ are plotted against the reconstructed by the Estimator for train (black) and test (red) points. The blue dashed line indicates $y=x$. (b-c) The true values of the Diffusion Maps coordinates $\psi_{1}$ and $\psi_{2}$ are plotted against the reconstructed $\hat{\psi}_{1}$, $\hat{\psi}_{2}$ coordinates by the autoencoder.

A.4 Forward Euler neural network

In this section, we describe how we identified the right-hand side of an ordinary differential equation directly from data. Let $\nu_{2}(t)$ be a state variable whose dynamics are governed by a $\sigma$-dependent ODE given by the general form

$$\dot{\nu}_{2}(t)=f(\nu_{2}(t);\sigma)$$ (14)

Our goal is to construct a neural network architecture inspired by numerical integrators of ODEs to estimate the right-hand-side $f(\nu_{2}(t);\sigma)$. To this end, we constructed a forward Euler residual neural network depicted in Figure 4. To train this network, we do not require long trajectories but only snapshots of the form $\mathcal{D}=\{\nu_{2}(t+h),\nu_{2}(t),\sigma,h\}$ where $\nu_{2}(t)$ is the state variable at time $t$, $\nu_{2}(t+h)$ is the state variable after a small time-step $h$. Given sampled data in the form of $\mathcal{D}$ we wish to approximate $f$ by using a neural network, with weights denoted as $\theta$.

To formulate the loss used to train the network $f_{\theta}$ we remind the reader that the forward Euler approximates the evolution of an ODE, by a small positive step $h$, as

$$\nu_{2}(t+h)=\nu_{2}(t)+hf(\nu_{2};\sigma).$$ (15)

In our case, without having access to $f$ but only data in the form $\mathcal{D}$ we wish to approximate the right-hand-side with the neural network $f_{\theta}$. This is achieved by performing for each pair of inputs ($\nu_{2}(t),\sigma$) one integration step of size $h$, estimating the evolved dynamics $\hat{\nu_{2}}(t+h)$ and minimizing the loss

$$\mathcal{L}(\theta|\nu_{2}(t),\nu_{2}(t+h),h,\sigma)=||\hat{\nu}_{2}(t+h)-\nu_{2}(t+h)||^{2}.$$ (16)

For our computations, the time step $h$ was kept constant but the overall approach can be easily extended to handle also varying time steps $h$.

A.4.1 Hyperparameter selection and training procedure

The implementation of the forward Euler neural networks was done with the Tensorflow/Keras Python libraries [46].

To train the forward Euler neural network we need to ensure our data is in the form of snapshots $\mathcal{D}$. To achieve that we used the Nyström extension formula (see Section A.2.1) to all the available sampled trajectories and obtained the corresponding trajectories in $\psi_{1}\psi_{2}$. We then evaluated the Encoder to get trajectories in terms of $\nu_{2}$. These two steps provided us with about $600,000$ snapshots. We then split the data into train$|$test$|$validation as 80$:$10$:$10. We then centered and whitened the data (based on the mean and variance of the training set) and applied the same transformation to the validation and test set.

The architecture consisted of two hidden layers with $10$ neurons each. The first hidden layer had a tanh(t) activation function and the second hidden layer had a linear activation function. We used ADAM to optimize this network. The mini-batch size was set to $32$, the number of total epochs to $100$, and the learning rate to $\eta=0.001$. The learning curves for the training and validation are shown in Figure 9. The ability of the network to fit the training set and generalize is shown in 9. Upon training of the network, the MSE on the train and test sets were $2.63e-03$ and $2.79e-05$ respectively. This neural network was used for the computations reported in III.3. We did not record the MSE for all the $5000$ networks used to check the robustness of our approach in identifying the critical transition.

A.4.2 Euler neural network: Additional Results

Figure 9: (a) The learning curves for the train and validation sets. (b) The true values of the state variable $\nu_{2}(t+h)$ plotted against the predicted by the Euler neural network $\hat{\nu}_{2}(t+h)$ for the train (black points) and test (red points).

In this section, we provide some additional results for the forward Euler neural network. The learning curves for the training and validation are shown in Figure 9. The ability of the network to fit the train set and generalize is shown in Figure 9. Note that the normalized value of the state variable $\nu_{2}$ is shown in Figure 4 and that the test points in this include values of the parameter $\sigma$ that are also in the training set.

A.5 Separation of time scales

In this Section, we discuss the separation of time scales for the reduced order models based on the equations for the moments proposed in [9]. As shown in Figure 10 the slowest eigenvalue $\lambda_{1}$ is at least one order of magnitude smaller than the second slowest $\lambda_{2}$. This suggests that the dynamics after a short transient are slaved in the direction of this slowest eigenvalue and are effectively one-dimensional.

Figure 10: The eigenvalues of the Jacobian computed across a range of parameter values $\sigma$ for the moments’ equations with six ODEs presented in [9]. The slowest eigenvalue (upper line) $\lambda_{1}$ crosses the real axis at $\sigma=1.89$ where the bifurcation occurs. The black square indicates the value of $\sigma$ where the (slowest) eigenvalue crosses zero.