Almost-Linear RNNs Yield Highly Interpretable Symbolic Codes in Dynamical Systems Reconstruction

We appreciate the referee’s overall positive assessment and the valuable feedback provided!

Weaknesses

W1 (usefulness of symbolic dynamics): First, please note that the symbolic encoding itself is not probabilistic, i.e. the symbolic seq. (as shown in Figs. 13, 14 or 17) are as deterministic as the underlying system itself. As new Fig. R5 in the provided PDF now further highlights, quantities like the topological entropy obtained from the symbolic encoding correlate highly with quantities obtained directly from the system dynamics, like the max. Lyapunov exp., but are much easier to compute. This further confirms that important topological properties of the underlying system can be inferred from the symbolic encoding we used. In general, symbolic dynamics has led to many powerful insights and results about the dynamics of certain systems (certain proofs about chaos or the number of unstable periodic orbits could only be derived symbolically, for instance; see Wiggins 1988, Guckenheimer & Holmes, 1983).

Hence, we think the referee’s point concerns more the representation of symbolic dynamics in the form of transition graphs. This type of graph represent. is quite standard in symbolic dynamics (see, e.g., textbook by Lind & Marcus 2021), where arrows between nodes usually represent admissible transitions, just like in graph representations of finite state machines or formal languages, for example. The graphs are meant to represent the set of all possible seq. (or ‘syntactically correct’ sentences). We agree, however, that one needs to be clear about the semantics of these graphs and their interpretation, which we will clarify.

W2 (guidance on model size): We chose $P$ according to a grid search as the min. value at which performance started to plateau (because we wanted to obtain topologically minimal representations), so its choice is not arbitrary but related to the kinks (or humps) in the curves in Fig. 3. Likewise, $M$ was determined by systematic grid search, see Appx. Fig. 9. However, we like the referee’s idea of determining an optimal number of linear subregions by regularization. We now implemented this, and find that the numbers of subregions determined this way agree well with those obtained by our prev. crit., see Fig. R6 in PDF.

W3 (teacher forcing): Sparse teacher forcing is indeed only applied during training, not during testing. Hence it is, correctly, not included as part of those equations. This is in fact crucial: DSR models are supposed to be generative models that, after training, can generate new data with the same geometrical and topological structure as those produced by the observed system. During test time, therefore, the once trained model cannot rely on actual observations, but follows solely its own dynamics. While TF is a broad term, sparse TF is different (Brenner et al. 2022, Mikhaeil et al. 2022) and has been introduced specifically in the context of DSR where it is SOTA (see also Tab. R1). We will clarify this.

W4 (reproducibility): Please note that we provided code for reproducibility with the org. submission here: https://anonymous.4open.science/r/SymbolicDSR-1955/README.md

The hyperpar. settings for all runs were further listed under A2: Training Method. In the rev. we will further clarify how exactly these optimal hyperpar. were determined.

W5 (Th.1 incorrect): The th. is correct, but we see where this misunderstanding comes from: For the referee’s example, the system would need to be non-hyperbolic, i.e. would need to have a center subspace. However, this case (req. conjugate pairs of eigenvalues to lie exactly on the unit circle) we explicitly ruled out in our theorems (it is a 0-measure set in par. space). We mentioned this in the pg. above the th. (since common to all of them), but for clarity will make explicit in the theorems themselves.

W6 (theoretical sects.): We are happy to reduce sect. 3.2 and move it partly to the Appx., the ref. is right that some of the concepts are not followed up upon. Others, like that of a shift operator, shift space, or topolog. partition, are however directly used in the theorems in sect. 4 (for the sake of the theorems, we prefer to be formally precise about certain concepts even if they may be intuitively clear).

W7 (undefined terms): N is the observation dim. of the data, and by ‘hyperbolic AL-RNN’ we mean an AL-RNN which is hyperbolic in each of its linear subregions, i.e., such that the transition matrices defined in eq. 2 have no eigenvalues on the unit circle. Will be clarified.

W8 (use of acronyms): Agreed, we will remove all acronyms which are not standard or are only rarely used.

Questions

Q1: Please note that our goal was not to introduce a novel SOTA method for DSR, but rather to introduce an approach for retrieving topologically minimal and symbolically interpretable representations from data. Still, one may ask whether the AL-RNN is at least on par with other SOTA methods for DSR. In Table R1 in the rebuttal PDF we answer this Q. We also included human EEG data & Lorenz-96 as further benchmarks. As can be seen, the AL-RNN even outperforms most current SOTA methods, which may be due to its simple design making training much more stable and robust (cf. also new Figs. R1-R3).

