We thank the reviewer for their constructive feedback. Many good points were raised that we address below.

For clarity, we numbered the questions to make them easier to follow. Please note that we moved Questions 2 and 4 to the end, where we also address points raised as weaknesses.

Q1: Equal compute and sequence lengths in the experiments

Yes, considering equal compute is becoming standard practice for fair comparison. In our experiments, the sequence lengths were determined by the characteristics of the datasets themselves, using their original length.

Q3: Synaptic parameters are constant matrices, how general...

What we meant is that the synaptic parameters, the learned values $a_{ji}, b_{ji}$, $g_{ji}$ and $k_{ji}$ are not state- or input-dependent, they remain constant after training and define the fixed connectivity between neurons (neuron $j$ and $i$). This is generally true for any NSSM formulation based on this structure.

Q5: Euler integration and number of iterations

We follow the integration scheme from LRCs, which requires only a single Euler integration step. As a result, the system can be solved immediately without the need for multiple Euler iterations.

Q6: Matrix diagonalisation, joint diagonalizer

Thank you for bringing this up, we will clarify it in the revision. We explicitly constrain the model to learn diagonal matrices (which can be viewed as vectors). Therefore, the need for a joint diagonalizer does not arise, as the matrices are inherently diagonal by design.

Weakness 1 and Q4: Motivation for LRC and model with constant capacitance

We considered models of non-spiking neurons, such as those based on electrical and chemical synapses. Electrical synapses can be modeled using continuous-time recurrent neural networks (CT-RNNs), but they tend to show degraded performance compared to models based on chemical synapses. In such models, part of the transition matrix is fixed, as in many state-space models (SSMs) before, and lacks input or state dependency, which offers nothing novel in this regard.

Instead, we focused on models of chemical synapses, the so called liquid neural networks, where the term liquid refers to the input- and state-dependent dynamics of the transition matrix. These dynamics are non-linear. Within this family, Liquid Time-Constant Networks (LTCs) do not have saturation (or gating) terms in their dynamics, which is desirable for achieving stable behavior. Fortunately, their saturated (gated) counterpart, Saturated LTCs (STCs), does include such terms.

We modified the model of STCs the same way as LRCs, achieving a diagonal Jacobian matrix (without the elastance term indicated in orange in the equations) and ended up StcSSMs. This model has constant capacitance. Models coming from modeling chemical synapses assume this.

We included a comparison of this model against LrcSSM across the datasets, using a fixed setup of 6 layers, 64-dimensional encoded input and 64-dimensional states. We found that LrcSSM outperform StcSSM, highlighting the need for the non-linear membrane capacitance - pointed out by the reviewer. We will include this observation in our ablation section.

Weakness 2 and Q2: Comparing diagonal and full matrix results

This is a valuable point that we wanted to address more thoroughly. To this end, we evaluated the LrcSSM model with a dense $A$ matrix as well, to compare performance and assess any potential trade-offs. We used the same architectural design as before to ensure a fair comparison.

These results show that empirically, we do not lose the expressive capabilities of LRCs by constraining $A$ to be diagonal. We hypothesize that this is because the strict structure encourages the model to learn more efficient representations, while also enabling exact (rather than approximate) computations due to the design.

This is another important part pointed out by the reviewer, which we will include in the paper.

Weakness 3: Computational costs

Let $T$ denote the input sequence length and $D$ the state dimension (i.e. the dimension of the diagonal of the matrix). Traditional sequential methods are difficult to parallelize and require $O(D)$ memory and $O(T D^2)$ computational work. The ELK technique [1] addresses this by enabling efficient parallelization through diagonal Jacobian computation, reducing both memory and computational complexity to $O(T D)$. Our approach achieves the same complexity: $O(T D)$ for both memory and computation. Since our model learns inherently diagonal matrices (i.e., vectors), Jacobians remain diagonal throughout, enabling efficient and scalable computation without additional overhead, making previously quasi-approximations exact.

Weakness 4: Large average variance in the datasets

We reported the variances in the submission to provide a more complete and transparent view of model performance, but actually, this is something omitted in all prior work using these datasets (only reporting mean values). The high variance is consistent across all models and reflects only the inherent variability in achievable accuracy on these datasets, which shifts the attention away from the more meaningful mean values. We decided now to only keep them for each dataset separately, because it is much more relevant there and shows the reliability of the models across 5 seeds.

