A machine learning approach to drawing phase diagrams of topological lasing modes

Identifying phases and analyzing the stability of dynamic states are ubiquitous and important problems which appear in various physical systems. Nonetheless, drawing a phase diagram in high-dimensional and large parameter spaces has remained challenging. Here, we propose a data-driven method to derive the phase diagram of lasing modes in topological insulator lasers. The classification is based on the temporal behaviour of the topological modes obtained via numerical integration of the rate equation. A semi-supervised learning method is used and an adaptive library is constructed in order to distinguish the different topological modes present in the generated parameter space. The proposed method successfully distinguishes the different topological phases in the Su-Schrieffer-Heeger (SSH) lattice with saturable gain.This demonstrates the possibility of classifying the topological phases without needing for expert knowledge of the system and may give valuable insight into the fundamental physics of topological insulator lasers via reverse engineering.


I. INTRODUCTION
Over the last few years, significant research efforts have been made on photonic topological insulator (PTI) lasers.While the efforts have been concentrated on the spatial stability of the topologically protected edge modes, namely on the existence of such topological edge modes in non-Hermitian PTIs [1][2][3][4][5][6][7][8][9], the temporal stability has not been the focus of interest so far [10][11][12][13][14][15].Due to the non-linear nature of PTI lasers, the temporal stability is an important characteristic to take into account for experimental demonstrations and real-life applications.Although the spatial stability of the topological modes, i.e., its robustness, may be guaranteed in active non-Hermitian PTIs, unstable behaviour may be present in the time domain.In this regard, the temporal dynamics of the topologically protected modes have been studied [10,12,13], mainly using linear stability analysis.It is, however, a challenging task to apply the same approach to more complex structures because of the highdimensional phase space and parameter space as well as the lack of analytical solutions [11].
Machine learning (ML) can be advantageous for the theoretical study of the stability of PTIs, which requires repetitive numerical simulations for several varying parameters.ML is a data-based method which can be implemented with different strategies, and the most appropriate ML strategy depends on the dataset under study.For instance, a supervised learning strategy relies on labelled data, a dataset of input-output pairs.This has been utilized in topological photonics to draw topological phase diagrams [16], calculate topological invariants [5], or explore topological band structures [17].On the other hand, an unsupervised learning strategy consists of extracting information from the dataset for which we do not have labels.This is used for dimensional reductions by keeping only the main features of the high-dimensional structure of the dataset or for clustering problems from which the data is divided into different types [18].For instance, this has been successful in obtaining the phase transition in the Ising model [19], and clustering Hamiltonians that belong to the same symmetry classes [20].
In the unsupervised learning strategy, modal decomposition is a common and successful method which reduces the analysis of very high-dimensional data to a set of relatively few modes.Among the modal decomposition methods, principal component analysis (PCA) is a method which derives the eigenmodes or the main features based on their variance in the data [21].These eigenmodes can then be utilized as a basis to represent the dataset [22].This reduced-order model method has been extended to identify distinct non-linear regimes [23][24][25][26][27][28][29] by constructing a library composed of representatives of these regimes: this is known as representation classification.Nevertheless, the preliminary identification of the regimes composing the library and the construction of the library is a manual process and requires expert knowledge of the complex system.
In this paper, we propose a representation classification method to study the spatio-temporal dynamics of non-linear topological systems.The results will be based on the phase diagram of the Su-Schrieffer-Heeger (SSH) lattice [30] with a domain wall and with saturable gain [12,13].To remove the necessity of the required expertise on the complex system, we present an algorithm which constructs an appropriate library of the different phases automatically.For this goal, we propose two approaches: a top-down approach in which the library has numerous phases that are merged into the equivalent phases, and a bottom-up approach in which the library is completed on the fly to get the most accurate classification.Via reverse engineering, our proposed method can be used as a tool to find novel topological lasing modes in more complicated settings.For given rate equations of a lasing system, one would only need to integrate the differential equations in the desired parameter space region and then apply the adaptive representation classification to obtain the phase diagram.