Q2: The number of linear units makes no difference to interpretability, since it does not affect the total number of linear subregions (hence neither the symbolic encoding nor the computation of fixed points or cycles). As shown in Fig. 9, more linear units can, however, still improve performance up to a certain level.

Q3: See W2 above.

Q4: See response to W1 above. We also would like to note that esp. in math. chaos theory probabilistic approaches to deterministic systems are indeed commonplace, e.g. in the def. of invariant measures (which are probability ms.) and ergodic theory (see, e.g., textbooks by Katok & Hasselblatt 1995, Alligood et al. 1996).

We thank the referee for the supportive and positive feedback, we are happy to hear the referee liked our work!

Weaknesses

In a sense this is the first study of its kind. We are not aware of any other work in the DSR field (and beyond) making this link to symbolic dynamics, and attempting to reduce model complexity in a way that allows for easy translations into topologically minimal representations. This approach to model interpretability is the major contribution of this work, so it is hard to compare to other methods.

However, one may still ask whether AL-RNNs can at least compete with other DSR methods in terms of the quality of the DS reconstructions they achieve (our examples demonstrate that they are good enough even with a very low number of linear subregions, which in itself is a major advantage for obtaining mechanistic insight into dynamics). New Table R1 which we now added to the rebuttal PDF confirms AL-RNNs are at least on par with - in fact even outperform - most other SOTA methods (which lack our method’s interpretability). We also added new benchmarks to this Table.

Questions

Q1: $\phi$ is the ReLU nonlinearity, sorry for this oversight. We will clarify this in the updated manuscript.

Q2: The ReLU ($\phi$) was absorbed into the $D_{\Omega(t)}$ matrix: Note that we can rewrite eq. 1 equivalently into eq. 2 by placing a 1 on the diagonal of $D_{\Omega(t)}$ for each state $z_{p,t}>0$ at time $t$, and a 0 otherwise (just by definition of the ReLU as $\max[0,z_{p,t}]$). This matrix depends on time $t$, however, because it is determined by the values of all the states at that time. $\Omega(t)$ denotes the subset of states for which we have $z_{p,t}>0$ at time $t$. This should also answer the referee’s last question: Since each entry on the diagonal of the $M \times M$ matrix $D_{\Omega(t)}$ can be either 0 or 1 at any time $t$, we have a total of $2^M$ possible configurations for $D_{\Omega}$. This may admittedly have been a bit hard to follow, especially as we didn’t make clear the meaning of $\phi$ in the equation above! We will clarify this whole section accordingly.