References:

[1] Xavier Gonzalez, Andrew Warrington, Jimmy Smith, and Scott Linderman. Towards scalable and stable parallelization of nonlinear rnns. Advances in Neural Information Processing Systems, 37:5817–5849, 2024

We would like to sincerely thank the reviewer for their thoughtful and detailed feedback. Many of the points raised were insightful and important, which we addressed in the following way:

Weakness 1 and Question on practical efficiency, runtime

Here we report the wall-clock times for the best models for each dataset, in seconds for 1000 training steps:

Number of iterations needed for convergence for LrcSSM on average per dataset in the same order: 3.4, 2.1, 4.4, 4.6, 3.0, 3.8.

Results for the other models are from [1]. As one can see, we do not lose significant runtime due to the extra iterations that are needed for the states to converge.

By efficient sequence modeling, we mean that we have taken a step toward scaling nonlinear RNNs to nonlinear state-space models (NSSMs), which can process long sequences in parallel. As the reviewer rightly pointed out, NSSMs require $k$ iterations to converge to the desired state dynamics, unlike linear SSMs which don't require this. However, compared to traditional nonlinear RNNs, which require sequential computation over the full sequence length $L$, our method overcomes this limitation through parallelization, at the cost of a small number of additional iterations.

In practice, this is highly favorable: the number of iterations $k$ is typically based on experiments very small (between 2 and 4), while sequence lengths $L$ range from 400 up to 18,000. We observed that the number of iterations $k$ is independent of both the sequence length $L$ and the state dimension $D$. Since $k$ is effectively a constant, we believe it is reasonable to omit it from the $O$-notation.

That said, we agree with the reviewer that this should be clearly discussed in the limitations or discussion section of the paper to make our claims more transparent.

Weakness 2. and 3. Performance gains and broader impact

The reviewer raised a critical point that we aim to address. Although in the submission we only demonstrated how to construct the model from LRC to LrcSSM, this approach can, in theory, be applied to any nonlinear RNNs of interest to obtain its SSM counterpart.

Sketch of the generalized method:

Any RNN can be reformulated by considering its underlying differential equations using an explicit first-order Euler integration scheme with a unit time step. We assume the non-linear RNN has the following general form:

$x_t = q_t(x_{t-1},u_t) x_{t-1} + s_t(x_{t-1},u_t)$,

where $q_t=\sigma(f_q(x_{t-1},u_t))$ and $s_t=\sigma(f_s(x_{t-1},u_t))$ are non-linear potentially state- and input-dependent gated functions, which we aim to express it as:

$\dot{x} = A(x,u)x+b(x,u)$.

To see this, consider:

$\dot{x} = (q(x,u)-1)x+s(x,u)$,

which yields (using Euler integration with $\Delta t=1.0$) the original equation:

$x_t = x_{t-1} + \dot{x} \Delta t = x_{t-1} + (q_t (x_{t-1},u_t) -1) x_{t-1} + s_t(x_{t-1},u_t) = q_t (x_{t-1},u_t) x_{t-1} + s_t(x_{t-1},u_t)$.

Thus, in this general formulation: $A(x,u)=(q(x,u)-1)$, $b(x,u) = s(x,u)$.

To obtain the desired efficient representation, we define:

$A(x,u) = \mathrm{diag}([\sigma(f_{q1}(x_1,u)) - 1, \sigma(f_{qi}(x_i,u)) - 1, \ldots, \sigma(f_{qm}(x_m,u)) - 1])$

$b(x,u) = [\sigma(f_{s1}(x_1,u)), \sigma(f_{si}(x_i,u)), ..., \sigma(f_{sm}(x_m,u))]$.

Example: GRU:

The GRU equations are:

$z_t = \sigma\left(W_z \left[ x_{t-1}, u_t \right] + b_z \right) = \sigma(f_z)$

$r_t = \sigma\left(W_r \left[ x_{t-1}, u_t \right] + b_r \right)= \sigma(f_r)$

$c_t = \tanh\left(W_h \left[ r_t x_{t-1}, u_t \right] + b_h \right) = \tau(f_p)$

$x_t = (1 - z_t) x_{t-1} + z_t c_t$

Rewriting the update using the differential form:

$\dot{x} = z(-x + c) = -z x + z c$, to check, we get back to the original form (using Euler integration with $\Delta t=1.0$)