A. Phase diagram of the SSH Model
As a toy model, we will consider the domain-wall-type SSH lattice with saturable gain [Fig.1(a)].The system has a domain wall, at the A site n = 0, which separates two SSH lattices, namely lattices composed of two sites per unit cell, A and B, and characterized by intraand inter-unit cell couplings, t intra and t inter , respectively.t intra = t 1 and t inter = t 2 (t intra = t 2 and t inter = t 1 ) for the lattice on the left (right) side of the domain wall, i.e., the sites with n < 0 (n ≥ 0).The dynamics of the system ψ(t) = (ψ −N (t), . . ., ψ −1 (t), ψ 0 (t), ψ 1 (t), . . ., ψ N (t)) ≡ x(t), with N s = 2N + 1 sites, reads, for n = −N, . . ., N : where ψ n is the amplitude of the n-th site.g n and γ n are the linear gain and linear loss at the n-th site, respectively.Using explicitly the amplitudes a p and b p of the A and B sites on the p-th unit cell, respectively, in ψ(t) = (. . ., a p (t), b p (t), . ..), Eq. ( 1) can be re-written as: where g σ and γ σ are the linear gain and linear loss at the site σ = A, B.
In the passive setting (g A = g B = 0, γ A = γ B = 0), this configuration, with t 1 > t 2 , is known to give a single topologically protected zero-energy (nonoscillating) mode localized at the domain wall and with non-vanishing amplitudes only on the A sites [31].This is due to the bulk-boundary correspondence at the domain wall between trivial and non-trivial topological SSH lattices.Indeed, an SSH lattice is topologically trivial (nontrivial) if the intra-unit cell coupling is lower (greater) than the inter-unit cell coupling.
In the active setting, it has been shown that the topological phase of the system depends on the gain and coupling parameters [12,13].As parameters, we use the values from Ref. [13] with t 1 = 1, t 2 = 0.7, g B = 0 and γ A = γ B ≡ γ AB .Figure 1 of the A (B) sites in Fig. 1(c) as in Ref. [12,13].Alternatively, more details can be understood by plotting the space-time dynamics of the topological modes as shown in Figs.1(d) and 1(e).The non-oscillating phase is similar to the zero-energy mode in the passive SSH lattice.We can see in Fig. 1(d) that the mode is localized at the interface and has the majority of its amplitudes on the A sublattice.On the other hand, the system with saturable gain exhibits a new topological phase with no counterpart in the passive setting.The new topological mode is characterized by an edge mode at the domain wall with an oscillating behaviour of the amplitudes on the A and B sites, as shown in Fig. 1(e).
The classification of the new topological phases in nonlinear systems requires, so far, an expert knowledge of the given non-linear systems, for example, the known results derived in Ref. [12].In fact, the phase diagram in Fig. 1(b) has been obtained solely by the fast Fourier transform of the time series in the parameter space.Thus, the main aim of this paper is to develop a tool to explore the topological phases of PTI lasers in more complicated settings for which we have little knowledge.
The phase diagram shown in Fig. 1(b) will serve as a reference for our proposed method.The dataset we will utilize throughout this paper is composed of about 1000 samples which are randomly generated from the same coupling and gain parameters' range as in Fig. 1(b).The coupled-mode equations [Eqs.(2) and (3)] are integrated using the fourth-order Runge-Kutta method and with a p (t = 0) = b p (t = 0) = 0.01, ∀p, as an initial condition.Although the integration has been performed using a fixed time step dt = 0.01 until a final time at t = 1400, only 2000 time snapshots are uniformly retrieved in order to keep the time series at a reasonable size.For the parameters given above, this sample rate leaves about 10 time steps per period for the oscillating regime case [Fig.1(c)].The phase diagram is then obtained solely from the time series of the states within the given parameter space.

