Oscillatory State-Space Models

## Abstract

Abstract We propose Linear Oscillatory State-Space models (LinOSS) for efficiently learning on long sequences. Inspired by cortical dynamics of biological neural networks, we base our proposed LinOSS model on a system of forced harmonic oscillators. A stable discretization, integrated over time using fast associative parallel scans, yields the proposed state-space model. We prove that LinOSS produces stable dynamics only requiring nonnegative diagonal state matrix. This is in stark contrast to many previous state-space models relying heavily on restrictive parameterizations. Moreover, we rigorously show that LinOSS is universal, i.e., it can approximate any continuous and causal operator mapping between time-varying functions, to desired accuracy. In addition, we show that an implicit-explicit discretization of LinOSS perfectly conserves the symmetry of time reversibility of the underlying dynamics. Together, these properties enable efficient modeling of long-range interactions, while ensuring stable and accurate long-horizon forecasting. Finally, our empirical results, spanning a wide range of time-series tasks from mid-range to very long-range classification and regression, as well as long-horizon forecasting, demonstrate that our proposed LinOSS model consistently outperforms state-of-the-art sequence models. Notably, LinOSS outperforms Mamba by nearly 2x and LRU by 2.5x on a sequence modeling task with sequences of length 50k.

## 1 Introduction

State-space models (Gu et al., ; Hasani et al., ; Smith et al., ; Orvieto et al., ) have recently emerged as a powerful tool for learning on long sequences. These models posses the statefullness and fast inference capabilities of Recurrent Neural Networks (RNNs) together with many of the benefits of Transformers (Vaswani, ; Devlin, ), such as efficient training and competitive performance on large-scale language and image modeling tasks. For these reasons, state-space models have been successfully implemented as foundation models, surpassing Transformer-based counterparts in several key modalities, including language, audio, and genomics (Gu & Dao, ).

Originally, state-space models have been introduced to modern sequence modelling by leveraging specific structures of the state matrix, i.e., normal plus low-rank HiPPO matrices (Gu et al., ; ), allowing to solve linear recurrences via a Fast Fourier Transform (FFT). This has since been simplified to only requiring diagonal state matrices (Gu et al., ; Smith et al., ; Orvieto et al., ) while still obtaining similar or even better performance. However, due to the linear nature of state-space models, the corresponding state matrices need to fulfill specific structural properties in order to learn long-range interactions and produce stable predictions. Consequently, these structural requirements heavily constrain the underlying latent feature space, potentially impairing the model’s expressive power.

In this article, we adopt a radically different approach by observing that forced harmonic oscillators, the basis of many systems in physics, biology, and engineering, can produce stable dynamics while at the same time seem to ensure expressive representations. Motivated by this, we propose to construct state-space models based on stable discretizations of forced linear second-order ordinary differential equations (ODEs) modelling oscillators. Our additional contributions are:

- •
we introduce implicit and implicit-explicit associative parallel scans ensuring fast training and inference.
- •
we show that our proposed state-space model yields stable dynamics and is able to learn long-range interactions only requiring nonnegative diagonal state matrix.
- •
we demonstrate that a symplectic discretization of our underlying oscillatory system conserves its symmetry of time reversibility.
- •
we rigorously prove that our proposed state-space model is a universal approximator of continuous and causal operators between time-series.
- •
we provide an extensive empirical evaluation of our model on a wide variety of sequential data sets with sequence lengths reaching up to 50k. Our results demonstrate that our model consistently outperforms or matches the performance of state-of-the-art state-space models, including Mamba, S4, S5, and LRU.

## 2 The proposed State-Space Model

Our proposed state-space model is based on the following system of forced linear second-order ODEs together with a linear readout,

$$ $\displaystyle{\bf y}^{\prime\prime}(t)$ $\displaystyle=-{\bf A}{\bf y}(t)+{\bf B}{\bf u}(t)+{\bf b},$ (1) $\displaystyle{\bf x}(t)$ $\displaystyle={\bf C}{\bf y}(t)+{\bf D}{\bf u}(t),$ $$

with hidden state ${\bf y}(t)\in\mathbb{R}^{m}$, output state ${\bf x}(t)\in{\mathbb{R}}^{q}$, time-dependent input signal ${\bf u}(t)\in\mathbb{R}^{p}$, weights ${\bf A}\in\mathbb{R}^{m\times m}$, ${\bf B}\in\mathbb{R}^{m\times p}$, ${\bf C}\in{\mathbb{R}}^{q\times m}$, ${\bf D}\in{\mathbb{R}}^{q\times p}$, and bias ${\bf b}\in{\mathbb{R}}^{m}$. Note that ${\bf A}$ is a diagonal matrix, i.e., with non-zero entries only on its diagonal. We further introduce an auxiliary state ${\bf z}(t)\in\mathbb{R}^{m}$, with ${\bf z}={\bf y}^{\prime}$. We can thus write (1) (omitting the bias ${\bf b}$ and linear readout ${\bf x}$) equivalently as,

$$ $\displaystyle{\bf z}^{\prime}(t)$ $\displaystyle=-{\bf A}{\bf y}(t)+{\bf B}{\bf u}(t),$ (2) $\displaystyle{\bf y}^{\prime}(t)$ $\displaystyle={\bf z}(t).$ $$

Figure: Figure 1: Schematic drawing of the proposed Linear Oscillatory State-Space model (LinOSS). The input sequences are processed through multiple LinOSS blocks. Each block is composed of a LinOSS layer followed by a nonlinear transformation, specifically a Gated Linear Units (Dauphin et al., ) (GLU) layer in our case. After passing through several LinOSS blocks, the latent sequences are decoded to produce the final output sequence.
Refer to caption: x1.png

### 2.1 Motivation and background

To demonstrate that the underlying ODE in (1) models a network of forced harmonic oscillators, we begin with the scalar case by setting $p=m=1$ in (1). Choosing ${\bf B}={\bf b}=0$, we obtain the classic ODE, ${\bf y}^{\prime\prime}+{\bf A}{\bf y}=0$, which describes simple harmonic motion with frequency ${\bf A}>0$, such as that of a simple pendulum (Guckenheimer & Holmes, ). Now, allowing ${\bf B}\neq 0$ introduces external forcing proportional to the input signal ${\bf u}(t)$, where ${\bf B}$ modulates the effect of the forcing. Finally, setting $p,m>1$ yields an uncoupled system of forced harmonic oscillators.

##### Neuroscience inspiration.

Our approach is inspired by neurobiology, where the periodic spiking and firing of action potentials in individual neurons can be observed and analyzed as oscillatory phenomena. Furthermore, entire networks of cortical neurons exhibit behavior-dependent oscillations, reviewed in Buzsaki & Draguhn (). Remarkably, despite their complexity, these neural oscillations share characteristics with harmonic oscillators, as described in Winfree (). Building on this insight, we distill the core essence of cortical dynamics and aim to construct a machine learning model following the motion of harmonic oscillations. Finally, we note that our focus is on the theoretical and empirical aspects of the proposed oscillatory state-space model in this manuscript. However, our model can be further used in the context of computational neuroscience emulating characteristic phenomena of brain oscillations, such as frequency-varying oscillations, transient synchronization or desynchronization of discharges, entrainment, phase shifts, and resonance, while at the same time being able to learn non-trivial input-output relations.

### 2.2 Comparison with related work

State-space models have seen continuous advancements since their introduction to modern sequence modeling in Gu et al. (). The original S4 model (Gu et al., ), along with its adaptations (Gu et al., ; Nguyen et al., ; Goel et al., ), utilized FFT to solve linear recurrences. More recently, a simplified variant, S5 (Smith et al., ), was introduced, which instead employs associative parallel scans to achieve similar computational speed. Another advancement to S4 have been Liquid Structural State-Space models (Hasani et al., ) which exchanged the static state matrix with an input-dependent state transition module. While all aforementioned models rely on the HiPPO parameterization (Gu et al., ), Linear Recurrent Units (LURs) (Orvieto et al., ) demonstrated that even simpler parameterization yield state-of-the-art performance.
Finally, selective state spaces have been introduced in Gu & Dao () in order to further close the performance gap between state-space models and Transformers utilized in foundation models.

Oscillatory dynamics as a neural computational paradigm has been originally introduced via Coupled Oscillatory RNNs (CoRNNs) (Rusch & Mishra, ). In this work, it has been shown that recurrent models based on nonlinear oscillatory dynamics are able to learn long-range interactions by mitigating the exploding and vanishing gradients problem (Pascanu, ). This approach was later refined for handling very long sequences in Rusch & Mishra (), which utilized uncoupled nonlinear oscillators for fast sequence modeling by adopting a diagonal hidden-to-hidden weight matrix. Interestingly, this method resembles modern state-space models but distinguishes itself by employing nonlinear dynamics. Since then, this concept, generally termed as neural oscillators (Lanthaler et al., ), has been extended to other areas of machine learning, e.g., Graph-Coupled Oscillator Networks (GraphCON) (Rusch et al., ) for learning on graph-structured data, and RNNs with oscillatory dynamics combined with convolutions leading to locally coupled oscillatory RNNs (Keller & Welling, ). Our proposed LinOSS model differs from all these methods by its explicit use of harmonic oscillators via a state-space model approach.

### 2.3 Discretization

Our aim is to approximately solve the linear ODE system (2) as fast as possible, while at the same time being able to guarantee stability over long time-scales. To this end, we suggest to leverage the following two time integration schemes.