$x_t = x_{t-1} + \dot{x} \Delta t = x_{t-1} + -z_t x_{t-1} + z_t c_t = (1 - z_t) x_{t-1} + z_t c_t$.

To express this in the form $\dot{x} = A(x,u)x+b(x,u)$, we set: $\dot{x} = A(x,u)x+b(x,u) =-z x + z c$, thus $A(x,u) = -z$, and $b(x,u)=z c$.

Using our proposed formulation for constructing $A(x,u)$ and $b(x,u)$:

$A(x,u) = \mathrm{diag}([-\sigma(f_{z1}(x_1,u)), -\sigma(f_{zi}(x_i,u)), ..., -\sigma(f_{zm}(x_m,u))]$,

$b(x,u) = [\sigma(f_{z1}(x_1,u))\tau(f_{p1}(x_1,u)), \sigma(f^*_{zi}(x_i,u))\tau(f_{pi}(x_i,u)), ..., \sigma(f_{zm}(x_m,u))\tau(f_{pm}(x_m,u))]$.

This reformulation results in a diagonal Jacobian matrix, enabling more efficient and exact (rather than quasi) computations for ELK.

In the line of this work, we extend our method to build up other nonlinear SSMs from Minimal Gated Units (MGU) [2], Gated Recurrent Units (GRU) [3], and LSTM [4], which we refer to as MguSSM, GruSSM and LstmSSM, respectively. We report their performance in the following table:

Test accuracy across these non-linear models:

Here, similar to our ablation studies, we used 6 layers with a 64-dimensional input encoding and 64-dimensional state. The results show that LrcSSM outperforms the other NSSMs built from MGU, GRU and LSTM using the same architecture.

Note: The stacked blocks (2,4 or 6 are part of the hyperparameter search grid) and the normalization blocks are also present in other models reported too.

Overall, with this, we show the performance gains in the same architecture across different non-linear models and we also show the broader impact of our paper, as this diagonalization technique can be applied to many RNNs.

Weakness 4: Model simplification

Again, this is a valuable point that we wanted to address more thoroughly. To this end, we evaluated the LrcSSM model with a dense $A$ matrix as well, to compare performance and assess any potential trade-offs. We used the same architectural design as before to ensure a fair comparison.

These results show that empirically, we do not lose the expressive capabilities of LRCs by constraining $A$ to be diagonal. We hypothesize that this is because the strict structure encourages the model to learn more efficient representations, while also enabling exact (rather than approximate) computations due to the design.

References:

[1] T. Konstantin Rusch and Daniela Rus. "Oscillatory State-Space Models". Thirteenth International Conference on Learning Representations (2025).

[2] Zhou, Guoxiang, Jianxin Wu, Chen-Lin Zhang and Zhi-Hua Zhou. “Minimal gated unit for recurrent neural networks.” International Journal of Automation and Computing 13 (2016): 226 - 234.

[3] Cho, Kyunghyun, Bart van Merrienboer, Çaglar Gülçehre, Dzmitry Bahdanau, Fethi Bougares, Holger Schwenk and Yoshua Bengio. “Learning Phrase Representations using RNN Encoder–Decoder for Statistical Machine Translation.” Conference on Empirical Methods in Natural Language Processing (2014).

[4] S. Hochreiter and J. Schmidhuber, "Long Short-Term Memory," in Neural Computation, vol. 9, no. 8, pp. 1735-1780, 15 Nov. 1997, doi: 10.1162/neco.1997.9.8.1735.

We thank the reviewer for their thoughtful feedback, which has helped us improve our work.

See our responses below:

appendix. Please move the definition forward or otherwise drop it.

Noted. Thank you.

Q2) Efficient computation of the diagonal of the jacobian

No, unfortunately, it is not straightforward.

Q3) Run times

Here are our wall-clock times for the best models for each dataset, in seconds for 1000 training steps:

Results for the other models are from [1]. We would like to point out that our efficiency claims are not against other linear models, because they are really efficient. We aimed to bring up non-linear models to the same/comparable efficiency. Iterations needed for convergence for LrcSSM on average per dataset in the same order: 3.4, 2.1, 4.4, 4.6, 3.0, 3.8. To make it explicit, this number is 1.0 for all other models across all datasets. As one can see, we do not gain significant runtime due to the extra iterations that are needed for the states to converge.

Q4) Discrete sequences