B. Representation classification method
To classify topological lasing modes based on their distinct non-linear regimes, we use a representation classification method [23,27,28].The general idea of representation classification relies on the assumption, and common situations, that the dynamics of a high-dimensional system evolves on a low-dimensional attractor such as fixed points or periodic orbits [32].The low-dimensional structure of the attractor allows for a reduced-order model that accurately approximates the underlying behaviour of the system: the dynamics of the complex system can thus be written using a basis that spans the lowdimensional space.Representation classification consists of constructing a library of appropriate basis, representative of the dynamical regimes of interest, and only then employ a filtering strategy to identify the regime corresponding to a given unknown time series.In the following, we will use the term "regime" to denote the different dynamical behaviours or the different topological phases in the non-linear SSH lattice with saturable gain.Besides, for convenience, we will plot only the total intensity on the A (I A = p |a p | 2 ) and B (I B = p |b p | 2 ) sublattices to represent the given regimes.Nevertheless, the time series of the complex amplitudes at each site will be considered for the construction of the library.
As is common in complex dynamical systems, the dynamics of a system close to an attractor lie in a lowdimensional space.This means that a given spatiotemporal dynamics, denoted by the vector x(t), can be approximately written in terms of a basis Φ = {φ i } i=1,...,D spanning the low D-dimensional space, namely: where β i are the weighted coefficients in the above linear combination of basis states φ i .Using the terminology used in the literature [23,27,28], x(t) will, in the following, be referred to the state measured at time t.However, one of the main characteristics of non-linear systems is the drastically different dynamical behaviours with respect to the system's parameters.Therefore, the reduced-order modelling strategy using a single representative basis, i.e., corresponding to a single regime, is bound to fail.Instead of finding a global basis, we here construct a set of local bases, i.e., construct a library composed of the bases of each non-linear regime of interest: where J is the number of regimes, Φ j 's are the bases of each of the dynamical regime j, and φ j,i 's are the corresponding basis states.This is the supervised learning part of the method, from which the data-driven method attempts to capture the dynamics of the system in the reduced-order model.Therefore, the library L contains the representative basis of each regime of interest, and corresponds to an overcomplete basis that approximates the dynamics of the system across the given parameter space.A better approximation of x(t), instead of using Eq. ( 4) for a single basis regime, then reads: where β j,i are the weighted coefficients in the above linear combination in the overcomplete basis library L. It is worth noting that the library modes φ j,i are not orthogonal to each other, but instead orthogonal in groups of modes for each different regime j.
Throughout this paper, the bases used for constructing the library L will be generated by using a timeaugmented dynamical mode decomposition (aDMD) method [28] to consider both the spatial and temporal behaviours (see Supplementary section SI for additional information).
Here, we use a classification scheme based on a simple hierarchical strategy [28].The regime classification approach is fundamentally a subspace identification problem, where each regime is represented by a different subspace.Given the state x(t i : t i+Nw ) measured within the time window [t i , t i+Nw ], with N w the time step window size, the correct regime j * is identified as the corresponding subspace in the library L closest to the measurement in the L2-norm sense [28].In other words, the classification strategy is to find the subspace that maximises the projection of the measurement onto the regime subspace: where P j is a projection operator given by: , and arg max the function that returns the index of the maximum value.
For the 1D system exemplified in this paper, the data collection is not too expensive because N s is reasonable.However, if N s is very large, sparse measurement might be desirable and a slight change in the methodology is then needed as explained in Supplementary section SII.

C. Phase Diagrams
In the following we will draw phase diagrams using different approaches to the representation classification.In the phase diagrams we will mark the different identified regimes by the color of the dots, where the (dark or light) blue dots always mark the oscillating regime, the green dots the non-oscillating regime, the red dots the transient regime and the orange dots the transition regime.The white and grey areas are overlays of the referenced phase diagram obtained in Fig. 1.The aDMD bases have been generated with N w = 25 from the time series starting at the 1800-th time step.choice.The remaining coloured dots in the plot represent the identified regime j * [Eq.( 7)] of each sample depicted by the dots in the parameter space.However, we can see that the phase diagram fails to correctly predict all the dynamical behaviours.Indeed, we observe that many time series are not correctly identified.Using a different choice of topological regimes in the parameter space to construct the library could be a solution to find a better phase diagram, but our attempts only showed marginal improvement of the agreement.Testing every possible choices in the dataset with no guarantee of finding the accurate phase diagram is therefore not a practical solution.Instead, using four bases in the library instead of two bases, or equivalently considering four regimes from the given parameter space, the phase diagram in Fig. 2(b) shows better results: The identified oscillating and nonoscillating regimes are more separated and have a better fitting with the referenced phase diagram, even though they belong to a distinct regime index j * .Merging the three oscillating regimes present in the library would then give a more accurate derived phase diagram.Therefore, Figure 2 suggests that increasing the number of bases in the library L and merging some of them might help to get closer to the desired phase diagram, as we will see in the later sections.