##### Implicit time integration.

We fix a timestep $0<{\Delta t}\leq 1$ and define our proposed state-space model hidden states at time $t_{n}=n{\Delta t}$ as the following implicit discretization of the first order system (2):

$$ $\displaystyle{\bf z}_{n}$ $\displaystyle={\bf z}_{n-1}+{\Delta t}(-{\bf A}{\bf y}_{n}+{\bf B}{\bf u}_{n}),$ $\displaystyle{\bf y}_{n}$ $\displaystyle={\bf y}_{n-1}+{\Delta t}{\bf z}_{n}.$ $$

This can be written in matrix form by introducing ${\bf x}_{n}=[{\bf z}_{n},{\bf y}_{n}]^{\top}$,

$$ ${\bf M}{\bf x}_{n}={\bf x}_{n-1}+{\bf F}_{n},$ $$

with

$$ ${\bf M}=\begin{bmatrix}\mathrm{I}&{\Delta t}{\bf A}\\ -{\Delta t}\mathrm{I}&\mathrm{I}\end{bmatrix},\qquad{\bf F}_{n}=\begin{bmatrix }{\Delta t}{\bf B}{\bf u}_{n}\\ \bf 0\end{bmatrix}.$ $$

We can now solve our discretized oscillatory system by simply inverting matrix ${\bf M}$ and computing the induced recurrence. Note that inverting matrices typically requires $\mathcal{O}(m^{3})$ operations (where $m$ is the number of rows or columns of the matrix) using methods like Gauss-Jordan elimination, making it computationally expensive. However, by leveraging the Schur complement, we can obtain an explicit form of ${\bf M}^{-1}$ that can be computed in $\mathcal{O}(m)$, thanks to the diagonal structure of ${\bf A}$. More concretely,

$$ ${\bf M}^{-1}=\begin{bmatrix}\mathrm{I}-{\Delta t}^{2}{\bf A}{\bf S}&-{\Delta t }{\bf A}{\bf S}\\ {\Delta t}{\bf S}&{\bf S}\end{bmatrix}=\begin{bmatrix}{\bf S}&-{\Delta t}{\bf A }{\bf S}\\ {\Delta t}{\bf S}&{\bf S}\end{bmatrix},$ (3) $$

with the inverse of the Schur complement ${\bf S}=(\mathrm{I}+{\Delta t}^{2}{\bf A})^{-1}$ which itself is a diagonal matrix and can thus be trivially inverted.
We point out that among other choices, a straightforward condition that ensures ${\bf S}$ is well-defined is ${\bf A}_{k}\geq 0$ for all $k=1,\dots,m$.
The recurrence of the proposed model is then given as,

$$ ${\bf x}_{n}={\bf M}^{\text{IM}}{\bf x}_{n-1}+{\bf F}_{n}^{\text{IM}},$ (4) $$

with ${\bf M}^{\text{IM}}={\bf M}^{-1}$, and ${\bf F}^{\text{IM}}_{n}={\bf M}^{-1}{\bf F}_{n}$.
As we will show in the subsequent section, this discretization leads to a globally asymptotically stable discrete dynamical system.

##### Implicit-explicit time integration (IMEX).

Another discretization yielding stable dynamics that, however, do not converge exponentially fast to a steady-state can be obtained by leveraging symplectic integrators. To this end, we fix again a timestep $0<{\Delta t}\leq 1$ and define our proposed state-space model hidden states at time $t_{n}=n{\Delta t}$ as the following implicit-explicit (IMEX) discretization of the first order system (2):

$$ $\displaystyle{\bf z}_{n}$ $\displaystyle={\bf z}_{n-1}+{\Delta t}(-{\bf A}{\bf y}_{n-1}+{\bf B}{\bf u}_{n }),$ (5) $\displaystyle{\bf y}_{n}$ $\displaystyle={\bf y}_{n-1}+{\Delta t}{\bf z}_{n}.$ $$

The only difference compared to the previous fully implicit discretization is the explicit treatment of the hidden state ${\bf y}$ (i.e., using ${\bf y}_{n-1}$ instead of ${\bf y}_{n}$) in the first equation of (5). As before, we can simplify this system in matrix form,

$$ ${\bf x}_{n}={\bf M}^{\text{IMEX}}{\bf x}_{n-1}+{\bf F}_{n}^{\text{IMEX}},$ (6) $$

with

$$ ${\bf M}^{\text{IMEX}}=\begin{bmatrix}\mathrm{I}&-{\Delta t}{\bf A}\\ {\Delta t}\mathrm{I}&\mathrm{I}-{\Delta t}^{2}{\bf A}\end{bmatrix},\qquad{\bf F }_{n}^{\text{IMEX}}=\begin{bmatrix}{\Delta t}{\bf B}{\bf u}_{n}\\ {\Delta t}^{2}{\bf B}{\bf u}_{n}\end{bmatrix}.$ $$

Interestingly, ODE system (2) represents a Hamiltonian system (Arnold, ), with Hamiltonian,

$$ ${\bf H}({\bf y},{\bf z},t)=\frac{1}{2}\sum_{k=1}^{m}{\bf A}_{k}{\bf y}_{k}^{2} +{\bf z}_{k}^{2}-2\left(\sum_{l=1}^{p}{\bf B}_{kl}{\bf u}(t)_{l}\right){\bf y} _{k}.$ (7) $$

The numerical approximation of this system using the previously described IMEX discretization is symplectic, i.e., it preserves a Hamiltonian close to the Hamiltonian of the continuous system. Thus, by the well-known Liouville’s theorem (Sanz Serna & Calvo, ), we know that the phase space volume of (2) as well as of its symplectic approximation (6) is preserved. This gives rise to invertible model architectures leading to memory efficient implementations of the backpropagation through time algorithm, similar as in (Rusch & Mishra, ). This denotes the most significant difference between the two different discretization schemes, i.e., the IMEX integration-based model is volume preserving, while IM integration-based model introduces dissipative terms. We will see in subsequent sections that both models have their own individual advantages depending on the underlying data.

### 2.4 Fast recurrence via associative parallel scans

Parallel (or associative) scans, first introduced in Kogge & Stone () and reviewed in Blelloch (), offer a powerful method for drastically reducing the computational time of recurrent operations. These scans have previously been employed to enhance the training and inference speed of RNNs (Martin & Cundy, ; Kaul, ). This technique was later adapted for state-space models in Smith et al. (), becoming a crucial component in the development of many state-of-the-art sequence models, including Linear Recurrent Units (LRUs) (Orvieto et al., ) and Mamba models (Gu & Dao, ).

A parallel scan operates on a sequence $[x_{1},\dots,x_{N}]$ with a binary associative operation $\bullet$, i.e., an operation satisfying $(x\bullet y)\bullet z=x\bullet(y\bullet z)$ for instances $x,y$, and $z$, to return the sequence $[x_{1},x_{1}\bullet x_{2},\dots,x_{1}\bullet x_{2}\bullet\dots\bullet x_{N}]$. Under certain assumptions, this operation can be performed in computational time proportional to $\lceil\log_{2}(N)\rceil$. This is in stark contrast to the computational time of serial recurrence that is proportional to $N$. It is straightforward to check that the following operation is associative:

$$ $({\bf a}_{1},{\bf a}_{2})\bullet({\bf b}_{1},{\bf b}_{2})=({\bf b}_{1}\mathbin {\vphantom{\circ}\text{\ooalign{\leavevmode\leavevmode\hbox{ \leavevmode\hbox to 5pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox to 0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{ 0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill {0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ } \nullfont\hbox to0.0pt{{}{}{}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{0.0 pt}{0.0pt}\pgfsys@lineto{5.00002pt}{0.0pt}\pgfsys@lineto{5.00002pt}{0.0pt} \pgfsys@closepath\pgfsys@clipnext\pgfsys@discardpath\pgfsys@invoke{ }{{{}{}{{} }{} {{}{{}}}{{}{}}{}{{}{}} {{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1. 0}{0.0pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\leavevmode\hbox{\hskip 0.0pt\hbox{ \set@color{$\bullet$}}}}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} } {}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope } \pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}\cr$\circ$\cr}}}{\bf a }_{1},{\bf b}_{1}\mathbin{\text{\ooalign{\raisebox{0.0pt}{\leavevmode \leavevmode\hbox{ \leavevmode\hbox to5pt{\vbox to4.44pt{\pgfpicture \makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox to0.0pt{\pgfsys@beginscope \pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0} \pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0} {0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to 0.0pt{{}{}{}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{0.0pt}{4.44444pt} \pgfsys@lineto{5.00002pt}{4.44444pt}\pgfsys@lineto{5.00002pt}{0.0pt} \pgfsys@closepath\pgfsys@clipnext\pgfsys@discardpath\pgfsys@invoke{ }{{{}{}{{} }{} {{}{{}}}{{}{}}{}{{}{}} {{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1. 0}{0.0pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\leavevmode\hbox{\hskip 0.0pt\hbox{ \set@color{$\bullet$}}}}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} } {}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope } \pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}\cr$\circ$\cr}}}{\bf a }_{2}+{\bf b}_{2}),$ $$

where $\mathbin{\vphantom{\circ}\text{\ooalign{\leavevmode\leavevmode\hbox{
\leavevmode\hbox to5pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt
\lower 0.0pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{
pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }
\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}
\pgfsys@invoke{ }\nullfont\hbox to0.0pt{{}{}{}{}\pgfsys@moveto{0.0pt}{0.0pt}
\pgfsys@lineto{0.0pt}{0.0pt}\pgfsys@lineto{5.00002pt}{0.0pt}\pgfsys@lineto{5.0
0002pt}{0.0pt}\pgfsys@closepath\pgfsys@clipnext\pgfsys@discardpath
\pgfsys@invoke{ }{{{}{}{{}}{}
{{}{{}}}{{}{}}{}{{}{}}
{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.
0}{0.0pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\leavevmode\hbox{\hskip 0.0pt\hbox{
\set@color{$\bullet$}}}}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
}
{}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }
\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}\cr$\circ$\cr}}},
\mathbin{\text{\ooalign{\raisebox{0.0pt}{\leavevmode\leavevmode\hbox{
\leavevmode\hbox to5pt{\vbox to4.44pt{\pgfpicture\makeatletter\hbox{\hskip 0.0
pt\lower 0.0pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{
pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }
\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}
\pgfsys@invoke{ }\nullfont\hbox to0.0pt{{}{}{}{}\pgfsys@moveto{0.0pt}{0.0pt}
\pgfsys@lineto{0.0pt}{4.44444pt}\pgfsys@lineto{5.00002pt}{4.44444pt}
\pgfsys@lineto{5.00002pt}{0.0pt}\pgfsys@closepath\pgfsys@clipnext
\pgfsys@discardpath\pgfsys@invoke{ }{{{}{}{{}}{}
{{}{{}}}{{}{}}{}{{}{}}
{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.
0}{0.0pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\leavevmode\hbox{\hskip 0.0pt\hbox{
\set@color{$\bullet$}}}}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
}
{}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }
\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}\cr$\circ$\cr}}}$ are the matrix-matrix and matrix-vector products. Note that both products can be computed in $\mathcal{O}(m)$ time leveraging $2\times 2$ block matrices with only diagonal entries in each block, e.g., matrices ${\bf M}^{\text{IM}}$ and ${\bf M}^{\text{IMEX}}$ in (4),(6). Clearly, applying a parallel scan based on this associative operation on the input sequence $[({\bf M},{\bf F}_{1}),({\bf M},{\bf F}_{2}),\dots]$ yields an output sequence $[({\bf M},{\bf F}_{1}),({\bf M}^{2},{\bf M}{\bf F}_{1}+{\bf F}_{2}),\dots]$ where the solution of the recurrent system ${\bf x}_{n}={\bf M}{\bf x}_{n-1}+{\bf F}_{n}$ (with initial value ${\bf x}_{0}=\bf 0$) is stored in the second argument of the elements of this sequence. Thus, we can successfully apply parallel scans to both discretizations (4)(6) of our
proposed system in order to significantly speed up computations. We refer to the application of this parallel scan to implicit formulations as implicit parallel scans. When applied to implicit-explicit formulations, we describe it as implicit-explicit parallel scans.

We term our proposed sequence model, which efficiently solves the underlying linear system of harmonic oscillators (2) using fast parallel scans, as Linear Oscillatory State-Space (LinOSS) model. To differentiate between the two discretization methods, we refer to the model derived from implicit discretization as LinOSS-IM, and the model based on the implicit-explicit discretization as LinOSS-IMEX. A schematic drawing of our proposed state-space model can be seen in Fig. 1.

## 3 Theoretical insights

### 3.1 Stability and learning long-range interactions

Hidden states of classical nonlinear RNNs are updated based on its previous hidden states pushed through a parametric function followed by a (usually) bounded nonlinear activation function such as tanh or sigmoid. This way, the hidden states are guaranteed to not blow up. However, this is not true for linear recurrences, where specific structures of the hidden-to-hidden weight matrix can lead to unstable dynamics. The following simple argument demonstrates that computing the eigenspectrum (i.e., set of eigenvalues) of the hidden-to-hidden weight matrix of a linear recurrent system suffices in order to analyse its stability properties. To this end, let us consider a general linear discrete dynamical system with external forcing given via the recurrence ${\bf x}_{n}={\bf M}{\bf x}_{n-1}+{\bf B}{\bf u}_{n}$, and initial value ${\bf x}_{0}=\bf 0$. Assuming ${\bf M}$ is diagonalizable (if not, one can make a similar argument leveraging the Jordan normal form), i.e., there exists matrix ${\bf S}$ such that ${\bf M}={\bf S}{\bf\Lambda}{\bf S}^{-1}$, where $\bf\Lambda$ is a diagonal matrix with the eigenvalues of ${\bf M}$ on its diagonal. The dynamics of the transformed hidden state $\bar{{\bf x}}_{n}={\bf S}^{-1}{\bf x}_{n}$ evolve according to $\bar{{\bf x}}_{n}={\bf S}^{-1}{\bf x}_{n}={\bf S}^{-1}{\bf M}{\bf x}_{n-1}+{
\bf S}^{-1}{\bf B}{\bf u}_{n}={\bf\Lambda}\bar{{\bf x}}_{n-1}+\bar{{\bf B}}{
\bf u}_{n}$ with initial value $\bar{{\bf x}}_{0}={\bf S}^{-1}{\bf x}_{0}$, where $\bar{{\bf B}}={\bf S}^{-1}{\bf B}$. Unrolling the dynamics (assuming ${\bf x}_{0}=\bf 0$) yields,

$$ $\bar{{\bf x}}_{1}=\bar{{\bf B}}{\bf u}_{1},\quad\bar{{\bf x}}_{2}={\bf\Lambda} \bar{{\bf B}}{\bf u}_{1}+\bar{{\bf B}}{\bf u}_{2},\quad...\quad\Rightarrow{\bf x }_{n}=\sum_{k=0}^{n-1}{\bf\Lambda}^{k}\bar{{\bf B}}{\bf u}_{n-k}.$ $$

Clearly, if all eigenvalues $\bf\Lambda$ have magnitude less or equal than 1 the hidden states will not blow up. In addition to stability guarantees, eigenspectra with unit norm allow the model to learn long-range interactions (Arjovsky et al., ; Gu et al., ; Orvieto et al., ) by avoiding vanishing and exploding gradients (Pascanu, ). Therefore, it is sufficient to analyse the eigenspectra of our proposed LinOSS models in order to understand their ability to generate stable dynamics and learn long-range interactions. To this end, we have the following propositions.

###### Proposition 3.1 .

Let ${\bf M}^{\text{IM}}\in{\mathbb{R}}^{m\times m}$ be the hidden-to-hidden weight matrix of the implicit model LinOSS-IM (4). We assume that ${\bf A}_{kk}\geq 0$ for all diagonal elements $k=1,\dots,m$ of ${\bf A}$, and further that ${\Delta t}>0$. Then, the complex eigenvalues of ${\bf M}^{\text{IM}}$ are given as,

$$ $\lambda_{j}=\frac{1}{1+{\Delta t}^{2}{\bf A}_{kk}}+i(-1)^{\lceil\frac{j}{m} \rceil}{\Delta t}\frac{\sqrt{{\bf A}_{kk}}}{1+{\Delta t}^{2}{\bf A}_{kk}}, \quad\text{for all }j=1,\dots,2m,$ $$

with $k=j\text{ mod }m$.
Moreover, the spectral radius $\rho({\bf M}^{\text{IM}})$ is bounded by $1$, i.e., $|\lambda_{j}|\leq 1$ for all $j=1,\dots,2m$.

###### Proof.

$$ $\displaystyle\det({\bf M}^{\text{IM}}-\lambda\mathrm{I})$ $\displaystyle=\begin{vmatrix}{\bf S}-\lambda\mathrm{I}&-{\Delta t}{\bf A}{\bf S }\\ {\Delta t}{\bf S}&{\bf S}-\lambda\mathrm{I}\end{vmatrix}=\begin{vmatrix}{\bf S }-\lambda\mathrm{I}&-{\Delta t}{\bf A}{\bf S}\\ {\bf 0}&{\bf S}-\lambda\mathrm{I}+{\Delta t}^{2}{\bf A}{\bf S}^{2}({\bf S}- \lambda\mathrm{I})^{-1}\end{vmatrix}$ $\displaystyle=\prod_{k=1}^{m}({\bf S}_{kk}-\lambda)\left({\bf S}_{kk}-\lambda+ \frac{{\Delta t}^{2}{\bf A}_{kk}{\bf S}_{kk}^{2}}{{\bf S}_{kk}-\lambda}\right) =\prod_{k=1}^{m}[({\bf S}_{kk}-\lambda)^{2}+{\Delta t}^{2}{\bf A}_{kk}{\bf S}_ {kk}^{2}].$ $$

Setting $k=j\text{ mod }m$, the eigenvalues of ${\bf M}^{\text{IM}}$ are thus given as,

$$ $\lambda_{j}={\bf S}_{kk}+i(-1)^{\lceil\frac{j}{m}\rceil}{\Delta t}{\bf S}_{kk} \sqrt{{\bf A}_{kk}},\quad\text{for all }j=1,\dots,2m.$ $$

In particular, assuming ${\bf A}_{kk}\geq 0$ for all $k=1,\dots,m$, the magnitude of the eigenvalues $\lambda_{j}$ are given as,

$$ $|\lambda_{j}|^{2}={\bf S}_{kk}^{2}+{\Delta t}^{2}{\bf S}_{kk}^{2}{\bf A}_{kk}= {\bf S}_{kk}^{2}(1+{\Delta t}^{2}{\bf A}_{kk})={\bf S}_{kk}\leq 1,$ $$

for all $j=1,\dots,2m$, with $k=j\text{ mod }m$.
∎

This proof reveals important insights into our proposed LinOSS-IM model. First, the magnitude of the eigenvalues at initialization can be controlled by either specific initialization of ${\bf A}$, or alternatively through a specific choice of the time-step parameter ${\Delta t}$. Moreover, LinOSS-IM yields asymptotically stable dynamics for any choices of positive parameters ${\bf A}$. This denotes a significant difference compared to previous first-order system, where the values of ${\bf A}$ (and thus its eigenvalues) have to be heavily constraint for the system to be stable. We argue that this flexibility in the parameterization of our model benefits the optimizer potentially leading to better performance in practice.

###### Remark 1 .

A straightforward adaptation of Proposition 3.1 can also be derived for the implicit-explicit version of the model, LinOSS-IMEX (6). In particular, a similar analysis shows that the eigenvalues $\lambda_{j}$ of ${\bf M}^{\text{IMEX}}$ in (6) all satisfy $|\lambda_{j}|=1$. This result underscores the key distinction between the two models: LinOSS-IM incorporates dissipative terms, introducing a forgetting mechanism into the state-space model, whereas LinOSS-IMEX allows for learning interactions over arbitrarily long sequences. As previously argued in the literature, forgetting is a crucial feature for many applications. Therefore, LinOSS-IM can be seen as a model that integrates forgetting, while LinOSS-IMEX is capable of learning sequential interactions over arbitrary time-scales.

###### Proposition 3.1 .

###### Proof.

###### Remark 1 .

##### Initialization and parameterization of weights.

Proposition 3.1 and Remark 1 reveal that both LinOSS models exhibit stable dynamics as long as the diagonal weights ${\bf A}$ in (1) are non-negative. This condition can easily be fulfilled via many different parameterizations of the diagonal matrix ${\bf A}$. An obvious choice is to square the diagonal values, i.e., a parameterization ${\bf A}=\hat{{\bf A}}\hat{{\bf A}}$, with diagonal matrix $\hat{{\bf A}}\in{\mathbb{R}}^{m\times m}$. Another straightforward approach is to apply the element-wise ReLU nonlinear activation function to ${\bf A}$, i.e., ${\bf A}=\operatorname{ReLU}(\hat{{\bf A}})$, with $\operatorname{ReLU}(x)=\max(0,x)$ and $\hat{{\bf A}}$ as before. The latter results in a LinOSS model where specific dimensions can be switched off completely by setting the corresponding weight in $\hat{{\bf A}}$ to a negative value. Due to this flexibility, we decide to focus on the ReLU-parameterization of the diagonal weights ${\bf A}$ in this manuscript.

After discussing the parameterization of ${\bf A}$, the next step is to determine an appropriate method for initializing its weights prior to training. Since we already know from Remark 1 that all absolute eigenvalues of the hidden matrix ${\bf M}^{\text{IMEX}}$ of LinOSS-IMEX are exactly one, the specific values of ${\bf A}$ are irrelevant for the model to learn long-range interactions. However, according to Proposition 3.1, ${\bf A}$ highly influences the eigenvalues of the hidden matrix ${\bf M}^{\text{IM}}$ of LinOSS-IM (4). As discussed in the previous section, high powers of the absolute eigenvalues of the hidden matrix are of particular interest when learning long-range interactions. To this end, we have the following Proposition concerning the expected powers of absolute eigenvalues for LinOSS-IM.

###### Proposition 3.2 .

Let $\{\lambda_{j}\}_{j=1}^{2m}$ be the eigenspectrum of the hidden-to-hidden matrix ${\bf M}^{\text{IM}}$ of the LinOSS-IM model (4). We further initialize ${\bf A}_{kk}\sim\mathcal{U}([0,A_{\text{max}}])$ with $A_{\text{max}}>0$ for all diagonal elements $k=1,\dots,m$ of ${\bf A}$ in (2). Then, the $N$-th moment of the magnitude of the eigenvalues are given as,

$$ $\mathbb{E}(|\lambda_{j}|^{N})=\frac{({\Delta t}^{2}A_{\text{max}}+1)^{1-\frac{ N}{2}}-1}{{\Delta t}^{2}A_{\text{max}}(1-\frac{N}{2})},$ (8) $$

for all $j=1,\dots,2m$, with $k=j\text{ mod }m$.

The proof, detailed in Appendix Section B.1, is a straight-forward application of the law of the unconscious statistician. Proposition 3.2 demonstrates that while the expectation of $|\lambda_{i}|^{N}$ might be much smaller than $1$, it is still sufficiently large for practical use-cases. For instance, even considering ${\Delta t}=A_{\text{max}}=1$ and an extreme sequence length of $N=100$k still yields $\mathbb{E}(|\lambda_{j}|^{N})=1/49999\approx 2\times 10^{-5}$. Based on this and the fact that the values of ${\bf A}$ do not affect the eigenvalues for LinOSS-IMEX, we decide to initialize ${\bf A}$ according to ${\bf A}_{kk}\sim\mathcal{U}([0,1])$ for both LinOSS models, while setting ${\Delta t}=1$.

###### Proposition 3.2 .

### 3.2 Universality of LinOSS within continuous and causal operators

While trainability is an important aspect of learning long-range interactions, it does not demonstrate why LinOSS is able to express complex mappings between general (i.e., not necessarily oscillatory) input and output sequences. Therefore, in this section we analyze the approximation power of our proposed LinOSS model. Following the recent work of Lanthaler et al. () we show that LinOSS is universal within the class of continuous and causal operators between time-series.
To this end, we consider a full LinOSS block,

$$ ${\bf z}(t)={\bf W}\sigma(\tilde{{\bf W}}{\bf y}(t)+\tilde{{\bf b}}),$ (9) $$

with weights ${\bf W}\in{\mathbb{R}}^{q\times\tilde{m}}$, $\tilde{\bf W}\in{\mathbb{R}}^{\tilde{m}\times m}$, bias $\tilde{{\bf b}}\in{\mathbb{R}}^{\tilde{m}}$, element-wise nonlinear activation function $\sigma$ (e.g., $\tanh$ or ReLU), and solution ${\bf y}(t)$ of the LinOSS differential equations (2).

Based on this, we are now interested in approximating operators $\Phi:C_{0}([0,T];{\mathbb{R}}^{p})\rightarrow C_{0}([0,T];{\mathbb{R}}^{q})$ with the full LinOSS model (9), where

$$ $C_{0}([0,T];{\mathbb{R}}^{p}):={\left\{{\bf u}:[0,T]\to{\mathbb{R}}^{p}\, \middle|\,t\mapsto{\bf u}(t)\text{ is continuous and }{\bf u}(0)=\bf 0\right\}}.$ $$

As pointed out in Lanthaler et al. (), the condition ${\bf u}(0)=\bf 0$ is not restrictive and can directly be generalized to the case of any initial condition ${\bf u}(0)={\bf u}_{0}\in{\mathbb{R}}^{p}$ simply by introducing an arbitrarily small warm-up phase $[-t_{0},0]$ with $t_{0}>0$ of the oscillators to synchronize with the input signal ${\bf u}$. In addition to these function spaces, we pose the following conditions on the underlying operator $\Phi$,

- 1.
$\Phi$ is *causal*, i.e., for any $t\in[0,T]$, if ${\bf u},{\bf B}\in C_{0}([0,T];{\mathbb{R}}^{p})$ are two input signals, such that ${\bf u}|_{[0,t]}\equiv{\bf B}|_{[0,t]}$, then $\Phi({\bf u})(t)=\Phi({\bf B})(t)$.
- 2.
$\Phi$ is *continuous* as an operator
$\Phi:(C_{0}([0,T];{\mathbb{R}}^{p}),\|{\,\cdot\,}\|_{L^{\infty}})\to(C_{0}([0,
T];{\mathbb{R}}^{q}),\|{\,\cdot\,}\|_{L^{\infty}}),$
with respect to the $L^{\infty}$-norm on the input-/output-signals.

###### Theorem 3.3 .

Let $\Phi:C_{0}([0,T];{\mathbb{R}}^{p})\to C_{0}([0,T];{\mathbb{R}}^{q})$ be a causal and continuous operator. Let $K\subset C_{0}([0,T];{\mathbb{R}}^{p})$ be compact. Then for any $\epsilon>0$, there exist hyperparameters $m$, $\tilde{m}$, diagonal weight matrix ${\bf A}\in{\mathbb{R}}^{m\times m}$, weights ${\bf B}\in{\mathbb{R}}^{m\times d}$, $\tilde{{\bf W}}\in{\mathbb{R}}^{\tilde{m}\times m}$, ${\bf W}\in{\mathbb{R}}^{q\times\tilde{m}}$ and bias vectors ${\bf b}\in{\mathbb{R}}^{m}$, $\tilde{{\bf b}}\in{\mathbb{R}}^{\tilde{m}}$, such that the output ${\bf z}:[0,T]\to{\mathbb{R}}^{q}$ of the LinOSS block (9) satisfies,

$$ $\sup_{t\in[0,T]}|\Phi({\bf u})(t)-{\bf z}(t)|\leq\epsilon,\quad\forall\,{\bf u }\in K.$ $$

The proof can be found in Appendix Section B.2. The main idea of the proof is to encode the infinite-dimensional operator $\Phi$ with a finite-dimensional operator that makes use of the structure of the LinOSS ODE system (1), and that can further be expressed by a (finite-dimensional) function.

###### Theorem 3.3 .

## 4 Empirical Results

In this section, we empirically test the performance of our proposed LinOSS models on a variety of challenging real-world sequential datasets, ranging from scientific datasets in genomics to practical applications in medicine. We ensure thereby a fair comparison to other state-of-the-art sequence models such as Mamba, LRU, and S5.

### 4.1 Learning Long-Range Interactions

In the first part of the experiments, we focus on a recently proposed long-range sequential benchmark introduced in Walker et al. (). This benchmark focuses on six datasets from the UEA Multivariate Time Series Classification Archive (UEA-MTSCA), selecting those with the longest sequences for increased difficulty. The sequence lengths range thereby from $400$ to almost $18$k. We compare our proposed LinOSS models to recent state-of-the-art sequence models, including state-space models such as Mamba and LRU. Table 1 shows the test accuracies averaged over 5 random model initialization and dataset splits of all six datasets for LinOSS as well as competing methods. We note that all other results are taken from Walker et al. (). Moreover, we highlight that we exactly follow the training procedure described in Walker et al. () in order to ensure a fair comparison against competing models. More concretely, we use the same pre-selected random seeds for splitting the datasets into training, validation, and testing parts (using $70/15/15$ splits), as well as tune our model hyperparameters only on the same pre-described grid. In fact, since LinOSS does not possess any model specific hyperparameters – unlike competing state-space models such as LRU or S5 – our search grid is lower-dimensional compared to the other models considered. As a result, the hyperparameter tuning process involves significantly fewer model instances.

**Table 1: Test accuracies averaged over $5$ training runs on UAE time-series classification datasets. All models are trained based on the same hyper-parameter tuning protocol in order to ensure fair comparability. The dataset names are abbreviations of the following UEA time-series datasets: EigenWorms (Worms), SelfRegulationSCP1 (SCP1), SelfRegulationSCP2 (SCP2), EthanolConcentration (Ethanol), Heartbeat, MotorImagery (Motor). The three best performing methods are highlighted in red (First), blue (Second), and violet (Third).**
|  | Worms | SCP1 | SCP2 | Ethanol | Heartbeat | Motor | Avg |
| --- | --- | --- | --- | --- | --- | --- | --- |
| Seq. length | 17,984 | 896 | 1,152 | 1,751 | 405 | 3,000 |  |
| #Classes | 5 | 2 | 2 | 4 | 2 | 2 |  |
| NRDE | $77.2\pm 7.1$ | $76.7\pm 5.6$ | $48.1\pm 11.4$ | $\bf\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{31.4\pm 4<br>.5}$ | $73.9\pm 2.6$ | $54.0\pm 7.8$ | $60.2$ |
| NCDE | $62.2\pm 2.2$ | $80.0\pm 2.0$ | $49.1\pm 6.2$ | $22.0\pm 1.0$ | $68.1\pm 5.8$ | $51.6\pm 6.2$ | $55.5$ |
| Log-NCDE | $82.8\pm 2.7$ | $82.1\pm 1.4$ | $54.0\pm 2.6$ | $\bf\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}{35.9\pm 6<br>.1}$ | $74.2\pm 2.0$ | $\bf\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}{57.2<br>\pm 5.6}$ | $\bf\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}{64.4}$ |
| LRU | $\bf\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{85.0\pm 6<br>.2}$ | $84.5\pm 4.6$ | $47.4\pm 4.0$ | $23.8\pm 2.8$ | $\bf\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}{78.1\pm 7<br>.6}$ | $51.9\pm 8.6$ | $61.8$ |
| S5 | $\bf\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}{83.9<br>\pm 4.1}$ | $\bf\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}{87.1<br>\pm 2.1}$ | $\bf\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}{55.1<br>\pm 3.3}$ | $25.6\pm 3.5$ | $73.9\pm 3.1$ | $53.0\pm 3.9$ | $63.1$ |
| Mamba | $70.9\pm 15.8$ | $80.7\pm 1.4$ | $48.2\pm 3.9$ | $27.9\pm 4.5$ | $\bf\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}{76.2<br>\pm 3.8}$ | $47.7\pm 4.5$ | $58.6$ |
| S6 | $\bf\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{85.0\pm 1<br>6.1}$ | $82.8\pm 2.7$ | $49.9\pm 9.4$ | $26.4\pm 6.4$ | $\bf\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{76.5\pm 8<br>.3}$ | $51.3\pm 4.7$ | $62.0$ |
| LinOSS-IMEX | $80.0\pm 2.7$ | $\bf\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{87.5\pm 4<br>.0}$ | $\bf\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}{58.9\pm 8<br>.1}$ | $\bf\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}{29.9<br>\pm 1.0}$ | $75.5\pm 4.3$ | $\bf\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{57.9\pm 5<br>.3}$ | $\bf\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{65.0}$ |
| LinOSS-IM | $\bf\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}{95.0\pm 4<br>.4}$ | $\bf\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}{87.8\pm 2<br>.6}$ | $\bf\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{58.2\pm 6<br>.9}$ | $\bf\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}{29.9<br>\pm 0.6}$ | $75.8\pm 3.7$ | $\bf\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}{60.0\pm 7<br>.5}$ | $\bf\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}{67.8}$ |