We are not sure if we get the question. The datasets we are using are discrete sequences.

Q5) Transformers

Good point. Before, we missed a recent publication [2] from last NeurIPS, where Transformers and Rough Transformers (RFormer) were evaluated on the same datasets. We report their results [2] here, which show that LrcSSM achieves better performance than these Transformer-based models. We will incorporate these new model comparisons into our paper.

Q6) Analogous discretizations for STC

Indeed, the same discretization approach can be applied to STCs. One might naturally ask why a model similar to LRC was not constructed for STCs. To address this, we present an alternative formulation of StcSSM, following the same steps used to achieve a diagonal Jacobian matrix (without the elastance term indicated in orange), in order to compare the performance of these two bio-inspired models.

Test accuracy across the same model architecture of 6 layers, 64-dimensional input encoding and 64-dimensional states:

These results show that our originally proposed LrcSSM model performs better than StcSSM, and this performance gain comes from the extra expressiveness of the underlying LRC model, which has an input- and state-dependent membrane capacitance in its bio-inspired model, which is a constant term in the case of the STC model.

Q7) PPG-Dalia

Thank you for pointing this out. We ran this additional experiment on the PPG-DaLiA dataset [4], and report the following results. As before, baselines are taken from [1]:

We ran the hyperparameter tuning following the same protocol as before. LrcSSM achieves performance comparable to other models, demonstrating its competitiveness in this benchmark.

Q8) Moving the pre-normalization block

Our goal was to apply a normalization layer before each SSM block. Moving it outside would only normalize the original input to the first layer, rather than each layer's input. This approach is consistent with other models, where normalization is similarly applied before the SSM layer.

Q9) Number of layers

We used the same hyperparamater tuning grid as it was used for the other models in [1],[3] to have full comparability. This grid search includes layers={2,4,6}.

Q10) Loss function for training

We used softmax cross-entropy for the classification tasks.

References:

[1] T. Konstantin Rusch and Daniela Rus. "Oscillatory State-Space Models". Thirteenth International Conference on Learning Representations (2025).

[2] Moreno-Pino, Fernando, Alvaro Arroyo, Harrison Waldon, Xiaowen Dong and Álvaro Cartea. “Rough Transformers: Lightweight and Continuous Time Series Modelling through Signature Patching.” Neural Information Processing Systems (2024).

[3] Benjamin Walker, Andrew Donald McLeod, Tiexin Qin, Yichuan Cheng, Haoliang Li, and Terry Lyons. "Log neural controlled differential equations: The lie brackets make a difference." Forty-first International Conference on Machine Learning (2024).

[4] Reiss, Attila, Ina Indlekofer, Philip Schmidt, and Kristof Van Laerhoven. 2019. "Deep PPG: Large-Scale Heart Rate Estimation with Convolutional Neural Networks" Sensors 19, no. 14: 3079. doi: 0.3390/s19143079

We thank the reviewer for their comments. In the responses below, we aim to clarify the concerns raised.

1. State-transition matrix

The reviewer is correct, this term usually refers to linear systems. To incorporate this feedback, we will explicitly highlight that we mean non-linear state transition matrix in the text when we talk about non-linear SSMs.

2. Comparison with nonlinear RNNs

To the best of our knowledge, there are currently no other NSSMs. However, the reviewer raises an excellent point. To enable comparison with other types of non-linear models, we formulated a general framework that allows scaling any non-linear RNN in the same way as we showed by constructing a diagonal Jacobian.

Using this approach, we extend our method to build NSSMs from Minimal Gated Units (MGU) [1], Gated Recurrent Units (GRU) [2], and LSTM [3], which we refer to as MguSSM, GruSSM and LstmSSM, respectively. We report their performance in the following table:

Test accuracy across non-linear models:

Here, similar to our ablation studies, we used 6 layers with a 64-dimensional input encoding and 64-dimensional state. The results show that LrcSSM outperforms the other NSSMs built from MGU, GRU and LSTM.

3. LinOSS relation to LSSMs

Good point. LinOSS is a special type of LSSM, which is based on forced linear second-order ODEs. Meanwhile, the other LSSMs have first-order linear dynamics.

4. Categorization

Thank you for pointing this out. We will go through the paper and clarify the terminology to ensure consistency.

For this and the previous point, we plan to include visualizations to illustrate how these models form categories. However, due to rebuttal restrictions requiring text-only responses, we are unable to provide them at the moment.