Top-down adaptive library
In the previous section, the construction of the library L was a manual process from which we already knew the different dynamical regimes.This, however, requires prior knowledge of the complex system considered.The strategy, here, is to adaptively construct the library based on the given data samples.Here, we employ a top-down approach in which we start with too many samples for the construction of the library, and then reduce the library size by merging some of them.Based on some measures in the decision process, this automated construction of the library thus removes the manual construction of the regimes.
The underlying assumption of the classification scheme is based on the dissimilarity between the subspace of different regimes.We thus propose to consider regimes that are similar as equivalent regimes.This would, for instance, help us to merge the three phases in the nonoscillating region in Fig. 2(b), and consider them as a single regime.In other words, the regimes i and j are said to be equivalent, denoted by i ∼ j, if the fraction of information retained after the projection onto each other, γ ij ∈ [0, 1], is high enough: where γ th ∈ [0, 1] is the hyper-parameter which decides the threshold value for merging different regimes and γ ij is the subspace alignment given by: with Importantly, the relation Eq. ( 9) is numerically computed in such a way that the transitivity property of the equivalence relation is satisfied, namely that if i ∼ j and j ∼ k then i ∼ k.The relation Eq. ( 9) is then indeed an equivalence relation because the reflexive (i ∼ i) and symmetric (i ∼ j ⇒ j ∼ i) property of the relation is automatically satisfied from the definition of γ ij [Eq.(10)].Supplementary section SIII gives more details on the top-down library generation principle.
The top-down representation classification strategy is to classify the time series according to a large library of bases, and only then merge the equivalent identified regimes via the calculated alignment subspace γ ij and the equivalence relation Eq. (9). Figure 3 shows the phase diagram obtained from the top-down algorithm that started with a library composed of J = 60 regimes randomly chosen, along with the representative of each phase.We observe in Fig. 3 However, we can see that the derived phase diagram is still failing in the low γ AB and low g A − γ AB region (bottom-left region of the present phase diagram), where some time series are interpreted as non-oscillating instead of oscillating regime.This shows the limitation of this method where the initially constructed library may lack some of the paths that may connect similar bases.For example, regimes i and k might not be similar enough to be considered as equivalent directly [Eq.( 9)] but are both equivalent to the regime k, i.e., i ∼ j and j ∼ k, which is missing in the library.The natural workaround would be to increase the initial library size and ensure that the regimes in the library have no missing paths, as we will see in the next section.
The hyper-parameter γ th is an important quantity in the algorithm since it dictates which regimes are equivalent or not.A low threshold γ th will easily merge regimes while a high γ th will barely reduce the size of the library as depicted in Fig. 4(a).The threshold is here arbitrarily chosen based on Fig. 4, and based on the refinement of the desired library.For example, we can see in Fig. 4(b) that the derived phase diagram with γ th = 0.55 has two different phases.For this coarse threshold, the transient regime is identified but the distinct non-oscillating and oscillating phases are merged together into a single phase.On the other hand, with the same library as in Fig. 4(b), Figure 4(c) displays the obtained phase diagram for a finer threshold value γ th = 0.95.The plot shows that the algorithm separates the parameter space into several regimes which can be grouped into four main regimes.In addition to the non-oscillating, the oscillating and the transient regimes, there is a regime corresponding to the transition between the two topological phases.Besides, this finer description allows us to see distinct sets of modes in the oscillating parameter space region [dark blue and light blue dots in Fig. 4(c)].Nevertheless, we observe again that the initial library misclassifies some of the non-oscillating time series most likely because of some missing paths, as said previously.

Bottom-up adaptive library
We propose an alternative and dual approach which considers fewer samples in the library.The core idea of this bottom-up approach is then to add samples on the fly during the classification of the given sample if the library is not good enough.
Here, the library is considered to be good enough if the maximal projection of the measurement onto the regimes' subspace is high enough.In other words, the library is said to be good enough if the worst relative reconstruction error, , is low enough: where th is the hyper-parameter which decides the threshold quality of the library and := max j=1,...,J with • 2 the L 2 -norm of a vector.Supplementary section SIV gives more details on the bottom-up library generation principle.The advantage of this bottom-up approach is the full exploration of the parameter space region and the automatic construction of a library based on its quality.This method does not suffer from the randomly chosen samples used to construct the library, and the library composition is not restricted to a narrow parameter space region.Using a good enough library quality, the algorithm should therefore be able to sort the missing paths issue in the top-down method.
The bottom-up representation classification scheme consists of classifying the time series according to a given library or adding this sample into the library if the library is not good enough, and only then merging the different phases obtained into groups of equivalent regimes using the top-down method.Figure 5 depicts the phase diagram derived from the bottom-up classification algorithm with a starting library composed of a single regime.Similarly to the top-down approach in Fig. 3(a), we observe three distinct regimes corresponding to the nonoscillating, oscillating and transient regimes [Fig.3(b)].Nevertheless, the obtained phase diagram now better predicts the regimes.The misclassifications of the nonoscillating and oscillating regimes, which were due to missing paths in the library, are now reduced, and only very few dots are not correctly identified due to being close to the topological transition boundary.Likewise, the transient points are indications of longer simulations needed because of the long transient time.
Along with the γ th hyper-parameter, the threshold hyper-parameter th is an important parameter since it tells us whether we want to add or not a given sample into the library.We observe in Fig. 6 that a low threshold th will add many samples to the library, whereas a high th will not add samples to the library at all.The threshold value th is, again, arbitrarily chosen but with a preference for a high-quality library, i.e., low th , in order to avoid missing paths.For example, we can see in Fig. 6(c) the phase diagram derived using th = 0.05 (and γ th = 0.95), namely with a library that gives less than 5% of the reconstruction error of the measurement.This set of hyper-parameter gives four main regimes that seem to correspond to the non-oscillating, oscillating, transition and transient regimes.Yet, there is some misidentification of the two topological phases most likely because of missing paths of the obtained library.On the other hand, with a better library quality, here th = 0.005 (and γ th = 0.95), the missing paths are retrieved and the derived phase diagram correctly predicts the topological phases [Fig.6(b)].Figure 6(b) shows that the different regimes, obtained previously with a lower library quality, are now better defined.The non-oscillating and oscillating regimes are well located in their respective parameter space region, and the transition points follow the transition boundary between the two topological phases.In addition, the bottom-up representation classification is predicting distinct oscillating modes [dark blue and light blue dots in Fig. 6(b)].
The presence of distinct oscillating modes is an example of new insights given by the data-driven classification method.Indeed, the complex values of the amplitudes of x(t) are, here, taken into account instead of solely the total intensity of each sublattice A and B as in Ref. [12,13].This allows for a finer description of the dynamic pattern based on the whole lattice with the relative phase difference of the sites or the absolute value of amplitudes.