We can see in Table 1 that on average both LinOSS models outperform any other model we consider here by reaching an average (over all six datasets) accuracy of $65.0\%$ for LinOSS-IMEX and $67.8\%$ for LinOSS-IM. In particular, the average accuracy of LinOSS-IM is significantly higher than the two next best models, i.e., Log-NCDE reaching an average accuracy of $64.4\%$, and S5 reaching an accuracy of $63.1\%$ on average. It is particularly noteworthy that LinOSS-IM yields state-of-the-art results on the two datasets with the longest sequences, namely EigenWorms and MotorImagery. Specifically, on the challenging EigenWorms dataset, LinOSS-IM improves the test accuracy from the previous state-of-the-art of $85\%$ to an impressive $95\%$.

### 4.2 Taking long ranges to the extreme

In this experiment, we test the performance of LinOSS in the case of taking long range interactions to the extreme. To this end, we consider the PPG-DaLiA dataset, a multivariate time series regression dataset designed for heart rate prediction using data collected from a wrist-worn device (Reiss et al., ). It includes recordings from fifteen individuals, each with approximately 150 minutes of data sampled at a maximum rate of 128 Hz. The dataset consists of six channels: blood volume pulse, electrodermal activity, body temperature, and three-axis acceleration. We follow Walker et al. () and divide the data into training, validation, and test sets with a 70:15:15 split for each individual. After splitting the data, a sliding window of length 49920 and step size 4992 is applied. As in previous experiments, we apply the exact same hyperparameter tuning protocol to each model we consider here to ensure fair comparison. The test mean-squared error (MSE) of both LinOSS models as well as other competing models are shown in Table 2. We can see that both LinOSS models significantly outperform all other models. In particular, LinOSS-IM outperforms Mamba by nearly a factor of $2$, as well as LRU by a factor of around $2.5$. This highlights the effectiveness of our proposed LinOSS models on sequential data with extreme length.