In our paper, we refer to NSSMs (non-linear state-space models) as those with non-linear transition matrices $A$, while LSSMs (linear state-space models) have linear transition matrices, which makes them inherently parallelizable. Our work focuses on the NSSM category, which has traditionally been avoided due to the belief that such models could not be solved in parallel and required sequential processing over the sequence length.

However, recent work [4] introduced the ELK algorithm as a way to parallelize these models efficiently and approximately. Building on this, we propose a model architecture with a diagonal Jacobian structure, enabling exact (not quasi) parallel computation for NSSMs.

Specifically, we explore LrcSSM, a model that leverages this structure alongside the new scaling technique of ELK. As shown in our experiments, LrcSSM performs strongly compared to other NSSMs (see above) and even many leading LSSMs (see the paper).

The key advantage of LrcSSM is that it allows for non-linear dynamics within a parallelizable framework, something not achieved before. While this introduces a small number of additional iterations per step, it does not significantly increase runtime and yields performance competitive with, or superior to, strong baselines like LRU (which argued non-linearity was unnecessary), S5, Mamba, and S6.

While LinOSS remains a strong baseline due to its rich second-order dynamics, which we plan to explore further in future work, our LrcSSM model introduces a novel direction by incorporating a non-linear, input- and state-dependent transition matrix $A$, which we believe adds a meaningful contribution to the field and opens a new category of exploring NSSMs.

References:

[1] Zhou, Guoxiang, Jianxin Wu, Chen-Lin Zhang and Zhi-Hua Zhou. “Minimal gated unit for recurrent neural networks.” International Journal of Automation and Computing 13 (2016): 226 - 234.

[2] Cho, Kyunghyun, Bart van Merrienboer, Çaglar Gülçehre, Dzmitry Bahdanau, Fethi Bougares, Holger Schwenk and Yoshua Bengio. “Learning Phrase Representations using RNN Encoder–Decoder for Statistical Machine Translation.” Conference on Empirical Methods in Natural Language Processing (2014).

[3] S. Hochreiter and J. Schmidhuber, "Long Short-Term Memory," in Neural Computation, vol. 9, no. 8, pp. 1735-1780, 15 Nov. 1997, doi: 10.1162/neco.1997.9.8.1735.

[4] Xavier Gonzalez, Andrew Warrington, Jimmy Smith, and Scott Linderman. Towards scalable and stable parallelization of nonlinear rnns. Advances in Neural Information Processing Systems, 37:5817–5849, 2024

We thank the reviewer for their highly constructive feedback. We appreciate it a lot and will incorporate it into our revision. Here are our replies:

1. LRA benchmarks and comparing linear vs. non-linear models

Thank you for the suggestion. We would like to point out that, particularly in recent work, researchers have increasingly shifted from the LRA benchmarks to the UEA-MTSCA archive. Our decision to focus on UEA-MTSCA aligns with this trend, as reflected in several recent top publications [1]-[3]. To further support this direction, we also cite here additional examples published on arXiv after the NeurIPS submission, in May/June [4]-[8].

That said, we completely agree that including LRA would provide valuable context and comparison. We are very interested in exploring how our non-linear LrcSSM performs on LRA, especially in relation to recent linear models such as LinOSS, which also don't report performance on these benchmarks. However, we believe such a broader investigation would be more out of the scope of the current paper.

Nonetheless, one can address this point in future work by exploring which types of datasets benefit more from Transformer-like, linear-SSM-like, or non-linear-SSM-like architectures.

For now, as also pointed out by other reviewers, we ran additional experiments on the PPG-DaLiA dataset [9], and can report the following results (baselines are taken from [3]):

We ran the hyperparameter tuning following the same protocol as before. LrcSSM achieves performance comparable to other models, demonstrating its competitiveness in this benchmark.

2. Compute-optimal scaling law regime

Thank you for raising this point. We agree that claims in the abstract should be clear and well-supported. Our intention was not to present a new empirical measurement of the scaling exponent $\beta$ for LrcSSM, but rather to argue based on properties of LrcSSM that it is expected to follow the same compute-optimal scaling regime. We will revise the abstract to clarify this and note that a full scaling study for LrcSSM is a valuable direction for future work.

3. Table 7 and input dependency