III. CONCLUSION
We have proposed a data-driven approach to identify and classify topological phases of dynamical systems.By utilizing the representation classification strategy based on the aDMD, we have successfully drawn the phase diagrams of the domain-wall-type SSH lattice with saturable gain.To avoid manual labelling in the classification, we have proposed two automatic construction schemes: top-down and bottom-up approaches that merge similar phases in a library or adaptively construct a library according to its quality, respectively.While the bottom-up method is preferred due to the missing paths issue, both methods allow exploration of parameter space without any expert knowledge of the complex non-linear system.Our approach is advantageous in doing the re-verse engineering to find new phases because only a small desired parameter space region is required for solving the laser rate equation, applying the bottom-up representation classification strategy and building a starting library.Our method opens the door for finding novel topological lasing modes by providing new insights into the dynamics of coupled lasers in more complicated settings.

SI. BASIS GENERATION METHODS
This supplementary section review few methods for generating a basis from a time series.There are many ways to generate the low-order model of a given dynamical behaviour.An important quantity is the so-called data matrix X built from the data at hand.The data matrix is a (N s × N t )-matrix that collects the N t data snapshots x(t i ) into columns: Here, the vector x(t i ) is chosen to be the complex-valued amplitudes of the modes at the A and B sites.Other "observables", such as the absolute values or the total intensity per sublattices, can be used.This may give different (better or worse) results and is left for a future study.The bases can then be constructed using dimensional reduction techniques on the data matrix.Here, we will cover different methods in order to highlight their importance and limitations for the classification scheme used.

A. Proper orthogonal decomposition
Proper orthogonal decomposition (POD) [24] is a commonly used tool for dimensional reduction of physical systems.This decomposition relies on the singular value decomposition (SVD) of the data matrix, given by: where U and V † are (N s × N s ) and (N s × N t ) unitary matrices, respectively.Σ is a diagonal (N 0 × N 0 )-matrix diag = (σ 1 , . . ., σ N0 ), with N 0 = min(N s , N t ).The diagonal entries of Σ are the so-called singular values and are ordered in ascendant order σ 1 > σ 2 > . . .> σ N0 ≥ 0.
The SVD gives us two orthonormal bases U and V † since the matrices U and V † are unitary matrices.The columns of U are ordered according to the variance σ i they capture in the data matrix and are called the singular vectors: these are the POD modes that are used in the basis Φ.Moreover, the POD modes are complex because of the complex data X.
For a low-dimensional attractor, the POD basis can be safely truncated at a cut-off value r while retaining the main information of the data matrix.Explicitly, the SVD reads: and keeping only the r highest terms in the decomposition [Eq.S3], we have the approximation: This is re-written, in a matrix form, as: where U r , Σ r and V † r are the truncated matrix of U , Σ and V † , respectively.Although the cut-off value r can be chosen based on different criteria [33], r is typically chosen so that the POD modes retain a certain amount of the variance (or energy) σ X in the data, namely: Figure S1 displays the POD method of the nonoscillating and oscillating regimes in the domain-wall SSH lattice with saturable gain [Fig.1].The truncation has been chosen such that σ X = 99% of the total variance is retained.We observe in Figs.S1(a) and S1(b) the normalised singular values, and that a single POD mode is retained for the zero-mode-like, whereas three POD modes are needed for the oscillating regime, as marked by the red dots.The real and imaginary parts of the field profile of the corresponding first few POD modes are plotted in Figs.S1(c) and S1(d) and Figs.S1(e) and S1(f) for the zero-mode-like non-oscillating and oscillating regimes, respectively.One can see that the main spatial feature of the zero-mode is captured in the single POD mode obtained after truncation, where the majority of its amplitudes are on the A sublattice.On the other hand, the POD modes of the oscillating dynamical regime also capture part of the information with some finite amplitudes on the A and B sublattices.
Importantly, in this decomposition [Eq.S3], SVD is implicitly doing a space-time separation of the data matrix, where the POD modes U contain the spatial information while V have the temporal information at each spatial grid point.Therefore, the POD modes give a static basis and do not explicitly model the temporal dynamics of the time series.This method will therefore most likely fail to identify the correct dynamical regime in the classification step.