**Table 2: Average test mean-squared error over $5$ training runs on the PPG-DaLiA dataset. All models are trained following the same hyper-parameter tuning protocol in order to ensure fair comparability. The three best performing methods are highlighted in red (First), blue (Second), and violet (Third).**
| Model | MSE $\times 10^{-2}$ |
| --- | --- |
| NRDE (Morrill et al., ) | $18.79\pm 2.69$ |
| Log-NCDE (Walker et al., ) | $\bf\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}{9.20<br>\pm 0.55}$ |
| LRU (Orvieto et al., ) | $15.64\pm 1.67$ |
| S5 (Smith et al., ) | $12.79\pm 0.72$ |
| S6 (Gu & Dao, ) | $12.88\pm 2.05$ |
| Mamba (Gu & Dao, ) | $10.65\pm 2.20$ |
| LinOSS-IMEX | $\bf\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{7.5\pm 0.<br>46}$ |
| LinOSS-IM | $\bf\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}{6.4\pm 0.<br>23}$ |

### 4.3 Long-horizon forecasting

Inspired by Gu et al. (), we test our proposed LinOSS on its ability to serve as a general sequence-to-sequence model, even with weak inductive bias. To this end, we focus on time-series forecasting which typically requires specialized domain-specific preprocessing and architectures. We do, however, not alter our LinOSS model nor incorporate any inductive biases. Thus, we simply follow Gu et al. () by setting up LinOSS as a general sequence-to-sequence model that treats forecasting as a masked sequence-to-sequence transformation. We consider a weather prediction task introduced in Zhou et al. (). In this task, several climate variables are predicted into the future based on local climatological data. Here, we focus on the hardest task in Zhou et al. () of predicting the future 720 timesteps (hours) based on the past 720 timesteps. Table 3 shows the mean absolute error for both LinOSS models as well as other competing models. We can see that both LinOSS models outperform Transformers-based baselines as well as the other state-space models.