In Table 7, we report the mean and standard deviation of accuracy across the 5 different dataset splits. We observed that performance can vary substantially from split to split, which naturally leads to higher standard deviations.

The deviations appear larger in Table 7 compared to earlier tables because these results are not fully tuned over the entire hyperparameter grid. Instead, we use a standardised configuration across ablations (6 layers, 64-dimensional input encoding, and 64 hidden states) to ensure comparability across the model variants.

Regarding statistical significance: yes, the standard deviations are relatively large, which makes the differences between variants less statistically robust. However, this is a common observation in Table 3 as well, where many baseline models show similar variance. This stems from the inherent variability in achievable accuracies across the datasets.

Despite the high variance, we focused on identifying configurations with the highest mean accuracy. Your interpretation is correct: the input-dependent variant, which converges in 1 iteration, but performs slightly worse than the state-dependent version. Given this trade-off, we choose the better performing model, as it does not impose a significant computational overhead, please check our reported results in the next point.

4. Wall-clock times

Wall clock time (in seconds) for 1000 training steps of the best performing (tuned) models for each dataset, same as it was reported in [3].

Iterations needed for convergence for LrcSSM on average per dataset in the same order: 3.4, 2.1, 4.4, 4.6, 3.0, 3.8. To make it explicit, this number is 1.0 for all other models across all datasets. As one can see, we do not lose significant runtime due to the extra iterations that are needed for the states to converge.

5. Liquid-S4 baseline

Instead of reporting results for Liquid-S4 directly, we include S6, which builds on the same core idea, which is an input-dependent $A$ matrix, but incorporates parallel scan operations for improved efficiency. In this sense, S6 can be viewed as a more efficient version of Liquid-S4, and we include it as a representative of that model class.

6. Mamba results

We would like to clarify that Mamba is an overall architecture design that uses S6 (or other SSM layers) as its core building block. Therefore, the more appropriate comparison is between S6 and LrcSSM. For the Eigenworms dataset, we report an accuracy of 85.0% for S6, which aligns closely with our results of the input-dependent version of LrcSSM in Table 7.

The difference in performance between S6 and the input-dependent-only LrcSSM likely stems from the linear nature of S6 vs the non-linear model of LrcSSM. This distinction can explain the performance gap observed despite both being input-dependent systems.

7. ELK

Thank you for pointing this out. We use (quasi)-ELK, which in our case is equivalent to ELK because our Jacobian matrices are diagonal. Thus, it is not an approximation but an exact computation for our model. We used the default parameter setting, with lambda=1/10^8. We will include this information explicitly in the paper.

8. Comparison against Transformers

Good point. Before, we missed this recent publication from last NeurIPS, which introduced a new Transformer model called Rough Transformer (RFormer) [1] and evaluated it on the same benchmarks. LrcSSM achieves better performance compared to the Transformer and competitive/better than the Rough Transformer baselines across the evaluated datasets. We will incorporate these results into our paper.

Limitations:

• What we meant by "parallel with a single prefix-scan" is that each iteration of our iterative algorithm processes the entire sequence of length $T$ simultaneously, rather than processing the sequence elements one at a time as in non-linear RNNs. We'll revise our text to make this clear.
• See point 4.

References:

[1] Moreno-Pino, Fernando, Alvaro Arroyo, Harrison Waldon, Xiaowen Dong and Álvaro Cartea. “Rough Transformers: Lightweight and Continuous Time Series Modelling through Signature Patching.” Neural Information Processing Systems (2024).

[2] Benjamin Walker, Andrew Donald McLeod, Tiexin Qin, Yichuan Cheng, Haoliang Li, and Terry Lyons. "Log neural controlled differential equations: The lie brackets make a difference." Forty-first International Conference on Machine Learning (2024).

[3] T. Konstantin Rusch and Daniela Rus. "Oscillatory State-Space Models". Thirteenth International Conference on Learning Representations (2025).

[4] Walker, Benjamin, Lingyi Yang, Nicola Muca Cirone, Cristopher Salvi, and Terry Lyons. "Structured Linear CDEs: Maximally Expressive and Parallel-in-Time Sequence Models." arXiv preprint arXiv:2505.17761 (2025).

[5] Boyer, Jared, T. Konstantin Rusch, and Daniela Rus. "Learning to Dissipate Energy in Oscillatory State-Space Models." arXiv preprint arXiv:2505.12171 (2025).