Q3: Yes, its interpretation is a bit similar to a decimal point (but note it’s not a decimal system of course): The symbolic sequences correspond to theoretically infinite length trajectories, and the point “separates past from future”, i.e. indicates which symbol corresponds to the current time $t$ along the trajectory. We will clarify this type of notation in the revision (it's indeed standard in symbolic dynamics), valid point, also for the parts above, thank you for your feedback on this!

We thank the referee for the enthusiastic support and appreciation of our work!

Weaknesses

W1 (consistency of fits/ model recovery): One crucial advantage of AL-RNNs is that they indeed consistently deliver the same model over many training repetitions. The errors in Figs. 5d-f are in fact very small: To put these numbers into context, we now normalized all 3 graphs by proper reference values (Fig. R1 in rebuttal PDF), with normalized numbers substantially below 1. Moreover, we now compared the consistency across training runs to those obtained with vanilla PLRNNs by evaluating the agreement in trajectory point distributions across linear subregions, revealing a much higher agreement among AL-RNN solutions (Fig. R2). The parsimony of AL-RNNs compared to other models leaves them with much less ‘wiggle room’ in finding different solutions, as can also be appreciated from the direct comparison in Fig. R3 (see also Figs. 20 and 21). Finally, as suggested, we also did model recovery experiments, finding that the recovered solutions are virtually identical across repetitions (original: $D_{stsp}=3.14$, $D_H=0.28$; recovered: $D_{stsp}=3.38+/-0.18$, $D_H=0.28+/-0.03$; 3 linear subregions in all cases).

W2 (degeneracies): While these models indeed do have the capacity to capture unobserved dimensions in their latent space (e.g. Brenner et al. 2024), we usually still would use delay embedding for this (as for the ECG). This, both, eases the training process itself and enables evaluation of agreement in attractor geometry. While it is possible to harvest the symbolic representation to measure the similarity in 2 reconstructions (overlap in symbolic codes and their graphs, e.g. agreement in adjacency matrices), for empirical data the topologically minimal representation is usually unknown, of course, and is precisely what we would like to infer through the AL-RNN (see Discussion). Hence, empirically, measures like $D_{stsp}$ or $D_H$ remain the methods of choice for initially evaluating the quality of DSR in delay-embedding spaces. But we checked this idea now for evaluating the consistency across training runs, and the symbolic graphs in fact remain identical across different runs.

Questions

Q1: Absolutely, thanks for pointing out this oversight!

Q2: Yes, true, the full $\mathbf{W}$ matrix is kind of a ‘historical quirk’ from previous formulations of the PLRNN and the resp. codebase. For model training and performance it makes no difference ($D_{stsp}$: $t(19)=-1.58, p=0.13$, $D_{H}$: $t(19)=0.68, p=0.51$), but for parsimony the diagonal of $\mathbf{W}$ should be removed for the linear units. We will comment on this in the revision.

Q3: Yes, thanks for catching!

Q4: With “hyperbolic AL-RNN” we mean the system is hyperbolic in each of its linear subregions, implying that none of the Jacobian matrices $\mathbf{A}+\mathbf{W} \mathbf{D}_{\Omega}$ has eigenvalues exactly on the unit circle (a measure-0 set in parameter space). Will be clarified!

Q5: Thanks for pointing out this source of misunderstandings, we will rephrase the theorems for clarity. The key is the term “corresponding” by which we meant the mapping from trajectories of the original system onto symbolic seq.. We will precisely define this mapping and the term “corresponding” in our rev..

Re Th.2: The notation here is correct, this is the standard def. for a p-cycle (otherwise it would just read $\mathbf{z}_k=\mathbf{z}_k$, but we wish to indicate the map returns after p iterations).

Q6: Some of the wiggling up and down of some of the curves after the initial kink (which we are looking for) is likely to be just noise across different training runs with different minima achieved. Also note that we are not directly optimizing for $D_{stsp}$ and $D_{H}$, so while in STF based algorithms MSE remains a good proxy in general (e.g. Hess et al. 2023), there is no guarantee the relation is strictly monotonic. However, the hump in the figure for the Lorenz-63 is a bit more suspicious, and so we dug a bit deeper: It seems that in this case the first minimum at P=2 indeed indicates the topologically minimal representation (see also Fig. R6), and that when surpassing this optimal point in fact performance first decays again. It is thus more a feature rather than a bug, and would make identifying the optimal representation even easier.

Q7: For the ECG data, the quantitative results were not included in Fig. 7 but in Fig. 3, so perhaps easy to miss, will be made clearer. For the fMRI data, as requested we now have produced the same type of graph, included as Fig. R4 in the rebuttal PDF. That an additional categ. decoder significantly improves DSR quality on these same fMRI data has been shown before in (Kramer et al. 2022; Brenner et al. 2024).

Limitations

Since the number of subregions grows exponentially with $P$ and $M$, $P<<M$ eases model analysis profoundly. One major purpose of DSR in our minds is to provide mechanistic insight into system dynamics (rather than just forecasting). The thorough understanding of the dynamical mechanisms of chaos in Fig.5, e.g., is only possible because we have less than a handful of linear subregions (likewise for the empirical examples). Besides this important visualization aspect facilitating human interpretation, there are also clear numerical benefits: For instance, in a brute-force search algorithm for fixed points and cycles, the number of iterations we would need to find such objects also increases exponentially with $M$. Hence reducing this to $P<<M$ is always a huge benefit that enables to dig much deeper into model mechanisms. Of course it is difficult to put an exact threshold on this. However, that we were able to capture even rather complex ECG and fMRI dynamics with just a few linear subregions ($\leq 8$) we find encouraging. Whether this will always or mostly be the case in empirical settings, is an interesting and open question, which will be discussed.

Dear Referees,

The discussion period is coming to a close, and we wondered whether we could satisfactorily address your points, or whether there are any issues remaining that may need further clarification? Thank you very much again for reviewing our work, and for your supportive and positive feedback so far!

Kindly, Authors