**Table 3: Mean absolute error on the weather dataset predicting the future 720 time steps based on the past 720 time steps. The three best performing methods are highlighted in red (First), blue (Second), and violet (Third).**
| Model | Mean Absolute Error |
| --- | --- |
| Informer (Zhou et al., ) | $0.731$ |
| Informer^† (Zhou et al., ) | $0.741$ |
| LogTrans (Li et al., ) | $0.773$ |
| Reformer (Kitaev et al., ) | $1.575$ |
| LSTMa (Bahdanau et al., ) | $1.109$ |
| LSTnet (Lai et al., ) | $0.757$ |
| S4 (Gu et al., ) | $\bf\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}{0.578}$ |
| LinOSS-IMEX | $\bf\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}{0.508}$ |
| LinOSS-IM | $\bf\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{0.528}$ |

## 5 Conclusion

In this paper, we introduce LinOSS, a state-space model based on harmonic oscillators. We rigorously show that LinOSS produces stable dynamics and is able to learn long-range interactions only requiring a nonnegative diagonal state matrix. In addition, we connect the underlying ODE system of the LinOSS model to Hamiltonian dynamics, which, discretized using symplectic integrators, perfectly preserves its symmetry of time reversibility. Moreover, we show that LinOSS is a universal approximator of continuous and causal operators between time-series. Together, these properties enable efficient modeling of long-range interactions, while ensuring stable and accurate long-horizon forecasting. Finally, we demonstrate that LinOSS outperforms state-of-the-art state-space models, such as Mamba, LRU, and S5.

#### Acknowledgments

The authors would like to thank Dr. Alaa Maalouf (MIT) and Dr. T. Anderson Keller (Harvard University) for their insightful feedback and constructive suggestions on an earlier version of the manuscript. This work was supported in part by the Postdoc.Mobility grant P500PT-217915 from the Swiss National Science Foundation, and the Department of the Air Force Artificial Intelligence Accelerator and was accomplished under Cooperative Agreement Number FA8750-19-2-1000. The views and conclusions contained in this document are those of the authors and should not be interpreted as representing the official policies, either expressed or implied, of the Department of the Air Force or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for Government purposes notwithstanding any copyright notation herein.

## Appendix A Training details

All experiments were conducted on Nvidia Tesla V100 GPUs, with the exception of the PPG experiment, which was run on Nvidia Tesla A100 GPUs due to higher memory demands. The hyperparameters of the models were optimized with the same grid search approach from Walker et al. () to ensure perfect comparability with competing methods, i.e., using the grid: learning rate $=\{0.00001,0.0001,0.001\}$, number of layers $=\{2,4,6\}$, number of hidden neurons $=\{16,64,128\}$, state-space dimension $=\{16,64,256\}$, include time dimension $\{\text{True},\text{False}\}$. The hyperparameters yielding the best performance can be seen in Table 4 for both LinOSS-IM and LinOSS-IMEX for each experiment in the main paper.

**Table 4: Best hyperparameters**
|  | Model | lr | hidden dim | state dim | #layers | include time |
| --- | --- | --- | --- | --- | --- | --- |
| Worms | IM | $0.001$ | $128$ | $64$ | $2$ | True |
| IMEX | $0.0001$ | $64$ | $16$ | $2$ | False |  |
| SCP1 | IM | $0.0001$ | $128$ | $256$ | $6$ | True |
| IMEX | $0.0001$ | $64$ | $256$ | $6$ | False |  |
| SCP2 | IM | $0.0001$ | $128$ | $64$ | $6$ | False |
| IMEX | $0.00001$ | $64$ | $256$ | $6$ | True |  |
| Ethanol | IM | $0.00001$ | $16$ | $16$ | $4$ | False |
| IMEX | $0.00001$ | $16$ | $256$ | $4$ | False |  |
| Heartbeat | IM | $0.001$ | $16$ | $16$ | $6$ | True |
| IMEX | $0.00001$ | $64$ | $16$ | $2$ | True |  |
| Motor | IM | $0.00001$ | $128$ | $16$ | $2$ | False |
| IMEX | $0.0001$ | $16$ | $256$ | $6$ | True |  |
| PPG | IM | $0.001$ | $16$ | $64$ | $6$ | True |
| IMEX | $0.0001$ | $64$ | $16$ | $2$ | True |  |

### A.1 Additional experiments

In the main paper, we have focused on parameterizing the state matrix ${\bf A}$ in (1) according to ${\bf A}=\operatorname{ReLU}(\hat{{\bf A}})$, with diagonal matrix $\hat{{\bf A}}\in{\mathbb{R}}^{m\times m}$. This was due to the fact that LinoSS requires the state matrix to be nonnegative in order to produce stable dynamics. However, another viable parameterization choice would be ${\bf A}=\hat{{\bf A}}\hat{{\bf A}}$. In this section, we test how the squared parameterization influences the performance on LinOSS on 6 long-range datasets taken from the main paper. Table 5 shows the average test accuracies of LinOSS-IM using a relu parameterization as well as using a squared parameterization together with the same baselines taken from the main paper. We can see that on average, the relu parameterizaton befores better. However, since the squared parameterization performs better on SCP2, Ethanol and Heartbeat (i.e., on 3 out of 6 datasets), we conclude that including the two parameterization choices in the hyperparameter optimization process will lead to even better performance.

**Table 5: Test accuracies averaged over $5$ training runs on UAE time-series classification datasets. All models are trained based on the same hyper-parameter tuning protocol in order to ensure fair comparability. The dataset names are abbreviations of the following UEA time-series datasets: EigenWorms (Worms), SelfRegulationSCP1 (SCP1), SelfRegulationSCP2 (SCP2), EthanolConcentration (Ethanol), Heartbeat, MotorImagery (Motor).**
|  | Worms | SCP1 | SCP2 | Ethanol | Heartbeat | Motor | Avg |
| --- | --- | --- | --- | --- | --- | --- | --- |
| Seq. length | 17,984 | 896 | 1,152 | 1,751 | 405 | 3,000 |  |
| #Classes | 5 | 2 | 2 | 4 | 2 | 2 |  |
| NRDE | $77.2\pm 7.1$ | $76.7\pm 5.6$ | $48.1\pm 11.4$ | $31.4\pm 4.5$ | $73.9\pm 2.6$ | $54.0\pm 7.8$ | $60.2$ |
| NCDE | $62.2\pm 2.2$ | $80.0\pm 2.0$ | $49.1\pm 6.2$ | $22.0\pm 1.0$ | $68.1\pm 5.8$ | $51.6\pm 6.2$ | $55.5$ |
| Log-NCDE | $82.8\pm 2.7$ | $82.1\pm 1.4$ | $54.0\pm 2.6$ | $35.9\pm 6.1$ | $74.2\pm 2.0$ | $57.2\pm 5.6$ | $64.4$ |
| LRU | $85.0\pm 6.2$ | $84.5\pm 4.6$ | $47.4\pm 4.0$ | $23.8\pm 2.8$ | $78.1\pm 7.6$ | $51.9\pm 8.6$ | $61.8$ |
| S5 | $83.9\pm 4.1$ | $87.1\pm 2.1$ | $55.1\pm 3.3$ | $25.6\pm 3.5$ | $73.9\pm 3.1$ | $53.0\pm 3.9$ | $63.1$ |
| Mamba | $70.9\pm 15.8$ | $80.7\pm 1.4$ | $48.2\pm 3.9$ | $27.9\pm 4.5$ | $76.2\pm 3.8$ | $47.7\pm 4.5$ | $58.6$ |
| S6 | $85.0\pm 16.1$ | $82.8\pm 2.7$ | $49.9\pm 9.4$ | $26.4\pm 6.4$ | $76.5\pm 8.3$ | $51.3\pm 4.7$ | $62.0$ |
| LinOSS-IM (ReLU) | $95.0\pm 4.4$ | $87.8\pm 2.6$ | $58.2\pm 6.9$ | $29.9\pm 0.6$ | $75.8\pm 3.7$ | $60.0\pm 7.5$ | $67.8$ |
| LinOSS-IM (squared) | $88.9\pm 2.5$ | $86.6\pm 1.8$ | $59.3\pm 7.8$ | $32.7\pm 6.2$ | $76.8\pm 2.2$ | $58.2\pm 8.4$ | $67.1$ |

## Appendix B Supplement to the theoretical insights

### B.1 Proof of proposition 3.2

###### Proposition .

Let $\{\lambda_{j}\}_{j=1}^{2m}$ be the eigenspectrum of the hidden-to-hidden matrix ${\bf M}^{\text{IM}}$ of the LinOSS-IM model (4). We further initialize ${\bf A}_{kk}\sim\mathcal{U}([0,A_{\text{max}}])$ with $A_{\text{max}}>0$ for all diagonal elements $k=1,\dots,m$ of ${\bf A}$ in (2). Then, the $N$-th moment of the magnitude of the eigenvalues are given as,

$$ $\mathbb{E}(|\lambda_{j}|^{N})=\frac{({\Delta t}^{2}A_{\text{max}}+1)^{1-\frac{ N}{2}}-1}{{\Delta t}^{2}A_{\text{max}}(1-\frac{N}{2})},$ (10) $$

for all $j=1,\dots,2m$, with $k=j\text{ mod }m$.

###### Proof.

By the law of the unconscious statistician together with the identity $|\lambda_{j}|=\sqrt{{\bf S}_{k}}$ from the proof of Proposition 3.1 it follows that,