[6] Pourcel, Guillaume, and Maxence Ernoult. "Learning long range dependencies through time reversal symmetry breaking." arXiv preprint arXiv:2506.05259 (2025).

[7] Nzoyem, Roussel Desmond, Nawid Keshtmand, Idriss Tsayem, David AW Barton, and Tom Deakin. "Weight-Space Linear Recurrent Neural Networks." arXiv preprint arXiv:2506.01153 (2025).

[8] Karuvally, Arjun, Franz Nowak, Anderson T. Keller, Carmen Amo Alonso, Terrence J. Sejnowski, and Hava T. Siegelmann. "Bridging Expressivity and Scalability with Adaptive Unitary SSMs." arXiv preprint arXiv:2507.05238 (2025).

[9] Reiss, Attila, Ina Indlekofer, Philip Schmidt, and Kristof Van Laerhoven. 2019. "Deep PPG: Large-Scale Heart Rate Estimation with Convolutional Neural Networks" Sensors 19, no. 14: 3079. doi: 0.3390/s19143079

I congratulate the authors on their great work and well-considered reviews. I am particularly heartened by their strengthened experimental evaluation. In particular, their wall clock time results are extremely exciting, indicating the speed of LrcSSM. Also exciting are the fast convergence rates of quasi-DEER. These results should definitely be included in the final paper.

Can the authors provide richer distributional statistics for the iterations needed for convergence for LrcSSM? That is, instead of just average, can they provide min, Q1, median, Q3, and max? (should be doable with one line of code, df.summary()). Doing so would give a better perspective as to possible worst case behavior during training and pathologies.

On the balance, this innovative paper has opened up new perspectives for me on the parallelization of nonlinear recurrent architectures. It has already begun to affect how I think about future research directions. Papers like these should definitely be discussed at conference like NeurIPS. I have strong conviction and confidence that this paper should be accepted for publication. I have raised my score to a "5."

I include some additional feedback aimed at further strengthening the paper. I have three minor points from this review, and a major point from Reviewer YHNE.

Table 7

I would recommend reporting the standard error instead of the standard deviation. This way, the average would illustrate a lower standard error (as it must), instead of a standard deviation that is larger than any particular experiment. Doing so would be more statistically accurate, and also would strengthen your story.

Mamba on Eigenworms

I am still underwhelmed by your discussion of the performance of Mamba on eigenworms. S6 is a part of Mamba. Presumably the authors of Mamba would not have made additions to from S6 to Mamba if they did not believe that improved performance. Can you please comment on the differences between S6 and Mamba, and why you believe that S6 is outperforming Mamba so dramatically here?

quasi-DEER or quasi-ELK

I looked through your very well organized code. It appears to me that quasi-DEER, and not quasi-ELK, is being called. Can you point me to the line in your zip code that calls quasi-ELK in your model training architecture?

Point from reviewer YHNE: what is a NSSM?

I am confused by your discussion of NSSM.

In our paper, we refer to NSSMs (non-linear state-space models) as those with non-linear transition matrices, while LSSMs (linear state-space models) have linear transition matrices, which makes them inherently parallelizable.

Both LrcSSM and Mamba are switching linear dynamical systems. Mamba's transition matrix is also a nonlinear function of its inputs. The difference between LrcSSM and Mamba is what their transition matrices are functions of. Mamba is only input dependent, which is why it can be parallelized in a single parallel scan. LrcSSM is state-dependent, which is why it requires quasi-DEER for parallelization.

I also do not think that reserving the term "SSM" for (switching) linear dynamical systems is a good decision. A "state space model" has classically referred to any stateful (Markovian) system, with arbitrary transition function. Under this more traditional definition, Reviewer YHNE is correct that xLSTM is a nonlinear state space model (because of the nonlinear dynamics in sLSTM). I agree with their recommendation that xLSTM should be added to the experiments (including wall clock time).

We thank the reviewer again for their thoroughness, including reading the other reviewers' comments, and for recognizing the value of our contribution by raising the score.

Below, we provide our responses:

Statistics for the iterations in LrcSSM

We will provide this in a boxplot format in the paper.

Table 7

Thank you for the suggestion. We agree that standard error is more appropriate, as it better reflects confidence in the mean than standard deviation.

Following your recommendation, we will revise Table 7 to report standard errors instead. For reference, the updated values are: 65.4 ± 8.02, 65.0 ± 8.04, 63.8 ± 8.37