B. Dynamical mode decomposition
Dynamical mode decomposition (DMD) [24,34,35] is an alternative to the POD method for learning the dynamics of non-linear systems.DMD is an explicitly temporal decomposition that takes the sequences of snapshots into account, and is able to derive the spatiotemporal patterns of the data matrix X.The dynamics of the system are taken into account by considering a linear matrix A that maps a data matrix X 1 , starting at some time steps (t 1 ), to the data matrix X 2 , starting at the next time step (t 2 ).The matrix A is thus defined as: and the corresponding data matrices are given by: and Interestingly, Equation S7 is similar to a linear stability analysis formulation for discrete maps if we think of the stability matrix as the linear matrix A. The DMD method thus consists of solving the following eigenvalues problem: where the columns of Φ are the DMD modes φ i and the corresponding DMD eigenvalues λ i are the diagonal entry of Λ.The DMD modes φ i give us the spatial eigenmodes while their corresponding eigenvalues λ i have their temporal information.Using a change of units from the data snapshots, observed at every ∆t, to units in time, the eigenvalues are complex-valued scalars: where µ i gives the growth (decay) rate if µ i > 0 (µ i < 0) and ω i the oscillation frequency of the DMD modes φ i .However, the size of the data matrix usually makes the eigendecomposition not feasible.The goal, here, is therefore to approximate the eigenvalues and eigenvectors of A, using only the data matrices X 1 and X 2 .The idea is to start by the truncated SVD of X 1 = U r Σ r V † r in which Eq. S7 becomes: Then the linear matrix A is reduced by considering its projection onto the truncated POD subspace: The eigenvalue problem for A r is solved with: from which we have the relation: The key feature of the DMD method is that it decomposes the data into a set of coupled spatio-temporal modes.The DMD resembles a mixture of the POD in the spatial domain and the discrete Fourier transform (DFT) in the time domain.Figure S2 shows the DMD results for the oscillating regime.We can, indeed, see in Fig. S2(a) that the singular values are similar to that of the POD.Besides, we observe in Figs.S2(c) and S2(d) that the field profile of the DMD modes closely resembles the POD modes in Fig. S1.The largest DMD modes not only look similar to the POD modes, but they also contain the oscillation frequencies from ω i , as in DFT.The DMD even goes beyond DFT by giving an estimate of the growth (decay) rate in time via µ i > 0 (µ i < 0).This can be seen by plotting the DMD modes, scaled by their contribution in the decomposition σ i , in the frequency plane of λ i .We can see in Fig. S2(b) that the dynamical regime has a single DMD mode with ω i = 0 akin to the offset of the oscillation amplitudes, and two DMD modes with opposite ω i = 0 corresponding to the oscillating behaviours.All the above three DMD modes have vanishing growth or decay rate µ i = 0.

C. Time-augmented dynamical mode decomposition
Although the DMD gives the temporal behaviours of the non-linear system, the temporal information is not fully incorporated into the DMD basis Φ since only the DMD modes are used.Exploiting the time evolution in the dynamical regime requires the use of DMD modes along with their eigenvalues.The idea is therefore to incorporate the dynamic information by augmenting the basis Φ [28].This time-augmented DMD will be denoted by aDMD in the remaining of this chapter.Using the defining relation of the eigenvalues λ i as similar to a time evolution operator, i.e. multiplying by λ is the same as shifting by one time step, we have for a given DMD mode φ i , its evolution given by λ Nw i φ i at N w time step ahead in time.Therefore, the time-augmented basis vector reads:  (c) and S3(d)] carry some temporal evolution information.In particular, we observe in the top panel of Figs.S3(c) and S3(d) the static behaviours of the aDMD mode with ω i = 0. On the other hand, we can see, in the bottom panel of Figs.S3(c) and S3(d), one of the first aDMD modes with ω i = 0 featuring some oscillating behaviour in time.The size of the basis mode is larger than the plain DMD, and can exhibit its time evolution.Nevertheless, the graphs do not exactly plot the temporal evolution of the DMD modes since the first N s entry of the basis state is for the N s sites; the next N s for again the N s sites but at the next time step, etc.