$$ $\displaystyle\mathbb{E}(|\lambda_{j}|^{N})$ $\displaystyle=\frac{1}{A_{\text{max}}}\int_{0}^{A_{\text{max}}}(1+{\Delta t}^{ 2}x)^{-\frac{N}{2}}dx=\frac{1}{{\Delta t}^{2}A_{\text{max}}}\int_{1}^{{\Delta t }^{2}A_{\text{max}}+1}u^{-\frac{N}{2}}du$ $\displaystyle=\frac{({\Delta t}^{2}A_{\text{max}}+1)^{1-\frac{N}{2}}-1}{{ \Delta t}^{2}A_{\text{max}}(1-\frac{N}{2})},$ $$

where we substituted $u={\Delta t}^{2}x+1$.
∎

###### Proposition .

###### Proof.

### B.2 Proof of Theorem 3.3

###### Theorem .

Let $\Phi:C_{0}([0,T];{\mathbb{R}}^{p})\to C_{0}([0,T];{\mathbb{R}}^{q})$ be a causal and continuous operator. Let $K\subset C_{0}([0,T];{\mathbb{R}}^{p})$ be compact. Then for any $\epsilon>0$, there exist hyperparameters $m$, $\tilde{m}$, diagonal weight matrix ${\bf A}\in{\mathbb{R}}^{m\times m}$, weights ${\bf B}\in{\mathbb{R}}^{m\times d}$, $\tilde{{\bf W}}\in{\mathbb{R}}^{\tilde{m}\times m}$, ${\bf W}\in{\mathbb{R}}^{q\times\tilde{m}}$ and bias vectors ${\bf b}\in{\mathbb{R}}^{m}$, $\tilde{{\bf b}}\in{\mathbb{R}}^{\tilde{m}}$, such that the output ${\bf z}:[0,T]\to{\mathbb{R}}^{q}$ of the LinOSS model (9) satisfies,

$$ $\sup_{t\in[0,T]}|\Phi({\bf u})(t)-{\bf z}(t)|\leq\epsilon,\quad\forall\,{\bf u }\in K.$ $$

###### Proof.

We begin by considering the simple forced harmonic oscillator,

$$ ${\bf y}^{\prime\prime}(t)=-A^{2}{\bf y}(t)+{\bf u}(t),$ (11) $$

where ${\bf u}(t)\in{\mathbb{R}}^{p}$ is an external forcing and we assume $A\neq 0$. We further introduce the time-windowed sine transform for the input signal ${\bf u}(t)$,

$$ ${\mathcal{L}}_{t}{\bf u}(A)=\int_{0}^{t}{\bf u}(t-\tau)\sin(A\tau)d\tau.$ $$

Then, a straightforward calculation shows that the solution ${\bf y}(t)$ to (11) computes (up to a constant) a time-windowed sine transform, i.e.,

$$ $y(t)=A^{-1}\int_{0}^{t}{\bf u}(\tau)\sin(A(t-\tau))d\tau.$ (12) $$

This can easily be verified by differentiating ${\bf y}$,

$$ $\displaystyle{\bf y}^{\prime}(t)$ $\displaystyle=\int_{0}^{t}{\bf u}(\tau)\cos(A(t-\tau))d\tau+A^{-1}[{\bf u}( \tau)\sin(A(t-\tau))]_{\tau=t}$ $\displaystyle=\int_{0}^{t}{\bf u}(\tau)\cos(A(t-\tau))d\tau.$ $$

Differentiating one more time yields,

$$ $\displaystyle{\bf y}^{\prime\prime}(t)$ $\displaystyle=-A\int_{0}^{t}{\bf u}(\tau)\sin(A(t-\tau))d\tau+[{\bf u}(\tau) \cos(A(t-\tau))]_{\tau=t}$ $\displaystyle=-A\int_{0}^{t}{\bf u}(\tau)\sin(A(t-\tau))d\tau+{\bf u}(t)$ $\displaystyle=-A^{2}{\bf y}(t)+{\bf u}(t).$ $$

We will now make use of the fundamental Lemma in Lanthaler et al. () that provides a finite-dimensional encoding of the operator $\Phi$ we wish to approximate, i.e., for any $\epsilon>0$, there exists $N\in\mathbb{N}$, weights $A_{1},\dots,A_{N}$ and a continuous mapping $\Psi:{\mathbb{R}}^{p\times N}\times[0,T^{2}/4]\to{\mathbb{R}}^{q}$, such that

$$ $|\Phi({\bf u})(t)-\Psi({\mathcal{L}}_{t}{{\bf u}}(A_{1}),\dots,{\mathcal{L}}_{ t}{{\bf u}}(A_{N});t^{2}/4)|\leq\epsilon,$ $$

for all ${\bf u}\in K$ and $t\in[0,T]$.

Using the result in (12) we can then construct the input vector $[{\mathcal{L}}_{t}{{\bf u}}(A_{1}),\dots,{\mathcal{L}}_{t}{{\bf u}}(A_{N}),t^{
2}/4]^{\top}\in{\mathbb{R}}^{pN}\times[0,T^{2}/4]$ based on the ODE system (1) underlying our LinOSS model:

$$ $\begin{bmatrix}{\mathcal{L}}_{t}{{\bf u}}(A_{1})\\ {\mathcal{L}}_{t}{{\bf u}}(A_{2})\\ \vdots\\ {\mathcal{L}}_{t}{{\bf u}}(A_{N})\\ t^{2}/4\end{bmatrix}=\tilde{{\bf A}}{\bf y}(t),$ (13) $$

where ${\bf y}(t)\in{\mathbb{R}}^{pN}$ solves the system,

$$ ${\bf y}^{\prime\prime}(t)=-{\bf A}^{2}{\bf y}+{\bf B}{\bf u}(t)+{\bf b},$ (14) $$

with

$$ $\displaystyle{\bf A}$ $\displaystyle=\operatorname{diag}([\underbrace{A_{1},\dots,A_{1}}_{p-\text{ times}},\underbrace{A_{2},\dots,A_{2}}_{p-\text{times}},\dots\dots,\underbrace {A_{N},\dots,A_{N}}_{p-\text{times}},0]),$ $\displaystyle{\bf B}$ $\displaystyle=[\underbrace{\mathrm{I}_{p},\dots,\mathrm{I}_{p}}_{N-\text{times }},0]^{\top},\quad{\bf b}=[0,\dots,0,1/2]^{\top}$ $$

where $\mathrm{I}_{p}\in{\mathbb{R}}^{p\times p}$ is the identity matrix and $\tilde{{\bf A}}$ equals ${\bf A}$, except that $\tilde{{\bf A}}_{pN,pN}=1$. Thus, $\tilde{{\bf A}}{\bf y}(t)$ computes exactly the input to the finite-dimensional operator $\Psi$.
By the universal approximation theorem for ordinary neural networks there exist weight matrices ${\bf W},\hat{{\bf W}}$ and bias $\tilde{{\bf b}}$, such that,

$$ $|\Psi({\mathcal{L}}_{t}{{\bf u}}(A_{1}),\dots,{\mathcal{L}}_{t}{{\bf u}}(A_{N} );t^{2}/4)-{\bf W}\sigma(\hat{{\bf W}}[{\mathcal{L}}_{t}{{\bf u}}(A_{1}),\dots ,{\mathcal{L}}_{t}{{\bf u}}(A_{N}),t^{2}/4]^{\top}+\tilde{{\bf b}})|<\epsilon.$ $$

Thus, for every ${\bf u}\in K$ we have,

$$ $\displaystyle|\Phi({\bf u}(t))-{\bf z}(t)|\leq$ $\displaystyle|\Phi({\bf u}(t))-\Psi({\mathcal{L}}_{t}{{\bf u}}(A_{1}),\dots,{ \mathcal{L}}_{t}{{\bf u}}(A_{N});t^{2}/4)|$ $\displaystyle+|\Psi({\mathcal{L}}_{t}{{\bf u}}(A_{1}),\dots,{\mathcal{L}}_{t}{ {\bf u}}(A_{N});t^{2}/4)-{\bf W}\sigma(\tilde{{\bf W}}{\bf y}(t)+\tilde{{\bf b }})|$ $\displaystyle<2\epsilon,$ $$

with $\tilde{{\bf W}}=\hat{{\bf W}}\tilde{{\bf A}}$.
Since $\epsilon>0$ was arbitrary, we conclude that for any causal and continuous operator $\Phi:C_{0}([0,T];{\mathbb{R}}^{p})\to C_{0}([0,T];{\mathbb{R}}^{q})$, compact set $K\subset C_{0}([0,T];{\mathbb{R}}^{p})$ and $\epsilon>0$, there exists a LinOSS model of form (9), which uniformly approximates $\Phi$ to accuracy $\epsilon$ for all ${\bf u}\in K$. This completes the proof.

∎

###### Theorem .

###### Proof.

## Peer Reviews

confidence: 5
contribution: 3
initial_rating: 5
presentation: 4
questions: - What are the runtimes of the different methods in the experiments? While the accuracy metrics are strong, one of the purported examples of SSMs is the efficiency. What are run times when enabling and disabling different introduced optimizations:
  - Run time when not using the equation for matrix inversion?
  - Run time when not using the parallel scan?