We will also apply this adjustment consistently across all relevant tables in the paper.

Mamba on EigenWorms

We agree that the reviewer raises a valid point and appreciate the opportunity to further share our thoughts. Our understanding is that the additions introduced in Mamba were designed to improve performance on the specific datasets evaluated in the Mamba paper. However, in the context of the tasks evaluated here, particularly the Eigenworms dataset, these enhancements do not appear to offer the same benefit.

It is important to note that the results we report for both S6 and Mamba come from a different paper [1], not the original Mamba paper [2]. While we did not reproduce these results ourselves, we assume they are valid and reported in good faith. Nonetheless, the observed performance gap between S6 and Mamba on this dataset remains a task for future investigation.

quasi-DEER or quasi-ELK

Thank you for catching this, the reviewer is absolutely right. This is a mistake on our part. We used quasi-DEER, which is in fact a specialization of quasi-ELK when the lambda values are close to zero. We will correct the terminology in the paper to accurately reflect this, and we appreciate you pointing it out.

NSSM discussion, dependency and (non)-linearity

We agree with the reviewer on the key distinction regarding dependency: Mamba is input-dependent, while LrcSSM is both input- and state-dependent in $A$, this is well understood.

We would like to clarify some aspects regarding parallel computations.

Mamba uses linear projections [2] and includes exponential terms for discretization for $\bar A$, where the time-varying component comes from input dependency. The Mamba authors were able to construct the computations in an associative manner, with a customized parallel scan. However, this approach does not generalize to arbitrary non-linear projections in $A$. In the general case, methods like ELK are needed to achieve parallelization.

Overall, when we refer to non-linear SSMs, the non-linear nature is a core part of $A$.

Regarding xLSTM [3], including experiments on it is a valuable suggestion. However, due to time constraints, we are unable to provide these during the rebuttal period. We hypothesize that its performance is reflected in our LstmSSM ablation study, which we have included in the rebuttal.

[1] Benjamin Walker, Andrew Donald McLeod, Tiexin Qin, Yichuan Cheng, Haoliang Li, and Terry Lyons. "Log neural controlled differential equations: The lie brackets make a difference." Forty-first International Conference on Machine Learning (2024).

[2] Gu, Albert, and Tri Dao. "Mamba: Linear-time sequence modeling with selective state spaces." arXiv preprint arXiv:2312.00752 (2023).

Statistics for the iterations in LrcSSM

These are excellent! Well done! I agree their inclusion makes the paper stronger.

quasi-DEER or quasi-ELK

Yes, quasi-DEER is quasi-ELK with $\lambda=0$ (no regularization). However, you are right to directly use the quasi-DEER code for numerical and computational reasons.

When you can parallelize a linear dynamical system (LDS)

We appear to have a fundamental misunderstanding about which LDSs can be parallelized. Luckily, the results of the paper hold as is. However, I believe that clearing up this misunderstanding is important for contextualization of this paper.

I believe that any LDS with transition matrix NOT dependent on state can be parallelized. This includes an LDS with transition matrices that depend on the inputs in nonlinear ways. I hope that we can discuss this point from first principles.

Blelloch ('91) has written an manuscript on the parallel associative scan, also known as the prefix sum operator. This parallel scan is implement in jax as jax.lax.associative_scan and is used by your code. The parallel scan parallelizes matrix multiplication with a binary tree like approach.

I can write in pseudocode a way to parallelize an LDS with $A_t$ depending nonlinearly on inputs $u_t$.

Say $A_t = f(u_t)$, where $f$ is an arbitrary (possibly nonlinear) function.

1. As = jax.vmap(f)(us), where us are all the inputs (must be known in advance, otherwise we could not parallelize. This operation is embarrassingly parallel. As is of shape (T,D,D), and is the transition matrices (depending nonlinearly on u_t).
2. jax.lax.parallel.scan(matmul_op, As) gives the output in a parallel manner, running in $O(\log (T))$ time.

This fact that any input dependent switching LDS is parallelizable with one applicaition of a parallel associative scan is best explained in Feng et al, 2024, "Were RNNs All we needed." See their beautiful algorithm blocks for the minGRU and minLSTM and their excellent background section, particularly section 2.3

I hope that the authors can clear up this misunderstanding, because doing so is important for contextualizing the paper!