D. Classification results from different decomposition methods
The classification results [Eq.( 7)] from the decomposition methods reviewed is shown in Fig. S4.We observe in Figs.S4(a) and S4(b) that the phase diagrams fail to correctly predict the distinct dynamical regimes.This is expected since the POD or DMD modes do not contain enough information about the temporal behaviours.Besides, the classification for these diagrams is based on a single snapshot (N w = 0).Thus, it is expected the classification fails to capture the correct dynamics since a single snapshot only relies on the spatial pattern of the regime.On the other hand, we can see in Figs.S4(c) and S4(d) that the derived phase diagrams have better accuracy when the bases are time-augmented, or equivalently when the classification scheme uses several snapshots.By increasing the time window in the classification, the derived phase diagram is even better as illustrated in Figs.S4(c) and S4(d) for N w = 5 and N w = 25, respectively.The phase diagram will get improved until the time window is large enough to capture the dynamic behaviour.

SII. SPARSE MEASUREMENT
This section detailed the slight change in the methodology in case of sparse sensing.Sparse sensing is often desirable since the measurement and the data collection can be expensive for a complex system if the space grid is too fine, i.e. if N s is very large.The compressed mea-surement y(t) is derived from the full-state measurement x(t) and the measurement matrix C: where C is a matrix of size (N p × N s ) with N p the number of measurements.Although the measurement matrix C can be represented by some advance and complex mapping [36], here we focus on point-wise measurements, namely the C ij entry in the matrix measurement corresponds to the i-th measurement at the j-th spatial grid point.Therefore the compressed basis is given by: where Φ is the basis obtained from the full-state data collection.The library of bases for the J distinct dynamical regimes is similarly re-written as: Nevertheless, the size of the current SSH lattice is, here, reasonable and allows us to choose the matrix measurement as the identity matrix C = 1 Ns .We will thus use the full-state instead of sparse measurements, but retain the Θ and y(t) notations to keep the general formalism.

SIII. TOP-DOWN LIBRARY GENERATION PRINCIPLE
This section supplement the top-down library generation principle.The general idea of the top-down construction of the library is to start with a library made of a high number of bases, generated from the time series randomly chosen in the given parameter space region, and then merge them into groups of equivalent regimes.Figure S5 illustrates the top-down algorithm.Starting with a library composed of J bases (here J = 60), the subspace alignment γ ij is computed [Fig.S5(a)] and then grouped into equivalent regimes according to Eq. ( 9) in the main text [Fig.S5(b)].Each representative of the regimes is then randomly selected within each group [Fig.S5(c)].We plot in Fig. S5(d) the time series of the representative of each regime.These regimes are the non-oscillating and oscillating regimes, as well as a third regime which may correspond to a transient regime.The vertical dashed line in the time series represents the initial time used for constructing the bases in the library.

SIV. BOTTOM-UP LIBRARY GENERATION PRINCIPLE
This section supplement the bottom-up library generation principle.Figure S6 illustrates the bottom-up methodology proposed.We start with a single sample in the library, randomly chosen in the parameter space region [Fig.S6(a)].The library is then adaptively constructed according to the relative reconstruction error (b) shows the phase diagram,

FIG. 1 .
FIG. 1. Phase diagram of the domain-wall-type SSH lattice with saturable gain.(a) Schematic of the domainwall-type SSH lattice considered.The vertical red line is a guide to the eye for the domain wall between the two SSH lattices, SSH 1 and SSH 2. (b) Phase diagram of the domainwall-type SSH lattice with saturable gain and linear loss on the A and B sublattices.The white and grey areas correspond to the non-oscillating and oscillating topological phases, respectively.(c) Representative total intensity IA (and IB) of the A (and B) sublattice in blue (and orange) for the nonoscillating and oscillating topological lasing mode on the top and bottom panel, respectively.Spatio-temporal dynamics of the (d) non-oscillating and (e) oscillating topological lasing modes.The non-oscillating and oscillating topological modes displayed are chosen at (γAB, gA − γAB) = (0.48, 0.06) and (0.16, 0.44), respectively.