- Line 74: Why is A diagonal? Is it only to induce the matrix inversion trick later?
- Some aspects of the architecture are unclear. How exactly are the LinOSS layers stacked into each other and within the network? 
  - In Figure 1, which parts of  the figure are u, y, z?
  - What effect do the nonlinearities have on the stability analysis?
  - How is the nonlinear layer in the figure a part of the model?
- In the experiments in Section 4, what are the exact hyperparameters and model graphs?
- The exploitation of Formula 3 to speed up the matrix is clever. How does this differentiate during training, though? Did you pass this formula through autodiff, or define a special differentiation rule?
- Does the parallel scans apply to only using first-order IM or IMEX integrators? Or does the ODE formula in general allow for the parallel scans with other integrators, such as a simple forward Euler, or a higher order IMEX?
- Table 1: s/UAE/UEA/g.
- Table 1: What does UEA stand for? Define abbreviations and add citations. (Is Walker 2024 supposed to be the citation?)
- Color highlighting in the tables is not colorblind friendly, nor BW printer/ereader friendly. Use symbols instead.
- Equation 7: Why is it important that equation 2 is a Hamiltonian system in this section? If the underlying ODE is indeed Hamiltonian, doesn’t that suggest that the IM discretization is not appropriate?
- Line 222: What are a, b supposed to be? Describe the specific case of operator & tuple shown here.
- Line 280: Why is it possible to assume that A_kk>0? Couldn’t they diverge from that assumption during training? What initialization is required?
- Line 289: s/constraint/constrained/
- Line 292: The steps in the proof are not obvious. More steps of proof should be presented. In the appendix would be sufficient.
- Just a remark: The title of section 4.2 is overly extravagant for the claims. These days, “extreme long ranges” would be context lengths of millions of tokens :)
- Appendix A: What are the parameter counts? What is the architecture of the “nonlinear layers”? What are the hyperparameters for the forecasting problem?
- Why does the PPG model have higher memory demands, when the models actually seem small?
- Appendix A: What ML library was used?
- Could you describe the loss functions and the complete architecture for the 3 different types of problems? How is time series classification grafted onto the LinOSS network?
- The importance of section 3.2 is unclear to me. It is nice to prove universality of a model, but what is special about LinOSS that other more general proofs of universality would not apply? I would not have questioned it. Is any aspect of the proof specific to LinOSS, or could it apply more broadly to more SSMs?
soundness: 3
strengths: The paper demonstrates a cross cutting expertise from dynamical systems analysis through implementation optimization. The design of the ODE is crafting three advantages at three different levels of abstraction simultaneously: enforce theoretically proven stabilization, allow for efficient matrix inversion, and allow for parallel scans of the sequence recurrence. The experiments on time series problems are a good set of problems, and the results of LinOSS stand out against the broad set of comparisons. The proofs are sound, however, it is not clear if they apply to what actually is going on: see below.
summary: This paper introduces a new continuous time recurrent network in the family of state space models. The architecture is proposed as an ODE that is discretized in two ways using first order implicit and implicit-explicit integrators. The terms of the ODE are further constrained to induce two different computational tricks to speed up computation; fast matrix inversion to make the implicit methods tractable, and parallel scans for faster sequence processing. The LinOSS architecture is empirically compared to other Neural ODEs and state-space models, and to transformers, on three different time domain problems for time series classification and prediction.
weaknesses: One issue with the paper is highlighted in the core claim of pre hoc controlling for “forgetting” versus stability by choosing between LinOSS-IM and IMEX (Line 298).  Line 205 claims to demonstrate different advantages between the two methods, but this is not actually evident in the experiments. What characteristics of the problems in Table 1 lead to IM vs. IMEX performing differently? If anything, there is no difference between IM and IMEX in all examples but Worms. Is there something special about Worms, or is it a fluke? The result of Table 2 seems to contradict the claim that IMEX is better at long ranges and IM is better at forgetting: why does IM perform better on a problem where memory over a long sequence should be important? 

It is odd to rely on the integrator choice to enforce stability vs. forgetting. Using an ODE framework, it would seem more natural to change the ODE with dissipative terms, for example, instead of changing the integrator. By relying on the integrator choice, it is unclear if those properties would actually hold after training, if the discretization is thought to add new properties that the ODE did not have. Given the unclear results of the experiments, is it possible that results of the theoretical analysis do not apply after the model is trained? Did you try inspecting if the parameters of A and S still hold the assumptions / initializations assumed in Section 3, after training?


Consider the case where the system being learned is actually stable, but the backward euler IM formula is applied. Backward Euler is dissipative even when the dynamical system does not want to dissipate. What happens when you try to learn a model that should be energy conserving using  the LinOSS-IM architecture? The authors could try this by just trying to learn to forecast a simple oscillatory system with LinOSS-IM. I would expect that the model would learn “through” the IM discretization, and converge to an parameterization of an “unstable” ODE that is stable after being discretized by IM. See Krishnapriyan, “Learning continuous models for continuous physics” for a discussion on overfitting on learning through ODE discretizations. 

One idea to start to tackle this problem: Try LinOSS-FE “forward euler”. The architecture would be similar and the tricks in parallel scans would still work. This would be an ablation that would illustrate why stabilization of the implicit and imex integrators is important. If there is no performance difference with the Forward euler discretization in the experiments, or if the IM method can forecast a stable system, then perhaps the theoretical results do not represent what actually happens after training.

---

confidence: 2
contribution: 4
initial_rating: 8
presentation: 3
questions: To the best of my understanding, the model cannot produce chaotic dynamics. What if the task in hand requires this? How does this contrast (if there is a contrast) with section 3.2?
soundness: 4
strengths: Provides strong theoretical together with intuitive explanations.
Contrasts their two proposed methods mathematically and also experimentally.
The experimental results are excellent and definitely contributes to the field significantly.
Supplementary material is comprehensive.
summary: The paper provides a novel state-space model. They have two different versions of it, in which they rigorously show the power of their algorithms and also experimentally verify it. The method outperforms SOTA methods in many tasks.
weaknesses: I think section 3.2 can be written more accessible.

I believe Figure 1 is very important but can be made more explanatory.

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: 1. Could the model be extended to incorporate coupled oscillations while maintaining stability guarantees?

2. What is the impact of the time step parameter Δt on model performance and stability?

3. How does the choice between IM and IMEX variants affect training dynamics and convergence?
soundness: 3
strengths: The paper demonstrates strong theoretical foundations by providing rigorous mathematical analysis of stability conditions, proving universal approximation capability, and establishing clear connections to Hamiltonian systems and symplectic integration. Implementation-wise, it offers a remarkably simple parameterization requiring only non-negative diagonal elements, achieves efficient computation through parallel scans, and presents two complementary variants with different preservation properties. The empirical results are great, showing strong performance on diverse tasks.
summary: This paper introduces Linear Oscillatory State-Space models (LinOSS), a novel approach to sequence modeling based on forced harmonic oscillators. The model comes in two variants: LinOSS-IM (implicit) and LinOSS-IMEX (implicit-explicit). The key innovation is using second-order ODEs with diagonal state matrices, offering stable dynamics while only requiring non-negative diagonal elements. The paper provides theoretical guarantees for stability and universal approximation, along with empirical validation showing significant improvements over state-of-the-art models like Mamba and LRU on long sequences.
weaknesses: 1. Limited Analysis of Model Interpretability:

While based on oscillatory dynamics, lacks discussion of learned frequencies

No analysis of how the model captures different timescales

2. Experiment:

No ablation studies on the impact of different initialization schemes

The implementation details are unclear. For instance, how does it compare to Mamba or S5 in terms of speed, training time, FLOPs, and memory usage? Discussing these aspects could enhance its practical utility.

---

confidence: 3
contribution: 2
initial_rating: 6
presentation: 3
questions: * L116-118. Since the oscillators are independent, can the proposed architecture model transient synchronization/desynchronization?
* L250. `Assuming M is diagonalizable [...]`. If $M$ has real eigenvalues with algebraic multiplicity $ > 1$, would that make the system unstable? Admittedly, norm growth would be sub-exponential, so perhaps in practice it's fine.
* L289 `constraint` $=>$ constrained
soundness: 3
strengths: The work is well motivated and the writing is clear and concise. The recurrence matrix afforded by the proposed approach has desirable stability properties and is less constrained compared to the typical SSM parametrizations. The expressivity is further supported by the theoretical analysis on universality of the proposed parametrization. The experimental results are thorough and demonstrate the efficacy of the proposed approach relative to baselines.
summary: In this work, the authors propose a state-space model (SSM) architecture, Linear-Oscillatory SSM, derived from the discretization of second-order linear ODEs that models a network of forced harmonic oscillators. They develop two discretization schemes for the approach, study the stability properties of the parametrization and show universal approximation of LinOSS for approximating general continuous, causal input-output mapping. Through empirical evaluations, they demonstrate that their approach outperforms other Linear SSMs on time series classification, prediction and long-horizon forecasting.
weaknesses: Given that the parametrization has been introduced in Rusch & Mishra (2021) and the universality results for the non-linear counterpart to this work have been shown in Lanthaler et al., 2024, the novelty of this work is somewhat limited. Still, I think that this is a useful contribution overall as it improves the ability of SSMs for learning long-range dependencies.

Meta-review: It is a pleasure when the reviewers all converge on a consistent recommendation. I see no reason to disagree with the universal praise for the paper, and strongly recommend acceptance, and believe the paper should be highlighted with a spotlight presentation.

Decision: Accept (Oral)