Figure 2 (FIG. 2 .
Figure 2(a) displays the phase diagram derived from a library basis made of the two topological regimes known from Fig. 1.These two topological modes, used for the construction of the library, are represented by the two magenta dots.They are randomly chosen from the known regimes' region in Fig. 1(b) and the details of the resulted phase diagrams is dependent from that random

FIG. 3 .
FIG. 3. Representation classification based on a topdown adaptive library.(a) Phase diagram obtained using the top-down classification strategy with an initial library composed of J = 60 regimes randomly chosen and γ th = 0.75.(b) Representative total intensity IA (and IB) of the A (and B) sublattice in blue (and orange) for the different regimes.The black vertical dotted line indicates the starting time from which the bases are generated.
(a) that the derived phase diagram is able to distinguish the non-oscillating [top panel of Fig. 3(b)] and oscillating regimes [middle panel of Fig. 3(b)].In addition, a third regime corresponding to a transient regime [bottom panel of Fig. 3(b)] is found close to the γ AB = 0 or g A − γ AB = 0 axis.This transient regime indicates that a longer simulation time might be needed to be considered either in the oscillating or non-oscillating regimes.

FIG. 5 .
FIG. 5. Representation classification based on a bottom-up adaptive library.(a) Phase diagram obtained using the bottom-up classification strategy with a starting library composed of a single regime randomly chosen, th = 0.005 and γ th = 0.75.(b) Representative total intensity IA (and IB) of the A (and B) sublattice in blue (and orange) for the different regimes.

FIG. 6 .
FIG. 6.th -hyper-parameter dependency.(a) Library size against th .The inset is a zoom-in of the plot.Phase diagram derived using the bottom-up classification strategy with γ th = 0.95 and (b) th = 0.005 and (c) th = 0.05.The initial library is composed of a single regime randomly chosen.

FIG. S1 .
FIG. S1.POD method.Singular values of the (a) nonoscillating and (b) oscillating regimes.The red dots correspond to the singular values accumulating 99% of the total variance of the data.(c) Real and (d) imaginary parts of the field profile of the first POD mode for the non-oscillating regime.(e) Real and (f) imaginary parts of the field profile of the (top) first and (bottom) second POD mode for the oscillating regime.The POD bases have been generated from the time series starting at the 1800-th time step.
FIG. S2.DMD method.(a) Singular values of the oscillating regime in the DMD.The red dots correspond to the singular values accumulating 99% of the total variance of the data.(b) Plot of the logarithm of the DMD values ln(λ) in the complex plane.The size of the open circle is proportional to their corresponding singular values.(c) Real and (d) imaginary parts of the field profile of the (top) first and (bottom) second DMD mode for the oscillating regime.The DMD basis has been generated from the time series starting at the 1800-th time step.

)
FIG. S3. aDMD method.(a) Singular values of the oscillating regime in the aDMD.The red dots correspond to the singular values accumulating 99% of the total variance of the data.(b) Plot of the logarithm of the aDMD values ln(λ) in the complex plane.The size of the open circle is proportional to their corresponding singular values.(c) Real and (d) imaginary parts of the "field profile" of the (top) first and (bottom) second aDMD mode for the oscillating regime.Here, by "Site n", we mean the n-th entry of the eigenvector.The aDMD basis has been generated with Nw = 25 from the time series starting at the 1800-th time step.
FIG. S4.Phase diagram derived using different decomposition methods.Phase diagrams obtained from a library composed of two regimes (one non-oscillating and one oscillating) from which the bases have been generated using the (a) POD, (b) DMD, (c) aDMD with Nw = 5 and (d) aDMD with Nw = 25.The green and blue dots correspond respectively to the identified non-oscillating and oscillating regimes.The magenta dots represent the regimes used for the construction of the library.These are located at (γAB, gA − γAB) = (0.48, 0.06) and (0.16, 0.44).The white and grey areas are overlays of the referenced phase diagram obtained in Fig. 1.The bases have been generated from the time series starting at the 1800-th time step.
FIG. S5.Principle of the adaptively top-down library generation.(a) Subspace alignment matrix for an initial library composed of J = 60 regimes.(b) Subspace alignment matrix from the same library as in (a) but grouped into equivalent regimes.(c) Subspace alignment matrix from the reduced library after the equivalent regimes are merged with γ th = 0.75 in (b).(d) Representative total intensity IA (and IB) of the A (and B) sublattice in blue (and orange) for the different regimes.The black vertical dotted line indicates the starting time from which the bases are generated.The aDMD bases have been generated with Nw = 25 from the time series starting at the 1800-th time step.

[
FIG. S6.Principle of the adaptively bottom-up library generation.(a) Total intensity IA (and IB) on the A (and B) sublattice in blue (and orange) for the single and randomly chosen regime composing the library.(b) Subspace alignment matrix from the adaptively constructed library with the single initial regime in (a) and with th = 0.01.(c) Subspace alignment matrix from the reduced library after the equivalent regimes are merged with γ th = 0.75 in (b).(d) Representative total intensity IA (and IB) of the A (and B) sublattice in blue (and orange) for the different regimes.The black vertical dotted line in (a) and (d) indicates the starting time from which the bases are generated.The aDMD bases have been generated with Nw = 25 from the time series starting at the 1800-th time step.