HOPE for a Robust Parameterization of Long-memory State Space Models

## Abstract

Abstract State-space models (SSMs) that utilize linear, time-invariant (LTI) systems are known for their effectiveness in learning long sequences. To achieve state-of-the-art performance, an SSM often needs a specifically designed initialization, and the training of state matrices is on a logarithmic scale with a very small learning rate.
To understand these choices from a unified perspective, we view SSMs through the lens of Hankel operator theory.
Building upon it, we develop a new parameterization scheme, called HOPE, for LTI systems that utilizes Markov parameters within Hankel operators.
Our approach helps improve the initialization and training stability, leading to a more robust parameterization.
We efficiently implement these innovations by nonuniformly sampling the transfer functions of LTI systems, and they require fewer parameters compared to canonical SSMs. When benchmarked against HiPPO-initialized models such as S4 and S4D, an SSM parameterized by Hankel operators demonstrates improved performance on Long-Range Arena (LRA) tasks. Moreover, our new parameterization endows the SSM with non-decaying memory within a fixed time window, which is empirically corroborated by a sequential CIFAR-10 task with padded noise.

## 1 Introduction

State-space models (SSMs) [] have gained popularity and success in sequence modeling.
Known for its excellent efficiency and capability of handling long sequences, an SSM leverages the continuous-time linear, time-invariant (LTI) systems.
These systems are often defined by four matrices $\Gamma=(\mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D})$ as

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

and they can be used to model the mappings from input time-series $\mathbf{u}(\cdot)$ to the output times-series $\mathbf{y}(\cdot)$, where $\mathbf{u}(t)\in\mathbb{R}^{m}$ and $\mathbf{y}(t)\in\mathbb{R}^{p}$ for every $t$.
The (hidden) states, which capture the latent dynamics, are denoted as $\mathbf{x}=\mathbf{x}(t)\in\mathbb{R}^{n}$.
The system matrices are of dimensions $\mathbf{A}\in\mathbb{C}^{n\times n}$, $\mathbf{B}\in\mathbb{C}^{n\times m}$, $\mathbf{C}\in\mathbb{C}^{p\times n}$, and $\mathbf{D}\in\mathbb{C}^{p\times m}$.
Often, the size $n$ of the state vector $\mathbf{x}$ is much larger than $m$ and $p$, which allows us to memorize information about the past inputs $\mathbf{u}|_{(-\infty,t]}$ in the state vector $\mathbf{x}(t)$ and retrieve it later to compute $\mathbf{y}$ via $\mathbf{C}$.

The so-called S4 [] and S4D [] models both set $m=p=1$, and they differ in the structural requirement of $\mathbf{A}$.
This framework was later generalized to the case where $m,p>1$ by the S5 model [] via the parallel scans.
Another line of research involves making the state transition rule $\mathbf{A}$ depend on the input $\mathbf{u}$, along which the two most notable models are Liquid-S4 [] and Mamba [], where the latter model achieves the state-of-the-art performance on large-scale real-world datasets.

However, SSMs typically need to be initialized and trained (very) carefully.
A randomly initialized SSM has suboptimal performance, but the so-called high-order polynomial projection operators (HiPPOs) [, , ] can be used to empirically improve it.
On the other hand, even a properly initialized SSM needs to be trained with care.
One often needs to set a smaller learning rate for the matrix $\mathbf{A}$ [], and the LTI systems require reparameterization to be trained stably []. To better understand these initialization and reparameterization efforts from a unified perspective, we analyze SSMs through the lens of Henkel Operator theory.
Specifically, we use the Hankel singular value decomposition (HSVD) to analyze an operator defined by $\Gamma$.
The decay of the Hankel singular values tells how “expressive” the LTI system is.
If the Hankel singular values of $\Gamma$ decay fast, then it informally means that our $\Gamma$ cannot capture the complex patterns in the input space $\{\mathbf{u}(\cdot)\}$; in fact, the theory of reduced-order modeling (ROM) says that $\Gamma$ can be well-approximated by a reduced system with a much smaller state-space dimension $k\ll n$ [].

We find that the decay of the Hankel singular values can be used for predicting the performance of an SSM, and that every previous effort in proposing a good initialization and training scheme can be viewed as an effort to avoid fast-decaying Hankel singular values.
This is reminiscent of the works by [] and [] that connect the singular values of the weight matrices to a deep neural network’s performance. Using the Hankel singular values as heuristics, we show that an S4D model is vulnerable to losing expressiveness during training. Moreover, even with a reparameterization, an LTI system is very sensitive to a perturbation of its parameters, $\mathbf{A}$, $\mathbf{B}$, and $\mathbf{C}$, impairing the training stability of an S4D model.

Figure: Figure 1: There are many equivalent ways to represent an LTI system. While most of the canonical SSMs use continuous LTI systems as their parameters, we propose to parameterize an SSM by the Markov parameters in its Hankel operator. The feedthrough matrix $\mathbf{D}$ is not shown in the diagram, but it is also a parameter of the LTI layers in both the canonical SSMs and our HOPE-SSM.
Refer to caption: x1.png

Based on these insights, we propose a completely different parameterization of the LTI systems. Instead of parameterizing the continuous-time systems by matrices $\mathbf{A}$, $\mathbf{B}$ and $\mathbf{C}$, we parameterize our LTI systems by the Markov parameters of the so-called discrete Hankel operator (HOPE).
A discrete Hankel operator is defined by a doubly infinite Hankel matrix, and is naturally associated with a discrete LTI system, and with a continuous-time system via the bilinear transform with $\Delta t=1$ [].
While this continuous-time LTI system acts on our sequential data, the optimization algorithms are applied to the Markov parameters of the Hankel matrix. (See Figure 1.)
We prove that unlike an LTI system parameterized by $(\mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D})$, one parameterized by the Markov parameters almost surely has slowly decaying singular values (see 3); moreover, it enjoys a global stability to perturbation (see 4). Hence, unlike a canonical SSM, our HOPE-SSM can be stably trained without reparameterization or reducing the learning rate, also reducing the need for hyperparameter tuning. We show that our HOPE-SSM can be implemented by nonuniformly sampling the transfer function, which shares the same computational complexity as the S4D model.
Moreover, it requires only $1/3$ the number of parameters of an LTI system in an S4D model to parameterize that in a HOPE-SSM.

The practical benefits of our novel parameterization are improved robustness (with respect to initialization and training stability) and performance (with respect to model quality and parameter count).
Moreover, we show that the memory of our HOPE-SSM does not decay (see eq. 8) in a fixed time window, making it possible to solve tasks that involve even longer-range dependency by tuning the sampling period $\Delta t$ at discretization.
This partially addresses the well-known issue that the LTI system of a canonical SSM suffers from exponentially decaying memory [].

Related Work. Initially proposed by [], the general idea of the HiPPO framework is to memorize the input by projecting it onto an orthogonal polynomial basis and storing the coefficients in the state vector $\mathbf{x}(t)$. This was later generalized to some different orthogonal polynomial bases in [, ]. The S4D model uses a slightly perturbed version of the HiPPO-LegS initialization, and the effect of the perturbation was studied in []. The initialization issue was also studied in [] in the discrete-time setting, which provides an alternative justification of HiPPO based on the spectrum of the state matrix. While a common way to reparameterize a diagonal SSM is by training $\text{Re}(\text{diag}(\mathbf{A}))$ on a logarithmic scale [], other stable reparameterizations were considered in []. A method that directly parameterizes the convolutional kernel of the discretized LTI system was presented in []. Compared to their work, in this paper, we adopt a compute-then-discretize strategy by still parameterizing the underlying continuous dynamics, making our SSM capable of handling sequences with varying lengths.

The decaying memory of RNNs and SSMs is analyzed and discussed by a wide literature [, , , ]. Ways to lift the memory capacity of an SSM were considered in [] and [] via either reparameterizing the state matrix or applying a spectral filter. We remark that while [] showed that a spectral SSM is exponentially close to an SSM, our HOPE-SSM is an actual SSM containing LTI systems, and it only differs from a canonical SSM in the parameterization.

Contributions.
Our main contributions are summarized as follows.

- 1.
We show that high-degree LTI systems (i.e., those with slow-decaying Hankel singular values) in an SSM lead to a good performance. We justify this causal relationship using ideas in reduced-order modeling (ROM).
We theoretically prove that expressive, high-degree LTI systems are scarce in the parameter space of $(\mathbf{A}$, $\mathbf{B},\mathbf{C})$.
Thus, they need careful designs and are numerically unstable.
This explains difficulties in initializing and training SSMs.
- 2.
We propose a new parameterization of LTI systems using the Hankel operator (HOPE), which can be implemented efficiently by nonuniformly sampling the transfer function and requires $1/3$ the number of parameters in an LTI system from an S4D model. We prove that our new parameter space is full of high-degree LTI systems.
Hence, our HOPE-SSM does not suffer from the lack of expressiveness during initialization and training; moreover, it can be stably trained and has non-decaying long memory.
- 3.
We empirically demonstrate that our HOPE-SSM is robust using the sCIFAR-10 task and that it has long-term memory using a noise-padded sCIFAR-10 dataset.
We test the performance of a full-scale HOPE-SSM on the Long-Range Arena and observe that its performance exceeds that of its S4 and S4D counterparts and many other SSMs for most tasks.

## 2 Preliminaries

Let $\Gamma=(\mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D})$ be a continuous-time LTI system defined in eq. 1. One can take a bilinear transform to obtain a discrete system $\overline{\Gamma}=(\overline{\mathbf{A}},\overline{\mathbf{B}},\overline{
\mathbf{C}},\overline{\mathbf{D}})$ so that the underlying dynamics is given by

$$ $\displaystyle{\mathbf{x}}_{k+1}=\overline{\mathbf{A}}{\mathbf{x}}_{k}+ \overline{\mathbf{B}}{\mathbf{u}}_{k},\qquad{\mathbf{y}}_{k}=\overline{\mathbf {C}}{\mathbf{x}}_{k}+\overline{\mathbf{D}}{\mathbf{u}}_{k}.$ (2) $$

The transfer functions of $\Gamma$ and $\overline{\Gamma}$ are rational functions on the complex plane $\mathbb{C}$ defined by

$$ $G(s)=\mathbf{C}(s\mathbf{I}-\mathbf{A})^{-1}\mathbf{B}+\mathbf{D}\qquad\text{ and}\qquad\overline{G}(z)=\overline{\mathbf{C}}(z{\mathbf{I}}-\overline{ \mathbf{A}})^{-1}\overline{\mathbf{B}}+\overline{\mathbf{D}},$ $$

respectively.
We usually care about the values of $G$ and $\overline{G}$ on the imaginary axis and the unit circle in $\mathbb{C}$, respectively. They “transfer” the inputs to the outputs in the frequency domain by multiplication:

$$ (continuous case) $\displaystyle\qquad\hat{\mathbf{y}}(s)=G(is)\hat{\mathbf{u}}(s),\qquad s\in \mathbb{R},$ (3) (discrete case) $\displaystyle\qquad\hat{{\mathbf{y}}}_{k}=\overline{G}({\bm{\omega}}_{k})\hat{ {\mathbf{u}}}_{k},\qquad 0\leq k\leq L-1,$ $$

where the hats on a function and on a vector mean the Fourier transform and the Fourier coefficients, respectively, and ${\bm{\omega}}_{k}=\text{exp}(i2\pi k/L)$ is the $k$th Fourier node of length $L$. Throughout this paper, we assume that our LTI systems are asymptotically stable, i.e., the eigenvalues of $\mathbf{A}$ have negative real parts, or equivalently, the eigenvalues of $\overline{\mathbf{A}}$ are contained in the open unit disk in the complex plane. In this paper, we also discuss a completely different notion of stability: the numerical stability of an LTI system. This refers to the system’s sensitivity to a perturbation of its parameters. We will clearly distinguish these two notions of stability by appending the adjectives “asymptotic” or “numerical”.

Given a discrete LTI system $\overline{\Gamma}$, its Hankel operator is defined by a doubly infinite Hankel matrix

$$ $\overline{\mathbf{H}}\in\mathbb{C}^{\infty\times\infty},\qquad\overline{ \mathbf{H}}:\ell^{2}(\mathbb{N})\rightarrow\ell^{2}(\mathbb{N}),\qquad \overline{\mathbf{H}}_{i,j}=\overline{\mathbf{C}}\overline{\mathbf{A}}^{i+j} \overline{\mathbf{B}},\qquad i,j\geq 0.$ (4) $$

This discrete Hankel operator on $\ell^{2}(\mathbb{N})$ is a bounded linear operator of rank $\leq n$, the number of latent states. We denote its singular values by $\sigma_{1}(\overline{\mathbf{H}})\geq\sigma_{2}(\overline{\mathbf{H}})\geq
\cdots\geq\sigma_{n}(\overline{\mathbf{H}})\geq 0$.
Analogously, one can define a continuous Hankel operator $\mathbf{H}:L^{2}([0,\infty))\rightarrow L^{2}([0,\infty))$ given a continuous-time LTI system $\Gamma$. If $\overline{\Gamma}$ is discretized from $\Gamma$ using the bilinear transform, then the singular values of $\mathbf{H}$ are equivalent to those of $\overline{\mathbf{H}}$, i.e., $\sigma_{j}(\mathbf{H})=\sigma_{j}(\overline{\mathbf{H}})$. (See Appendix B for more details.) In this paper, we consider the decay of the Hankel singular values. To quantitatively measure how “small” a singular value is, we introduce the numerical rank of an LTI system. Given a small number $\epsilon\geq 0$, we define the $\epsilon$-rank of $\Gamma$ to be the number of relative Hankel singular values larger than $\epsilon$:

$$ $\text{rank}_{\epsilon}(\Gamma)=\max\{j\mid\sigma_{j}(\mathbf{H})/\sigma_{1}( \mathbf{H})>\epsilon\}.$ $$

When $\epsilon=0$, we reproduce the exact rank of $\mathbf{H}$, which is rarely $<n$ in the floating-point arithmetic. A small positive $\epsilon$ allows us to eliminate the small but perhaps positive singular values.

## 3 Unravel a Mystery: Hankel Singular Values in Initialization and Training

Figure: Figure 2: Test accuracy of the SSMs on the sCIFAR task. The LTI systems are initialized in three different ways and are either trained or untrained. We notice that when the LTI systems are initialized with $\texttt{init}_{1}$ (red), training the LTI system together with other model parameters is impairing the model accuracy. This is in contrast to SSMs initialized with $\texttt{init}_{3}$ (blue), where assigning the LTI system a small positive learning rate is helping the performance.

We begin our exploration by presenting a mystery. We use it to elicit and support the use of the Hankel singular values as protocols for explaining and predicting the success of SSMs, without stressing it as a fine-grained study of different initializations. We train an S4D model to learn the sCIFAR-10 image classification task [, ]. We consider three initialization schemes of the LTI systems in the SSM: $\texttt{init}_{1},\texttt{init}_{2}$, and $\texttt{init}_{3}$. While $\texttt{init}_{3}$ is the HiPPO-LegS initialization, we treat the others as black boxes in this paper, leaving the details to Appendix C. Instead, we later explain their success or failure by measuring their Hankel singular values. For an SSM initialized using a certain $\texttt{init}_{j}$ ($1\leq j\leq 3$), we train it in two different ways: either by freezing $\mathbf{A},\mathbf{B}$, and $\mathbf{C}$ and only training the other model parameters or by assigning $\mathbf{A},\mathbf{B}$, and $\mathbf{C}$ a small learning rate. The three initializations and two learning rate assignments comprise a total of six combinations, summarized in Figure 2. Besides the natural question of why we see a general difference between models initialized differently, Figure 2 raises a more intriguing mystery:
{adjustwidth}15pt15pt
For an SSM initialized by $\texttt{init}_{1},\texttt{init}_{2}$, or $\texttt{init}_{3}$, why does training the LTI systems impair, level, or improve the performance of the model, respectively?
To answer these questions, we study the Hankel singular values of the systems, but why are they relevant? The reason is that the Hankel singular values measure the “complexity” of an LTI system. The easiest way to see this is via ROM of the system $\Gamma$ [, ]:
{adjustwidth}15pt15pt
For any $k<n$, there exists a reduced-order approximation $\tilde{\Gamma}=(\tilde{\mathbf{A}},\tilde{\mathbf{B}},\tilde{\mathbf{C}},
\tilde{\mathbf{D}})$ with $\tilde{\mathbf{A}}\in\mathbb{C}^{k\times k}$, such that $\|G-\tilde{G}\|_{\infty}\leq\sum_{j=k+1}^{n}\sigma_{j}(\mathbf{H})\leq(n-k)
\sigma_{k+1}(\mathbf{H})$, where $\tilde{G}$ is the transfer function of $\tilde{\Gamma}$ and $\|\cdot\|_{\infty}$ is the $L$-infinity norm over the imaginal axis.
In particular, if the sum of the trailing Hankel singular values $\sum_{j=k+1}^{n}\sigma_{j}(\mathbf{H})$ is small, then by eq. 3, $\Gamma$ and $\tilde{\Gamma}$ produce similar outputs on any input.
Hence, ROM says that if the Hankel singular values decay fast, then we can replace the LTI system with a much smaller one without too much loss; in other words, most states in $\mathbf{x}(t)$ are not properly used to memorize the input.

Figure: Figure 3: The distribution of all relative Hankel singular values $\sigma_{j}(\mathbf{H})/\sigma_{1}(\mathbf{H})$ of the LTI systems in an SSM. For each initialization, the distribution is shown both at initialization and after the SSM is trained for $10$ epochs. Note that the second row only applies when the LTI systems are not frozen.

In our experiment, each SSM has $4$ layers and $128$ channels. These comprise $4\times 128=512$ different copies of single-input/single-output LTI systems.
Every LTI has $n=64$ latent states. These make up $512\times 64=32768$ relative Hankel singular values $\sigma_{j}(\mathbf{H})/\sigma_{1}(\mathbf{H})$ to consider. In Figure 3, we show the histograms of all these relative Hankel singular values at two different stages of training. Note that the three histograms on the second row only apply when the LTI systems are trained with a small learning rate; when they are frozen, the distributions always stay the same as those at initialization. We can explain the behaviors of the three SSMs using their Hankel singular values:

- 1.
The systems in a model initialized by $\texttt{init}_{1}$ initially have high numerical ranks. Hence, when the systems are untrained, they define random mappings that capture the rich content in the input data, yielding the rest of the work to other model parameters in the SSM. However, when the systems are trained, their Hankel singular values start to decay rapidly with only $27.87\%$ of the singular values satisfying that $\sigma_{j}(\mathbf{H})/\sigma_{1}(\mathbf{H})>0.01$, and the systems can no longer handle a variety of distinct inputs. This makes it harder to parse complicated images in the sCIFAR-10 dataset.
- 2.
No matter whether trained or not, the Hankel singular values of the systems in the model initialized by $\texttt{init}_{2}$ decay very fast, which means the systems essentially lack expressiveness and cannot capture the complicated patterns in the input space. Hence, the SSMs with both trained and untrained LTI systems could not learn the task effectively.
- 3.
The systems in the model initialized by $\texttt{init}_{3}$, trained or not, have high numerical ranks. In particular, over $87.82\%$ of the singular values satisfy that $\sigma_{j}(\mathbf{H})/\sigma_{1}(\mathbf{H})>0.01$ after $10$ epochs.
In this case, training the LTI systems allows us to accelerate the optimization.

We just established the Hankel singular values as a protocol for evaluating and predicting the performance of an SSM on tasks that involve complicated long-range dependencies. We can further derive theory to explain the patterns we observed in the histograms in Figure 3, which brings out potential weaknesses of the S4D models and motivates the development of our HOPE-SSM.

### 3.1 Many LTI Systems Have Low Ranks

We saw that LTI systems with high numerical ranks make an SSM thrive, but do we have an abundance of them? We approach this question from a random matrix theory perspective. We randomly sample a discrete LTI system. We assume that every diagonal entry $a_{j}$ of $\overline{\mathbf{A}}=\text{diag}(a_{1},\ldots,a_{n})$ are sampled i.i.d. from a distribution $F_{{a}}$ with

$$ $\mathbb{P}\left(\left|a_{j}\right|>(1-\rho)\right)=\mathcal{O}(\rho^{\alpha}) \qquad\text{ as }\qquad\rho\rightarrow 0^{+}$ $$

for some $\alpha>0$. For example, if $a_{j}$ is uniformly distributed on the open unit disk, then we have $\alpha=1$. Moreover, assume that $\overline{\mathbf{B}}\circ\overline{\mathbf{C}}^{\top}$ is a random vector with each entry sampled i.i.d. from a normal distribution $\mathcal{N}(0,1)$, where $\circ$ is the Hadamard product. The skip connection matrix $\overline{\mathbf{D}}$ has no effect on the Hankel singular values. With these assumptions, we formally show that the system has a low $\epsilon$-rank with high probability.

###### Theorem 1 .

Given any $\epsilon>0$, $0<\alpha\leq 1$, and $0<\delta\leq 1$, with probability at least $1-\delta$, the $\epsilon$-rank of $\overline{\Gamma}=(\overline{\mathbf{A}},\overline{\mathbf{B}},\overline{
\mathbf{C}},\overline{\mathbf{D}})$ with $a_{j}\sim F_{a}$ i.i.d. and $b_{j}c_{j}\sim\mathcal{N}(0,1)$ i.i.d. is $\mathcal{O}(\ln\!\left(\delta^{-3/2}\epsilon^{-1}n\right)n^{\beta})$, where

$$ $\beta\leq\frac{1}{1+\alpha}+\frac{\ln(2+\sqrt{\ln(1/\delta)/2})}{\ln(n)}$ $$

and the constant in $\mathcal{O}$ is universal.

We defer the proof to Appendix D. What 1 says is that if we ignore the logarithmic factors in the $\mathcal{O}$-notation, then the $\epsilon$-rank of $\overline{\Gamma}$ scales like $n^{\beta}$ as $n\rightarrow\infty$, where $\beta<1$. For example, when $a_{j}$ are uniformly distributed on the unit disk, we have $\beta=1/2$ plus a small number bounded by $\ln(2+\sqrt{\ln(1/\delta)/2})/\ln(n)$. In practice, we find that this term can almost be ignored (see Appendix G). The importance of 1 is twofold: from an expository point of view, it theoretically verifies, using our Hankel operator framework, that random initializations of LTI systems could lead to poor model performance, as observed in [] and []; from a practical perspective, it suggests high-rank systems are only scarce in the space of S4D model parameters. Hence, even when an LTI system is initialized with slow-decaying Hankel singular values, when $\mathbf{A}$, $\mathbf{B}$, and $\mathbf{C}$ are perturbed during training, it is at risk of losing numerical ranks and thus expressiveness. This is indeed observed in Figure 3, even for $\texttt{init}_{3}$.

###### Theorem 1 .

### 3.2 LTI Systems are Numerically Unstable under Perturbations

In this section, we perform a sensitivity analysis of an S4D model, which suggests a numerical stability issue in training the model.

###### Theorem 2 .

Let $\Gamma=(\mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D})$ be a stable continuous-time LTI system, where $\mathbf{A}=\text{diag}(a_{1},\ldots,a_{n})$ is diagonal. Let $\tilde{\Gamma}=(\tilde{\mathbf{A}},\tilde{\mathbf{B}},\tilde{\mathbf{C}},{
\mathbf{D}})$ be a perturbed stable system with $\tilde{\mathbf{A}}=\text{diag}(\tilde{a}_{1},\ldots,\tilde{a}_{n})$. Assume there exist $\Delta_{\mathbf{A}},\Delta_{\mathbf{B}}>0$ such that $\left|a_{j}-\tilde{a}_{j}\right|\leq\Delta_{\mathbf{A}}\leq\min_{j}|\text{Re}(
a_{j})|/2$ and $|b_{j}c_{j}-\tilde{b}_{j}\tilde{c}_{j}|\leq\min(|b_{j}c_{j}|,\Delta_{\mathbf{B
}})$ for all $j=1,\ldots,n$. Let $G$ and $\tilde{G}$ be the transfer functions of $\Gamma$ and $\tilde{\Gamma}$, respectively. Then, the following statements hold.

- (a)
We have
$\|G-\tilde{G}\|_{\infty}\leq 4n\Delta_{\mathbf{A}}\max_{j}\frac{\left|b_{j}c_{
j}\right|}{\left|\text{Re}(a_{j})\right|^{2}}+n\Delta_{\mathbf{B}}\max_{j}
\frac{1}{|\text{Re}(a_{j})|}.$
- (b)
The upper bound is tight up to a factor of $n$. That is, given any $\Gamma$, $\Delta_{\mathbf{A}}\leq\min_{j}|\text{Re}(a_{j})|/2$, and $\Delta_{\mathbf{B}}\leq\min_{j}|b_{j}c_{j}|$, there exist two systems $\tilde{\Gamma}_{\mathbf{A}}$ and $\tilde{\Gamma}_{\mathbf{B}}$, with transfer functions $\tilde{G}_{\mathbf{A}}$ and $\tilde{G}_{\mathbf{B}}$, respectively, that satisfy the above perturbation conditions and have
$\|G-\tilde{G}_{\mathbf{A}}\|_{\infty}\geq\Delta_{\mathbf{A}}\max_{j}\frac{
\left|b_{j}c_{j}\right|}{\left|\text{Re}(a_{j})\right|^{2}},\qquad\|G-\tilde{G
}_{\mathbf{B}}\|_{\infty}\geq\Delta_{\mathbf{B}}\max_{j}\frac{1}{|\text{Re}(a_
{j})|}.$

The proof can be found in Appendix E. 2 says that when a diagonal LTI system is trained, its numerical stability depends on two things: the proximity of its poles $a_{j}$ to the imaginary axis and the magnitudes of $\mathbf{B}$ and $\mathbf{C}$. We defer the discussion of $a_{j}$ to the end of this subsection. One may imagine that $\|G\|_{\infty}$ is related to $|b_{j}c_{j}|$, and hence, $|b_{j}c_{j}|$ plays no role in the relative numerical stability (i.e., $\|G-\tilde{G}\|_{\infty}/\|G\|_{\infty}$) of the system. However, this is not true as a system with a small $\|G\|_{\infty}$ may have an arbitrarily large $\mathbf{B}\circ\mathbf{C}^{\top}$ (see Figure 4 (Left)). Working with such a system can be numerically hazardous even at inference time due to the so-called cancellation errors []. Later, we show that our HOPE parameterization does not suffer from these issues.

Combined with the spectral information of the state matrices $\mathbf{A}$, 2 lets us explain why systems initialized by $\texttt{init}_{1}$ are much more sensitive to perturbation than those initialized by $\texttt{init}_{3}$. In Figure 4, we show the locations of the poles $a_{j}$ and their associated magnitudes $|b_{j}c_{j}|$. Compared to the systems from $\texttt{init}_{3}$, the poles of the systems from $\texttt{init}_{1}$ are much closer to the imaginary axis whereas the values of $|b_{j}c_{j}|$ are also much larger. Hence, by 2, they are much more sensitive to training, and therefore lose numerical ranks easily. This is indeed seen in Figure 4.

2 delivers a disturbing message. As stated in [], “[the Hankel singular values] decay more rapidly the farther the $\Lambda(\mathbf{A})$ falls in the left half of the complex plane,” where $\Lambda(\mathbf{A})$ is the spectrum of $\mathbf{A}$, which is equivalent to $\{a_{1},\ldots,a_{n}\}$ in the diagonal case. Hence, many diagonal LTI systems with a high $\epsilon$-rank, i.e., those that we earlier found necessary for an SSM to capture the long-range dependency, would have eigenvalues $a_{j}$ close to the imaginary axis, making the system more sensitive to perturbation and thus the training less numerically stable.

Figure: Figure 4: A random perturbation to the imaginary part of $\mathbf{A}$ is added to a system from $\texttt{init}_{1}$ and a HiPPO-LegS system from $\texttt{init}_{3}$. The magnitude of the perturbation is set to $0.1\%$ and $1\%$ of the original matrix $\mathbf{A}$. For each system, on the left, we show the relative Hankel singular values $\sigma_{j}/\sigma_{1}$ of the original and perturbed systems; on the right, we plot the location of each $a_{j}$ in the complex plane and use the color to indicate the magnitude of its associated $|b_{j}c_{j}|$.

###### Theorem 2 .

## 4 HOPE-SSM: A Rankful, Stable, and Long-memory Parameterization

Given the potential issues of an SSM discussed in section 3, we propose an entirely different parameterization of the LTI systems called HOPE. We first describe the details of our HOPE-SSM and then explain how it resolves the low-rank and numerical instability issues of an LTI system. In addition, it also benefits from long-term memory. Our strategy is to use a Hankel matrix defined in eq. 4 to parameterize an LTI system. Instead of having $\mathbf{A}$, $\mathbf{B}$, $\mathbf{C}$, $\mathbf{D}$, and $\Delta t$ as the model parameters, we now have a vector $\mathbf{h}$ of length $n$, the skip connection $\mathbf{D}$, and $\Delta t$ as our model parameters. Hence, we use $n$ complex parameters from $\mathbf{h}$ to replace the $3n$ (resp. $4n$) complex parameters from $\mathbf{A}$, $\mathbf{B}$, and $\mathbf{C}$ in S4D (resp. S4); (^1^11There are different variants of S4D models that combine $\mathbf{B}$ and $\mathbf{C}$ into a single complex vector or use fewer copies of $\mathbf{A}$ than the number of channels. Here, we compare HOPE to a vanilla S4D model.) We remark that, as mentioned in the introduction, our model only modifies the LTI systems in an SSM. That is, it requires $1/3$ the number of parameters in an S4D model for each LTI system. Since there are other components in an SSM (e.g., the encoder and the decoder), we do not compress the number of parameters in the entire model by a factor of $1/3$. The finite Hankel matrix $\overline{\mathbf{H}}\in\mathbb{C}^{n\times n}$ is then defined by the Markov parameters in $\mathbf{h}$:

$$ $\overline{\mathbf{H}}_{i,j}=\mathbbm{1}_{\{i+j<n\}}\mathbf{h}_{i+j}.$ (5) $$

In an S4D model, we start with $\mathbf{A}$, $\mathbf{B}$, and $\mathbf{C}$ so that $\mathbf{h}_{i+j}=\overline{\mathbf{C}}\overline{\mathbf{A}}^{i+j}\overline{
\mathbf{B}}$; in a HOPE-SSM, we start with $\mathbf{h}_{i+j}$ and the matrix $\overline{\mathbf{H}}$ corresponds to a discrete LTI system $\overline{\Gamma}$, which is further associated with a continuous-time system $\Gamma$ via the bilinear transform with $\Delta t=1$. Notably, in our HOPE-SSM, it is the continuous system $\Gamma$ that takes the role of $(\mathbf{A},\mathbf{B},\mathbf{C})$ in a canonical SSM, e.g., S4 and S4D. In particular, $\Gamma$ is then discretized with a trainable sampling period $\Delta t$ for a discrete sequential input.

If we set $\Delta t=1$ and discretize $\Gamma$, the convolutional kernel is exactly $\mathbf{h}$ padded with $L-n$ zeros, where $L$ is the length of the sequential input. However, we usually want to discretize with a different $\Delta t$. We can do so via resampling the transfer functions of $\Gamma$ and $\overline{\Gamma}$, which equal []

$$ $G(s)=\sum_{j=0}^{n-1}\mathbf{h}_{j}\left((1+s)/(1-s)\right)^{-j-1}\qquad \Leftrightarrow\qquad\overline{G}(z)=\sum_{j=0}^{n-1}\mathbf{h}_{j}z^{-j-1}.$ (6) $$

Let $\bm{\omega}^{(L)}=\begin{bmatrix}\omega_{0}^{(L)}&\cdots&\omega_{L-1}^{(L)}
\end{bmatrix}^{\top}$ be the vector of the $L$th roots of unity. Given an input $\mathbf{u}\in\mathbb{C}^{L\times 1}$ of length $L$, when $\Delta t=1$, the outputs can be evaluated as

$$ $\mathbf{y}=\texttt{iFFT}(\texttt{FFT}(\mathbf{u})\circ\overline{G}(\bm{\omega} ^{(L)})),\qquad\overline{G}(\bm{\omega}^{(L)})=\begin{bmatrix}\overline{G}( \omega_{0}^{(L)})&\cdots&\overline{G}(\omega_{L-1}^{(L)})\end{bmatrix}^{\top}.$ $$

For a different $\Delta t$, one way to think of it is that we have compressed or dilated the time domain of $\mathbf{u}$ by a factor of $\Delta t$. Hence, the frequency domain of its Fourier transform $\hat{\mathbf{u}}$ is dilated or compressed by a factor of $1/\Delta t$. That is, we should relocate our samplers in the frequency domain as

$$ $\bm{\omega}^{(L,\Delta t)}=\begin{bmatrix}\omega_{0}^{(L,\Delta t)}&\cdots& \omega_{L-1}^{(L,\Delta t)}\end{bmatrix}^{\top},\qquad\omega_{j}^{(L,\Delta t) }=\frac{1+s^{(L)}_{j}/\Delta_{t}}{1-s^{(L)}_{j}/\Delta_{t}},\qquad s^{(L)}_{j} =\frac{\omega_{j}^{(L)}-1}{\omega_{j}^{(L)}+1},$ $$

where $s_{j}^{(L)}/\Delta t$ and $\omega_{j}^{(L,\Delta t)}$ are the scaled samplers in the time and angular domain, respectively. Then, the output sequence $\mathbf{y}$ can be computed, from the nonuniform samples $\overline{G}(\bm{\omega}^{(L,\Delta t)})$, as

$$ $\mathbf{y}=\texttt{iFFT}(\texttt{FFT}(\mathbf{u})\circ\overline{G}(\bm{\omega} ^{(L,\Delta t)})),\qquad\overline{G}(\bm{\omega}^{(L,\Delta t)})=\begin{ bmatrix}\overline{G}(\omega_{0}^{(L,\Delta t)})&\cdots&\overline{G}(\omega_{L- 1}^{(L,\Delta t)})\end{bmatrix}^{\top}.$ (7) $$

We provide the pseudocode in Algorithm 1 and a detailed derivation of eq. 7 in Appendix H.

Figure: Algorithm 1 Computing the output of an LTI system parameterized by its Hankel matrix.

With $L$ processors, we can sample $\overline{G}$ independently at $L$ locations in parallel, each of which takes $\mathcal{O}(n)$ time. Computing the FFT and iFFT takes $\mathcal{O}(L\log L)$. Overall, the evaluation of the Hankel-parameterized LTI system takes $\tilde{\mathcal{O}}(L+n)$ time and $\mathcal{O}(L)$ space, which agree with the complexities of the S4 and S4D models. To make the model recurrent during the inference time, one can either identify a system $(\mathbf{A},\mathbf{B},\mathbf{C})$ whose Hankel operator is $\mathbf{H}$ and compute as if it is an S4D model [, ], or directly compute the convolutional kernel using the iFFT of $\mathbf{g}$ in Algorithm 1.

Advantage I: A HOPE-SSM has a High Numerical Rank. The properties of a random Hankel matrix have been studied in [] and []. We can build upon their work to prove the following result.

###### Theorem 3 .

Let $\mathbf{h}_{1},\mathbf{h}_{2},\ldots$ be a sequence of i.i.d. random variables with mean 0 0 and variance $1$ that have finite third moments. We almost surely have that for any $\epsilon>0$, the $(\epsilon/\sqrt{\ln(n)})$-rank of $\overline{\mathbf{H}}_{n}$ is $\Omega(n)$, where $\overline{\mathbf{H}}_{n}$ is the $n\times n$ Hankel matrix defined in eq. 5.

Hence, if we ignore the logarithmic factor of $\sqrt{\ln(n)}$, then the numerical rank of a random Hankel matrix should be proportional to $n$ as $n\rightarrow\infty$. That means unlike the S4D model, we can always randomly initialize a high-rank LTI system with HOPE; even better, a system parameterized by HOPE does not lose rank during training (see Figure 5), because the low-rank systems are themselves rare in the space of the Markov parameters $\mathbf{h}$ (but not in the space of $(\mathbf{A},\mathbf{B},\mathbf{C})$).

Advantage II: A HOPE-SSM is Numerically Stable under Perturbation. Unlike the LTI system, a Hankel matrix $\overline{\mathbf{H}}$ is very numerically stable under perturbation, as shown in the following theorem.

###### Theorem 4 .

Let $\mathbf{h}\in\mathbb{C}^{n\times 1}$ be a vector and $G$ be the transfer function defined in eq. 6. Suppose we perturb $\mathbf{h}$ to $\tilde{\mathbf{h}}$ and let $\tilde{G}$ be its transfer function. Then, we have $\|G-\tilde{G}\|_{\infty}\leq\sqrt{n}\|\mathbf{h}-\tilde{\mathbf{h}}\|_{2}$.

Therefore, the numerical stability of the HOPE-SSM does not depend on the parameter $\mathbf{h}$ itself. Consequently, our HOPE-SSMs are trained without having to rescale the Markov parameters $\mathbf{h}$ (see []).

Advantage III: A HOPE-SSM’s Memory Does Not Fade. While many LTI systems are tailored for long-term memory [, ], an asymptotically stable system must inevitably suffer from an exponential decay in its memory []. This can be manifested by the Hankel matrix: since $\begin{bmatrix}\mathbf{y}_{0}^{\top}&\mathbf{y}_{1}^{\top}&\cdots\end{bmatrix}
^{\top}=\overline{\mathbf{H}}\begin{bmatrix}\mathbf{u}_{-1}^{\top}&\mathbf{u}_
{-2}^{\top}&\cdots\end{bmatrix}^{\top}$, we can write

$$ $\mathbf{y}_{j}=\sum_{k=1}^{\infty}\overline{\mathbf{H}}_{j,k-1}\mathbf{u}_{-k} =\sum_{k=1}^{\infty}\overline{\mathbf{H}}_{0,k+j-1}\mathbf{u}_{-k},\qquad j \geq 0.$ (8) $$

That means the effect of the input $\mathbf{u}_{-k}$ on the output $\mathbf{y}_{j}$ depends on the magnitude of $\overline{\mathbf{H}}_{0,k+j-1}$. Hence, the decay of $|\overline{\mathbf{H}}_{0,t}|$ gives us a measurement of the decay of the memory. For a discrete LTI system, we have by definition that $\overline{\mathbf{H}}_{0,t}=\overline{\mathbf{C}}\overline{\mathbf{A}}^{t}
\overline{\mathbf{B}}$. Since the spectrum of $\overline{\mathbf{A}}$ for an S4D model, regardless of $\Delta t$, is contained in the open unit disk, $\overline{\mathbf{A}}$ is a contraction operator and the “memory” of the system decays as time goes by. Asymptotically in time, this decay of $|\overline{\mathbf{H}}_{0,t}|$ is exponential, and for most systems, it starts at the very beginning when $t$ increases from 0 0.(^2^22In particular, this is the case for HiPPO-LegS (see Figure 6).) However, if we parameterize the LTI system by the Markov parameters in the Hankel matrix, then $|\overline{\mathbf{H}}_{0,t}|=\mathbf{h}_{t}$, which does not have to decay as long as $t<n$. We acknowledge that when $t\geq n$, we have $\overline{\mathbf{H}}_{0,t}=0$, which means the LTI system has essentially no memory after time $n$. We remark, however, that our HOPE-SSM utilizes the continuous-time LTI system associated with a Hankel matrix $\overline{\mathbf{H}}$. That is, from a continuous-time perspective, unlike S4D, our system’s memory does not decay for $t\in[0,n]$ (see Figure 6). Hence, even with a long sequence whose length is much larger than $n$, by setting the sampling period $\Delta t$ to be small enough, our HOPE-SSM could still enjoy non-decaying memory.

###### Theorem 3 .

###### Theorem 4 .

## 5 Experiments and Discussions

Figure: Figure 5: The test accuracy of the HOPE-SSM on the sCIFAR-10 task and the evolution of the Hankel singular values of the model. The plots are to be compared with Figure 2 and Figure 3.

Experiment I: Singular Values of a HOPE-SSM. In this section, we implement a randomly initialized HOPE-SSM to learn the sCIFAR-10 task. We use the same model hyperparameters as the S4D models in section 3. In particular, the Hankel matrices in this model are $64$-by-$64$. We randomly initialize the Hankel matrix and do not set a smaller learning rate for the Hankel matrix entries $\mathbf{h}$, i.e., all model parameters except for $\Delta t$ have the same learning rate. In Figure 5, we show the test accuracy and the Hankel singular values of the HOPE-SSM. Compared to Figure 2, a random HOPE-SSM can be trained to a high accuracy. In addition, training does not reduce the numerical rank of a Hankel matrix. This corroborates our findings in 3 and 4.

Figure: Figure 6: Left: the test accuracies of the S4D model and our HOPE-SSM on the noisy-sCIFAR task. Right: a unit impulse at $t=0$ is acted on the LTI system. The plot shows the decay of $|\mathbf{y}(t)|$ as $t$ increases. Data are collected for all $512$ LTI systems in a trained model. The dark curve shows the median of $|\mathbf{y}(t)|$ over the $512$ numbers. The darkly shaded region is from the first quartile to the third quartile. The lightly shaded region is from the minimum to the maximum.

Experiment II: HOPE-SSMs Have Long Memory. In section 4, we claim that a HOPE-SSM benefits from non-decaying memory. We show this experimentally in this section. To do so, for each flattened picture in the sCIFAR-10 dataset, which contains $1024$ vectors of length $3$, we append a random sequence of $1024$ vectors of length $3$ to the end of it. The goal is still to classify an image by its first $1024$ pixels. We call this task noise-padded sCIFAR-10. This task obviously requires a long memory so that the earlier data can still be retrieved after the noises are taken. We train both an S4D model and our HOPE-SSM to learn this task, using a common sampling period of $\Delta t=0.1$ (see Section I.2 for different values of $\Delta t$). We see in Figure 6 that the HOPE-SSM significantly outperforms the S4D model on this task. To understand this gap, we give a unit impulse at $t=0$ to the LTI systems in both the trained S4D model and our HOPE-SSM. We watch how the impulse response $\mathbf{y}(t)$ decays as time goes by. In Figure 6, we see that the memory of the trained S4D model decays exponentially while that of the trained HOPE-SSM does not decay at all for $t\in[0,64]$, where $n=64$ is the size of the Hankel matrix. This is exactly what we expected in eq. 8.

Experiment III: Performance in the Long-Range Arena. Finally, we present the performance of HOPE-SSM on large datasets. We use the same model architecture as in the S4D paper [], except that we replace the LTI blocks with our HOPE blocks. The specific model and training hyperparameters are reported in Appendix I. We show the performance of our model in Table 1, where we see that our HOPE-SSM outperforms most sequential models on many tasks.

**Table 1: Test accuracies in the Long-Range Arena of our HOPE-SSM and other models. We report the median and the standard deviation of executions with $5$ random seeds. The bold (resp. underlined) numbers indicate the best (resp. second best) performance on a task. An entry is left blank if no result is found. We use the same model sizes as those in the S4D model.**
| Model | ​​​ListOps​​​ | Text | ​​​Retrieval​​​ | Image | ​​​​Pathfinder​​​​ | Path-X | ​​​Avg. |
| --- | --- | --- | --- | --- | --- | --- | --- |
| DSS [] | $57.60$ | $76.60$ | $87.60$ | $85.80$ | $84.10$ | $85.00$ | $79.45$ |
| S4++ [] | $57.30$ | $86.28$ | $84.82$ | $82.91$ | $80.24$ | - | - |
| Reg. S4D [] | ${61.48}$ | $88.19$ | $91.25$ | $88.12$ | $94.93$ | $95.63$ | $86.60$ |
| Spectral SSM [] | $60.33$ | ${89.60}$ | $90.00$ | - | $\underline{95.60}$ | $90.10$ | - |
| Liquid S4 [] | $\bm{62.75}$ | $89.02$ | $91.20$ | $\bm{89.50}$ | $94.80$ | $96.66$ | $87.32$ |
| S5 [] | $62.15$ | $\underline{89.31}$ | $\underline{91.40}$ | $88.00$ | $95.33$ | $\bm{98.58}$ | $\underline{87.46}$ |
| S4 [] | $59.60$ | $86.82$ | $90.90$ | ${{88.65}}$ | $94.20$ | $96.35$ | $86.09$ |
| S4D [] | ${{60.47}}$ | $86.18$ | $89.46$ | $88.19$ | $93.06$ | $91.95$ | $84.89$ |
| HOPE-SSM | $\underline{{62.60}}$ | ${\bm{89.83}}$ | ${\bm{91.80}}$ | $\underline{{88.68}}$ | $\bm{{95.73}}$ | $\underline{{98.45}}$ | $\underline{\bm{87.85}}$ |
|  | $\qquad\pm 0.92$ | $\qquad\pm 0.37$ | $\qquad\pm 0.11$ | $\qquad\pm 0.44$ | $\qquad\pm 0.28$ | $\qquad\pm 0.17$ |  |

## 6 Conclusion

In this paper, we presented a new theory based on the Hankel singular values to understand the difficulties in initializing and training an SSM. We proposed a new parameterization scheme, called HOPE, that is based on the Markov parameters in a discrete Hankel operator. We proved that a HOPE-SSM can be robustly initialized and trained and has a long memory.

#### Acknowledgments

AY was supported by the SciAI Center, funded by the Office of Naval Research under Grant Number N00014-23-1-2729. NBE would like to acknowledge NSF, under Grant No. 2319621, and the U.S. Department of Energy, under Contract Number DE-AC02-05CH11231 and DE-AC02-05CH11231, for providing partial support of this work. We would also like to thank Alex Townsend, Anil Damle, and Sarah Dean for some inspiring discussions.

## Appendix

This is the collection of appendices for the paper titled “HOPE for a Robust Parameterization for Long-memory State Space Models.” It is organized as follows. In Appendix A, we survey the background of the bilinear transform of the LTI systems, the Mobius transform of the transfer function domains, and the Hankel operators more thoroughly. In Appendix B, we survey the basic properties of Hankel singular values and different strategies to compute them. In Appendix C, we explain the three different initialization schemes introduced in the paper. The proofs of 1, 2, and 3-4 are given in Appendix D, Appendix E, and Appendix F, respectively. In Appendix G, we conduct some numerical experiments to corroborate the four theorems presented in the main text. In Appendix H, we carefully derive how we can use the nonuniform samples of the transfer function to enable a different discretization step size $\Delta t$. Finally, the details of the experiments shown in the main text are given in Appendix I.

## Appendix A More background on the LTI system and Hankel operators

The purpose of this section is to provide a more thorough exposition of section 2, which leaves out some concepts that are unnecessary in order to understand the main text. Let $\Gamma=(\mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D})$ be a continuous-time LTI system defined in eq. 1. One can take a bilinear transformation to obtain a discrete LTI system:

$$ $\overline{\mathbf{A}}\!=\!(\mathbf{I}\!+\!\mathbf{A})(\mathbf{I}\!-\!\mathbf{A })^{-1},\quad\overline{\mathbf{B}}\!=\!(\mathbf{I}\!+\!\overline{\mathbf{A}}) \mathbf{B}/\sqrt{2},\quad\overline{\mathbf{C}}\!=\!\mathbf{C}(\mathbf{I}\!+\! \overline{\mathbf{A}})/\sqrt{2},\quad\overline{\mathbf{D}}\!=\!\mathbf{D}+ \overline{\mathbf{C}}(\mathbf{I}\!-\!\overline{\mathbf{A}})^{-1}\overline{ \mathbf{B}}.$ $$

The bilinear transformation is invertible and defines a one-to-one correspondence between $\Gamma$ in eq. 1 and $\overline{\Gamma}=(\overline{\mathbf{A}},\overline{\mathbf{B}},\overline{
\mathbf{C}},\overline{\mathbf{D}})$ given by

$$ $\displaystyle\overline{\mathbf{x}}_{k+1}$ $\displaystyle=\overline{\mathbf{A}}\overline{\mathbf{x}}_{k}+\overline{\mathbf {B}}\overline{\mathbf{u}}_{k},$ $\displaystyle\overline{\mathbf{y}}_{k}$ $\displaystyle=\overline{\mathbf{C}}\overline{\mathbf{x}}_{k}+\overline{\mathbf {D}}\overline{\mathbf{u}}_{k}.$ $$

The transfer functions of $\Gamma$ and $\overline{\Gamma}$ are

$$ $G(s)=\mathbf{C}(s\mathbf{I}-\mathbf{A})^{-1}\mathbf{B}+\mathbf{D}\qquad\text{ and}\qquad\overline{G}(z)=\overline{\mathbf{C}}(z{\mathbf{I}}-\overline{ \mathbf{A}})^{-1}\overline{\mathbf{B}}+\overline{\mathbf{D}},$ $$

respectively. The two transfer functions are equivalent via a Mobius transformation:

$$ $G(s)=\overline{G}\left((1+s)/(1-s)\right)\qquad\Leftrightarrow\qquad\overline{ G}(z)=G\left((z-1)/(z+1)\right).$ $$

The importance of the transfer functions is that they bring the inputs to the outputs in the Laplace domain by multiplication, which reduces to the Fourier domain if we assume that $\mathbf{u}$ is bounded and compactly supported:

$$ $\hat{\mathbf{y}}(s)=G(is)\hat{\mathbf{u}}(s),\qquad\hat{\overline{\mathbf{y}}} =\overline{G}({\bm{\omega}})\hat{\overline{\mathbf{u}}},$ (9) $$

where the hat of a function and a vector means the Fourier transform and the Fourier coefficients, respectively, and ${\bm{\omega}}$ is the vector of roots of unity.

Given a continuous-time LTI system $\Gamma$, one can define its Hankel operator by

$$ $\mathbf{H}:L^{2}(0,\infty)\rightarrow L^{2}(0,\infty),\qquad(\mathbf{H}\mathbf {v})(t)=\int_{0}^{\infty}\mathbf{C}\,\text{exp}((t+\tau)\mathbf{A})\mathbf{B} \mathbf{v}(\tau)d\tau.$ $$

The Hankel operator maps the past inputs to the future outputs, i.e., if $\mathbf{u}(t)=0$ for $t\geq 0$ and we define $\mathbf{v}(t)=\mathbf{u}(-t)$, then we have $\mathbf{y}(t)=(\mathbf{H}\mathbf{v})(t)$. Analogously, the Hankel matrix of a discrete LTI system $\overline{\Gamma}$ is a doubly infinite matrix defined by

$$ $\overline{\mathbf{H}}:\ell^{2}\rightarrow\ell^{2},\qquad\overline{\mathbf{H}}_ {i,j}=\overline{\mathbf{C}}\overline{\mathbf{A}}^{i+j}\overline{\mathbf{B}}, \qquad i,j\geq 0.$ $$

The Hankel matrix has the similar physical interpretation: if $\overline{\mathbf{u}}_{k}=0$ for all $k\geq 0$, then we have $\begin{bmatrix}\overline{\mathbf{y}}_{0}^{\top}&\overline{\mathbf{y}}_{1}^{
\top}&\cdots\end{bmatrix}^{\top}=\overline{\mathbf{H}}\begin{bmatrix}\overline
{\mathbf{u}}_{-1}^{\top}&\overline{\mathbf{u}}_{-2}^{\top}&\cdots\end{bmatrix}
^{\top}$. Both $\mathbf{H}$ and $\overline{\mathbf{H}}$ are bounded linear operators of rank $\leq n$, the number of latent states. In fact, the singular values $\sigma_{1}(\mathbf{H})\geq\sigma_{2}(\mathbf{H})\geq\cdots\geq\sigma_{n}(
\mathbf{H})\geq 0$ of $\mathbf{H}$ are equivalent to those of $\overline{\mathbf{H}}$, and they are called the Hankel singular values of $\Gamma$ and $\overline{\Gamma}$.

## Appendix B More background on Hankel singular values

In section 2, we define the Hankel singular values $\sigma_{1},\ldots,\sigma_{n}$ to be the singular values of the finite-rank bounded linear Hankel operator on a separable Hilbert space. Throughout the paper, we treated them as black boxes and used their distributions to understand the performance of SSMs. The goal of this section is to open the black boxes and introduce more useful properties of the Hankel singular values. Note that the proof of 1 to be presented in Appendix D heavily relies on the background discussed in this section. The concepts in this section can be presented on a continuous-time LTI system $(\mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D})$ or analogously on a discrete LTI system $(\overline{\mathbf{A}},\overline{\mathbf{B}},\overline{\mathbf{C}},\overline{
\mathbf{D}})$. We focus mainly on a continuous-time system for cleanliness. Since $\mathbf{D}$ has no effect on the Hankel singular values, we further assume that $\mathbf{D}=\bm{0}$ and write only $\Gamma=({\mathbf{A}},{\mathbf{B}},{\mathbf{C}})$.

### B.1 Hankel singular values from balanced realization

One of the most popular ways to compute the Hankel singular values is via the so-called balanced realization. Assume $({\mathbf{A}},{\mathbf{B}},{\mathbf{C}})$ is asymptotically stable, i.e. the spectrum of $\mathbf{A}$ is contained in the open left half-plane. Then, there exist two Hermitian and positive semi-definite matrices $\mathbf{P}\in\mathbb{C}^{n\times n}$ and $\mathbf{Q}\in\mathbb{C}^{n\times n}$ such that

$$ $\displaystyle\mathbf{A}\mathbf{P}+\mathbf{P}\mathbf{A}^{*}+\mathbf{B}\mathbf{B }^{*}$ $\displaystyle=\bm{0},$ (10) $\displaystyle\mathbf{A}^{*}\mathbf{Q}+\mathbf{Q}\mathbf{A}+\mathbf{C}^{*} \mathbf{C}$ $\displaystyle=\bm{0}.$ $$

The equations in eq. 10 are called the Lyapunov equations and $\mathbf{P}$ and $\mathbf{Q}$ are called the controllability Gramian and the observability Gramian, respectively. They can be explicitly expressed by the following matrix integrals:

$$ $\displaystyle\mathbf{P}=\int_{0}^{\infty}\text{exp}(\mathbf{A}t)\mathbf{B} \mathbf{B}^{*}\text{exp}(\mathbf{A}^{*}t)dt,\qquad\mathbf{Q}=\int_{0}^{\infty} \text{exp}(\mathbf{A}^{*}t)\mathbf{C}^{*}\mathbf{C}\text{exp}(\mathbf{A}t)dt.$ (11) $$

The controllability Gramian $\mathbf{P}$ is positive definite if and only if the system is controllable, i.e., for any $\mathbf{x}_{0},\mathbf{x}_{1}\in\mathbb{C}^{n}$ and any $T>0$, there exists an input $\mathbf{u}$ on $[0,T]$ that makes $\mathbf{x}(T)=\mathbf{x}_{1}$ when $\mathbf{x}(0)=\mathbf{x}_{0}$. Likewise, $\mathbf{Q}$ is positive definite if and only if the system is observable, i.e., for any $T>0$, the initial state $\mathbf{x}(0)$ can be determined by the input $\mathbf{u}$ and the output $\mathbf{y}$ on $[0,T]$ [].

The singular values of $\mathbf{P}\mathbf{Q}$ turn out to be exactly the squares of the Hankel singular values, i.e., $\sigma_{1}^{2},\ldots,\sigma_{n}^{2}$, of the system $\Gamma$. In general, $\mathbf{P}$ and $\mathbf{Q}$ are dense matrices. However, one can use the so-called balanced realization algorithm [] to compute an equivalent system $\Gamma_{b}=({\mathbf{A}}_{b},{\mathbf{B}}_{b},{\mathbf{C}}_{b})=(\mathbf{V}^{-
1}{\mathbf{A}}\mathbf{V},\mathbf{V}^{-1}{\mathbf{B}},{\mathbf{C}}\mathbf{V})$, where $\mathbf{V}\in\mathbb{C}^{n\times n}$ is an invertible matrix, so that its Gramians $\mathbf{P}_{b}$ and $\mathbf{Q}_{b}$ are equivalent and diagonal, i.e.,

$$ $\mathbf{P}_{b}=\mathbf{Q}_{b}=\text{diag}(\sigma_{1},\ldots,\sigma_{n}).$ $$

Balanced realization is an algebraic method that is based on singular value decompositions (SVDs). In practice, it is a very popular method to compute the Hankel singular values of a system.

### B.2 Hankel singular values as a rational approximation problem

The balanced realization gives us a good way to compute the Hankel singular values in practice. However, it does not offer too much insight in theoretically analyzing them. The theory of Hankel singular values is usually derived via function approximation. Here, we introduce how a Hankel singular value can be reframed as the solution to a rational approximation problem. To this end, we let $H^{\infty}_{+}$ (resp. $H^{\infty}_{-}$) be the Hardy space with functions $h:\mathbb{C}\rightarrow\mathbb{C}^{m\times p}$ that are bounded and analytic in the open right (resp. left) half-plane. The non-tangential limit of $h$ exists almost everywhere on the imaginary axis, and by the Maximum Modulus Principle, we must have

$$ $\|h\|_{\infty}:=\text{ess\,sup}_{\text{Re}(s)=0}\|h(s)\|_{2}=\sup_{\text{Re}(s )>0}\|h(s)\|_{2},\qquad h\in H^{\infty}_{+}$ $$

and

$$ $\|h\|_{\infty}:=\text{ess\,sup}_{\text{Re}(s)=0}\|h(s)\|_{2}=\sup_{\text{Re}(s )<0}\|h(s)\|_{2},\qquad h\in H^{\infty}_{-}.$ $$

The Adamyan–Arov–Krein (AAK) Theory [] says the following:

$$ $\sigma_{k+1}=\inf_{R_{k},F}\|G-R_{k}-F\|_{\infty},\qquad 0\leq k\leq n-1,$ (12) $$

where $R_{k}$ ranges over all rational functions with at most $k$ poles, all contained in the open left half-plane, and $F$ ranges over $H^{\infty}_{-}$. In fact, as discussed in section 2, $R_{k}$ is the transfer function of a reduced system $\tilde{\Gamma}=(\tilde{\mathbf{A}},\tilde{\mathbf{B}},\tilde{\mathbf{C}},
\tilde{\mathbf{D}})$ with $\tilde{\mathbf{A}}\in\mathbb{C}^{k\times k}$. Hence, eq. 12 tells us that if the transfer function of $\Gamma$ can be well-approximated by rational functions, then it has fast-decaying Hankel singular values. The other direction is not true, due to the existence of $F$. That is, if the Hankel singular values decay fast, then it does not necessarily mean that $G$ can be well-approximated by rational functions. However, section 2 shows that this is true if the sum of the tails of the Hankel singular values decay rapidly.

## Appendix C Three Different Initialization Schemes

In this section, we explain the details of $\texttt{init}_{1}$, $\texttt{init}_{2}$, and $\texttt{init}_{3}$ in the main text. We stress again that we do not claim that these initialization schemes are representative. Indeed, we only use them to elicit the story of the Hankel singular values. As mentioned in the main text, $\texttt{init}_{3}$ is the HiPPO-LegS initialization scheme []. To sample a random with $\texttt{init}_{1}$, we create random samples of the transfer function $\{(is_{j},G_{j})\}_{j=1}^{N}$, where $s_{j}\in\mathbb{R}$ and $G_{j}\in\mathbb{C}$. We then identify a system $\Gamma$ whose transfer function $G$ satisfies that $G(is_{j})\approx G_{j}$. The way we identify this model is via the so-called AAA algorithm [, ]. For $\texttt{init}_{2}$, we sample a discrete system by assuming the diagonal entries of $\overline{\mathbf{A}}=\text{diag}(a_{1},\ldots,a_{n})$ are uniformly sampled on the unit disk and $\overline{\mathbf{B}}\circ\overline{\mathbf{C}}^{\top}$ is a random vector with each entry sampled i.i.d. from a normal distribution $\mathcal{N}(0,1)$. We then compute the corresponding continuous-time system from the discrete system using the bilinear transform.

## Appendix D Proof of 1

Let $\overline{G}$ be the transfer function of a random system $\overline{\Gamma}_{2}$, i.e.,

$$ $\overline{G}(z)=\overline{\mathbf{C}}(z\mathbf{I}-\overline{\mathbf{A}})^{-1} \overline{\mathbf{B}}+\overline{\mathbf{D}}.$ $$

By the AAK theory, the Hankel singular values of $\overline{\Gamma}$ can be studied via the rational approximation of $\overline{G}$. Since the matrix $\overline{\mathbf{D}}$ does not affect the Hankel singular values of a system, we assume, without loss of generality, that $\overline{\mathbf{D}}=\bm{0}$. We let $\sigma_{1},\ldots,\sigma_{n}$ be the random variables that are equal to the singular values of the random LTI system $\overline{\Gamma}_{2}$. We approach 1 in three steps:

- 1.
We separate out the poles of the transfer function $G$ of $(\overline{\mathbf{A}},\overline{\mathbf{B}},\overline{\mathbf{C}},\bm{0})$ that are close to the boundary of the unit disk. This breaks $\overline{G}$ into two low-rank systems $\overline{G}_{1}$ and $\overline{G}_{2}$, where $\overline{G}_{1}$ is low-rank because it has few poles and $\overline{G}_{2}$ is low-rank because its poles are far away from the boundary.
- 2.
We then estimate the decay of the singular values of $\overline{G}_{2}$ using the information of the maximum moduli of its poles.
- 3.
Finally, we control $\sigma_{1}$ in probability. This gives a control on the relative singular values.

We will make these three steps into three lemmas and use them to derive the result at the end.

###### Lemma 1 .

Given $\gamma>0$, with probability at least $1-\delta$, there are at most $n^{\beta}$ poles of $\overline{G}(z)$ outside the disk $D(0,1-n^{-\gamma})$, where

$$ $\beta=1+\log_{n}\left(n^{-\gamma\alpha}+\sqrt{\frac{\ln(1/\delta)}{2n}}\right).$ $$

###### Proof.

Let $Z$ be the number of poles inside $D(0,1-n^{-\gamma})$. Then, $Z$ has a binomial distribution

$$ $Z\sim B(n,1-n^{-\gamma\alpha}).$ $$

From Hoeffding’s inequality, we have

$$ $\displaystyle\mathbb{P}(Z\leq n-n^{\beta})\leq\text{exp}\left(-2n\left(1-n^{- \gamma\alpha}-\frac{n-n^{\beta}}{n}\right)^{2}\right)=\text{exp}\left(-2n\left (-n^{-\gamma\alpha}+n^{\beta-1}\right)^{2}\right)$ $$

Set

$$ $\beta=1+\log_{n}\left(n^{-\gamma\alpha}+\sqrt{\frac{\ln(1/\delta)}{2n}}\right) =1-\gamma\alpha+\log_{n}\left(1+\sqrt{\frac{\ln(1/\delta)}{2n^{1-2\gamma\alpha }}}\right).$ $$

Then, we have

$$ $n^{\beta-1}=n^{-\gamma\alpha}+\sqrt{\frac{\ln(1/\delta)}{2n}}$ $$

so that

$$ $\mathbb{P}(Z\leq n-n^{\beta})\leq\text{exp}\left(-2n\frac{\ln(1/\delta)}{2n} \right)=\delta.$ $$

This finishes the proof.
∎

###### Lemma 2 .

For any $\gamma>0$, let the random variable $k$ be the number of poles of $\overline{G}(z)$ inside $D(0,1-n^{-\gamma})$. Let $\kappa>\gamma$ be given. Then, with conditional probability (given $k$) at least $1-\delta$, we have

$$ $\sigma_{(n-k)+n^{\kappa}+2}\leq\mathcal{O}\left(\!\sqrt{n^{\gamma+1}(\ln(1/ \delta)\!+\!n)}\right)e^{-(n^{(\kappa-\gamma)})},$ $$

where the constant in $\mathcal{O}$ is universal.

###### Proof.

Let $z_{1},\ldots,z_{k}$ be the poles of $\overline{G}(z)$ inside $D(0,1-n^{-\gamma})$. Assume $G(z)$ can be written as

$$ $\overline{G}(z)=\overline{G}_{1}(z)+\sum_{i=1}^{k}\frac{c_{i}}{z-z_{i}},$ $$

where $\overline{G}_{1}(z)$ is a degree-$(n-k)$ rational function with poles inside the annulus of inner radius $1-n^{-\gamma}$ and outer radius $1$ and $c_{i}$’s are i.i.d. random variables with distribution $\mathcal{N}(0,1)$. We can further write

$$ $\sum_{i=1}^{k}\frac{c_{i}}{z-z_{i}}=\sum_{i=1}^{k}c_{i}\sum_{j=0}^{\infty}z_{i }^{j}z^{-j-1}=\sum_{j=1}^{\infty}z^{-1-j}\left(\sum_{i=1}^{k}c_{i}z_{i}^{j} \right).$ $$

By the AAK theory, for any $K>0$, we have

$$ $\sigma_{(n-k)+K}\leq\sup_{\left|z\right|=1}\left|\sum_{j=K-2}^{\infty}z^{-1-j} \left(\sum_{i=1}^{k}c_{i}z_{i}^{j}\right)\right|\leq\sum_{j=K-2}^{\infty}\| \mathbf{c}\;\mathbf{z}^{j}\|_{1}\leq\|\mathbf{c}\|_{2}\sum_{j=K-2}^{\infty}\| \mathbf{z}^{j}\|_{2},$ $$

where $\mathbf{c}=[c_{1},\ldots,c_{k}]^{\top}$ and $\mathbf{z}^{j}=[z_{1}^{j},\ldots,z_{k}^{j}]^{\top}$, and the last step follows from Hölder’s inequality. Since $\|\mathbf{c}\|_{2}^{2}$ follows the $\chi^{2}_{k}$ distribution, for $M>n$, we have that

$$ $P(\|\mathbf{c}\|_{2}^{2}>M)\leq\text{exp}\left(-\frac{n}{2}\left(\frac{M}{n}-1 -\ln\left(\frac{M}{n}\right)\right)\right).$ $$

Hence, there exists a universal constant $C>0$ such that when $M\geq C(\ln(1/\delta)+n)$(^3^33Here, $n$ is to guarantee that $M>n$), we have

$$ $P(\|\mathbf{c}\|_{2}^{2}>M)\leq\text{exp}\left(-\frac{n}{2}\frac{M}{n}\right) \leq\delta.$ $$

Moreover, since $\left|z_{i}^{j}\right|\leq(1-n^{-\gamma})^{j}$ for all $1\leq i\leq k$, we have

$$ $\|\mathbf{z}^{j}\|_{2}\leq\sqrt{n}(1-n^{-\gamma})^{j}.$ $$

That is, with probability at least $1-\delta$, we have

$$ $\sigma_{(n-k)+K}\leq\mathcal{O}\!\left(\!\sqrt{\ln(1/\delta)\!+\!n}\right)\! \sqrt{n}\!\!\sum_{j=K-2}^{\infty}(1-n^{-\gamma})^{j}=\mathcal{O}\left(\!\sqrt{ n^{\gamma+1}(\ln(1/\delta)\!+\!n)}\right)(1\!-\!n^{-\gamma})^{K-2}.$ $$

Suppose $K\geq n^{\kappa}+2$. Then, we have

$$ $(1\!-\!n^{-\gamma})^{K-2}\leq(1\!-\!n^{-\gamma})^{(n^{\kappa})}=\left((1\!-\!n ^{-\gamma})^{(n^{\gamma})}\right)^{(n^{\kappa-\gamma})}.$ $$

Since $(1\!-\!n^{-\gamma})^{n^{\gamma}}\rightarrow e^{-1}$ as $n\rightarrow\infty$, if $\kappa>\gamma$, then $(1\!-\!n^{-\gamma})^{K-2}$ decays faster than any negative power of $n$. Hence, we have

$$ $\sigma_{(n-k)+K}\leq\mathcal{O}\left(\!\sqrt{n^{\gamma+1}(\ln(1/\delta)\!+\!n) }\right)e^{-(n^{(\kappa-\gamma)})},$ $$

as desired.
∎

###### Lemma 3 .

With probability at least $1-\delta$, the leading Hankel singular value $\sigma_{1}$ satisfies

$$ $\sigma_{1}\geq\mathcal{O}(\sqrt{n}\delta),$ $$

where the constant in $\mathcal{O}$ is universal.

###### Proof.

The leading Hankel singular value $\sigma_{1}$ is equivalent to the spectral norm of the Hankel matrix

$$ $\begin{bmatrix}\overline{\mathbf{C}}\overline{\mathbf{B}}&\overline{\mathbf{C} }\overline{\mathbf{A}}\overline{\mathbf{B}}&\overline{\mathbf{C}}\overline{ \mathbf{A}}^{2}\overline{\mathbf{B}}&\cdots\\ \overline{\mathbf{C}}\overline{\mathbf{A}}\overline{\mathbf{B}}&\overline{ \mathbf{C}}\overline{\mathbf{A}}^{2}\overline{\mathbf{B}}&\cdots&\cdots\\ \overline{\mathbf{C}}\overline{\mathbf{A}}^{2}\overline{\mathbf{B}}&\vdots& \ddots&\vdots\\ \vdots&\vdots&\cdots&\ddots\end{bmatrix}.$ $$

Hence, we have

$$ $\sigma_{1}\geq\left|\overline{\mathbf{C}}\overline{\mathbf{B}}\right|,$ $$

where $\overline{\mathbf{C}}\overline{\mathbf{B}}/\sqrt{n}\sim\mathcal{N}(0,1)$. Hence, if we set $M=C\sqrt{n}\delta$ for some universal constant $C>0$, we have

$$ $P(\sigma_{1}\leq M)\leq P\big{(}\left|\overline{\mathbf{C}}\overline{\mathbf{B }}\right|\leq M\big{)}=P\big{(}\left|\overline{\mathbf{C}}\overline{\mathbf{B} }\right|/\sqrt{n}\leq M/\sqrt{n}\big{)}\leq\delta.$ $$

This completes the proof.
∎

Now, let’s assemble our ultimate statement.

###### Proof of 1 .

By 2 and 3, with probability at least $1-\delta/2$, we have that

$$ $\frac{\sigma_{(n-k+n^{\kappa}+2)}}{\sigma_{1}}\leq\mathcal{O}\left(\!\sqrt{n^{ \gamma}(\ln(1/\delta)\!+\!n)}\delta^{-1}\right)e^{-(n^{(\kappa-\gamma)})}.$ $$

Hence, if we set $\kappa$ so that

$$ $n^{\kappa}\geq n^{\gamma}\ln\left(C\delta^{-1}\epsilon^{-1}\sqrt{n^{\gamma}( \log(1/\delta)+n)}\right)=n^{\gamma}\ln\left(\mathcal{O}\left(\delta^{-3/2} \epsilon^{-1}n^{(\gamma+1)/2}\right)\right)$ $$

for a sufficiently large universal constant $C>0$, then we guarantee that

$$ $\frac{\sigma_{(n-k+n^{\kappa}+2)}}{\sigma_{1}}\leq\epsilon.$ $$

By 1, we have that with probability at least $1-\delta/2$,

$$ $n-k\leq n^{\beta},\qquad\beta=1+\log_{n}\left(n^{-\gamma\alpha}+\sqrt{\frac{ \ln(2/\delta)}{2n}}\right).$ $$

Set $\gamma=1/(1+\alpha)$. Then, we have

$$ $\displaystyle\beta$ $\displaystyle=1+\log_{n}\left(n^{-\gamma\alpha}\left(1+\sqrt{\frac{\ln(2/ \delta)}{2n^{1-2\gamma\alpha}}}\right)\right)=1-\gamma\alpha+\log_{n}\left(1+ \sqrt{\frac{\ln(2/\delta)}{2n^{1-2\gamma\alpha}}}\right)$ $\displaystyle\leq\frac{1}{1+\alpha}+\log_{n}(1+\sqrt{\ln(2/\delta)/2})<\frac{1 }{1+\alpha}+\frac{\ln(2+\sqrt{\ln(1/\delta)/2})}{\ln(n)}.$ $$

Since

$$ $n^{\kappa}+2=\mathcal{O}\left(n^{\gamma}\ln\left(\delta^{-3/2}\epsilon^{-1}n \right)\right)\leq\mathcal{O}\left(n^{\beta}\ln\left(\delta^{-3/2}\epsilon^{-1 }n\right)\right).$ $$

The claim is proved.
∎

###### Lemma 1 .

###### Proof.

###### Lemma 2 .

###### Proof.

###### Lemma 3 .

###### Proof.

###### Proof of 1 .

## Appendix E Proof of 2

In this section, we prove 2. Our proof focuses on the worst-case perturbation by construction. As a consequence, it simultaneously proves the sharpness of the result. Intuitively, consider a rational function

$$ $s\mapsto\frac{bc}{s-a},$ $$

since we only care about its values on the imaginary axis, the closer the pole $a$ is to the imaginary axis, the less stable it is. On the other hand, it is obvious that the size of $bc$ also controls the (absolute) conditioning of the rational function. We state the rigorous proof below.

###### Proof of 2 .

Without loss of generality, we assume that $\mathbf{B}=\tilde{\mathbf{B}}=\begin{bmatrix}1&1&\cdots&1\end{bmatrix}^{\top}$. (^4^44Otherwise, we can redefine $\mathbf{B}=\tilde{\mathbf{B}}=\begin{bmatrix}1&1&\cdots&1\end{bmatrix}^{\top}$, $\mathbf{C}=\mathbf{B}^{\top}\circ\mathbf{C}$, and $\tilde{\mathbf{C}}=\tilde{\mathbf{B}}^{\top}\circ\tilde{\mathbf{C}}$, and the redefined transfer functions are unchanged.) The transfer functions of $\Gamma$ and $\tilde{\Gamma}$ are

$$ $G(s)=\sum_{j=1}^{n}\frac{c_{j}}{s-s_{j}}\qquad\text{and}\qquad\tilde{G}(s)= \sum_{j=1}^{n}\frac{\tilde{c}_{j}}{s-\tilde{s}_{j}},$ $$

respectively. Then, for any $s$ on the imaginary axis, we have

$$ $\displaystyle\left|\frac{c_{j}}{s-s_{j}}-\frac{\tilde{c}_{j}}{s-\tilde{s}_{j}}\right|$ $\displaystyle=\left|\frac{c_{j}(s-\tilde{s}_{j})-\tilde{c}_{j}(s-s_{j})}{(s-s_ {j})(s-\tilde{s}_{j})}\right|=\left|\frac{c_{j}s-c_{j}\tilde{s}_{j}-\tilde{c}_ {j}s+\tilde{c}_{j}s_{j}+\tilde{c}_{j}\tilde{s}_{j}-\tilde{c}_{j}\tilde{s}_{j}} {(s-s_{j})(s-\tilde{s}_{j})}\right|$ $\displaystyle\leq\frac{\left|c_{j}-\tilde{c}_{j}\right|\left|s-\tilde{s}_{j} \right|+\left|\tilde{c}_{j}\right|\left|s_{j}-\tilde{s}_{j}\right|}{\left|s-s_ {j}\right|\left|s-\tilde{s}_{j}\right|}=\frac{\left|c_{j}-\tilde{c}_{j}\right| }{\left|s-s_{j}\right|}+\frac{\left|\tilde{c}_{j}\right|\left|s_{j}-\tilde{s}_ {j}\right|}{\left|s-{s}_{j}\right|\left|s-\tilde{s}_{j}\right|}$ $\displaystyle\leq\frac{\Delta_{\mathbf{B}}}{|\text{Re}(s_{j})|}+\frac{2\left|c _{j}\right|\Delta_{\mathbf{A}}}{|\text{Re}(s_{j})|^{2}/2}.$ $$

Hence, we have

$$ $\displaystyle\|G-\tilde{G}\|_{\infty}$ $\displaystyle\leq\sum_{j=1}^{n}\left(\frac{\Delta_{\mathbf{B}}}{|\text{Re}(s_{ j})|}+4\frac{\left|c_{j}\right|\Delta_{\mathbf{A}}}{|\text{Re}(s_{j})|^{2}}\right)$ $\displaystyle\leq n\Delta_{\mathbf{B}}\max_{j}\frac{1}{|\text{Re}(s_{j})|}+4n \Delta_{\mathbf{A}}\max_{j}\frac{\left|c_{j}\right|}{\left|\text{Re}(s_{j}) \right|^{2}}.$ $$

This proves the upper bound. To prove the lower bound, let $j_{1}$ be an index that maximizes $1/|\text{Re}(s_{j})|$ and $j_{2}$ be an index that maximizes $|c_{j}|/|\text{Re}(s_{j})|^{2}$. Define $\tilde{\Gamma}_{\mathbf{B}}$ by perturbing $c_{j_{1}}$ to $c_{j_{1}}+\Delta_{\mathbf{B}}$. Then, we have

$$ $\|G-\tilde{G}_{\mathbf{B}}\|_{\infty}=\left\|\frac{\Delta_{\mathbf{B}}}{s-s_{j _{1}}}\right\|_{\infty}=\Delta_{\mathbf{B}}\max_{j}\frac{1}{|\text{Re}(s_{j})|}.$ $$

Define $\tilde{\Gamma}_{\mathbf{A}}$ by perturbing $s_{j_{2}}$ to $s_{j_{2}}+\Delta_{\mathbf{A}}$. Then, we have

$$ $\displaystyle\|G-\tilde{G}_{\mathbf{A}}\|_{\infty}$ $\displaystyle=\left|c_{j_{2}}\right|\left\|\frac{1}{s-s_{j_{2}}}-\frac{1}{s-s_ {j_{2}}+\Delta_{\mathbf{A}}}\right\|_{\infty}=\left|c_{j_{2}}\right|\Delta_{ \mathbf{A}}\left\|\frac{1}{(s-s_{j_{2}})(s-s_{j_{2}}+\Delta_{\mathbf{A}})} \right\|_{\infty}$ $\displaystyle\geq\left|c_{j_{2}}\right|\Delta_{\mathbf{A}}\frac{1}{\left|\text {Re}(s_{j_{2}})\right|^{2}}=\Delta_{\mathbf{A}}\max_{j}\frac{\left|c_{j}\right |}{\left|\text{Re}(s_{j})\right|^{2}}.$ $$

This proves the sharpness of the theorem.
∎

###### Proof of 2 .

## Appendix F Proof of 3 and 4

The proof of 3 is a straightforward assembly of two results in random matrix theory. The first result, due to [], controls the Hankel norm $\sigma_{1}(\overline{\mathbf{H}}_{n})$ of a random Hankel matrix, whereas the second result by [] studies the distribution of all absolute singular values $\sigma_{j}(\overline{\mathbf{H}}_{n})$ of a random Hankel matrix. Our study of the relative Hankel singular values is achieved by taking the quotient of the subjects of the two prior works.

###### Proof of 3 .

By [], with probability $1$, we have that

$$ $\sigma_{1}(\overline{\mathbf{H}}_{n})=\|\overline{\mathbf{H}}_{n}\|=\mathcal{O }(\sqrt{n\ln n}).$ $$

Define $\overline{\mathbf{K}}_{n}=\overline{\mathbf{H}}_{n}(1\!:\!\lceil n/2\rceil,1\!
:\!\lceil n/2\rceil)$. Then, by [], with probability $1$, we have that $\mu(\overline{\mathbf{K}}_{n}/\sqrt{n})$ converges in distribution to a fixed probability measure, where

$$ $\mu(\overline{\mathbf{K}}_{n}/\sqrt{n})=\frac{1}{\lceil n/2\rceil}\sum_{j=1}^{ \lceil n/2\rceil}\delta_{\lambda_{j}(K_{n}/\sqrt{n})}$ $$

is the spectral measure of $\overline{\mathbf{K}}_{n}/\sqrt{n}$. Since $\overline{\mathbf{K}}_{n}$ is symmetric, the singular values of $\overline{\mathbf{K}}_{n}$ are the moduli of the eigenvalues of $\overline{\mathbf{K}}_{n}$. Hence, fix some $\epsilon>0$, we have that

$$ $\left|\{j\mid\sigma_{j}(\overline{\mathbf{K}}_{n})/\sigma_{1}(\overline{ \mathbf{H}}_{n})>\epsilon/\sqrt{\ln(n)}\}\right|=\Omega(n).$ $$

Since $\overline{\mathbf{K}}_{n}$ is a submatrix of $\overline{\mathbf{H}}_{n}$, we have $\sigma_{j}(\overline{\mathbf{H}}_{n})\geq\sigma_{j}(\overline{\mathbf{K}}_{n})$ for all $1\leq j\leq\lceil n/2\rceil$. Hence, its $(\epsilon/\sqrt{\ln(n)})$-rank can be controlled as

$$ $\left|\{j\mid\sigma_{j}(\overline{\mathbf{H}}_{n})/\sigma_{1}(\overline{ \mathbf{H}}_{n})>\epsilon/\sqrt{\ln(n)}\}\right|\geq\left|\{j\mid\sigma_{j}( \overline{\mathbf{K}}_{n})/\sigma_{1}(\overline{\mathbf{H}}_{n})>\epsilon/ \sqrt{\ln(n)}\}\right|=\Omega(n).$ $$

This finishes the proof.
∎

We remark that in the statement of 3, we study the $\epsilon/\sqrt{\ln(n)}$-rank of the Hankel matrix instead of the $\epsilon$-rank. The reason is that, unlike the singular values of a random matrix, the spectral measure of a normalized random Hankel matrix $\mathbf{H}_{n}/\sqrt{n}$ has unbounded support. In order to pull the $1/\sqrt{\ln(n)}$ factor out and make it into the $\Omega(n)$ bound, we need to study the distribution of the spectral measure of $\mathbf{H}_{n}/\sqrt{n}$. As pointed out by [], however, this seems to be a hard problem. Nevertheless, we can empirically test the statement by numerical experiments.

Next, we provide the proof of 4, which is almost immediate from Hölder’s inequality.

###### Proof of 4 .

By Hölder’s inequality, we have

$$ $\|G-\tilde{G}\|_{\infty}=\|\overline{G}-\overline{\tilde{G}}\|_{\infty}\leq \max_{|z|=1}\sum_{j=0}^{n-1}|\mathbf{h}_{j}-\tilde{\mathbf{h}}_{j}||z|^{-j-1} \leq\|\mathbf{h}-\tilde{\mathbf{h}}\|_{2}\sqrt{n}.\qed$ $$

###### Proof of 3 .

###### Proof of 4 .

## Appendix G Some numerical experiments on Hankel matrices

### G.1 Numerical ranks of random LTI systems and random Hankel matrices in practice

While the theoretical part of our paper focuses on the $\epsilon$-rank of an LTI system, in the main text, we showed the distribution of all Hankel singular values of different LTI systems. The main reason for showing the histograms instead of a single number (i.e., the $\epsilon$-rank) is that the histogram gives us more information while the $\epsilon$-rank is merely a cutoff. In this section, we empirically compute the $\epsilon$-rank to verify the two theorems (i.e., 1 and 3).

In this experiment, we always set $\epsilon=0.01$. For every $n$ in our experiment, we first randomly initialize a random $n\times n$ Hankel matrix

$$ $\overline{\mathbf{H}}_{n}=\begin{bmatrix}{h}_{0}&{h}_{1}&{h}_{2}&\cdots&{h}_{n -1}\\ {h}_{1}&{h}_{2}&\cdots&{h}_{n-1}&0\\ {h}_{2}&\cdots&{h}_{n-1}&0&0\\ \vdots&\iddots&\iddots&\vdots&\vdots\\ {h}_{n-1}&0&\cdots&0&0\end{bmatrix},$ $$

where $h_{j}$ are i.i.d. random Gaussian variables with variance of $1$. We compute its $\epsilon$-rank and we repeat the experiment for $1000$ trials. Similarly, for every $n$, we randomly initialize an LTI system with

$$ $\mathbf{A}=\text{diag}(a_{1},\ldots,a_{n}),\qquad a_{j}\sim\text{Uniform}( \mathbb{D}),$ $$

where $\mathbb{D}$ is the open unit disk in the complex plane and the elements of $\mathbf{B}\circ\mathbf{C}^{\top}$ are sampled i.i.d. from $\mathcal{N}(0,1)$. We compute its $\epsilon$-rank and also repeat the experiment for $1000$ trials. From Figure 7, we see that a random LTI system has a low rank, whereas a random Hankel matrix has a high rank in the sense that it is about proportional to $n$. This observation aligns with our theory in 1 and 3,

Figure: Figure 7: The $\epsilon$-rank of a random LTI system of size $n$, where $\epsilon=0.01$. The random systems are parameterized by a Hankel matrix or by the matrices $\mathbf{A},\mathbf{B}$, and $\mathbf{C}$. The lines are the average rank and the shaded regions indicate the $10\%$-$90\%$ range over $1000$ trials.

### G.2 Hankel matrices are stable to perturbation in practice

4 predicts that the Hankel matrices are very stable when being perturbed. In this section, we run experiments in parallel with those in Figure 4, where we perturb a Hankel matrix $\overline{\mathbf{H}}$. As in Figure 4, we also set the size of the random perturbation to be $1\%$ and $0.1\%$, respectively, of the original system. 4 is corroborated by Figure 8, where we see that a small perturbation has a minimal effect on the Hankel singular values.

Figure: Figure 8: A random perturbation is added to $\overline{\mathbf{H}}$. The magnitude of the perturbation is set to $0.1\%$ and $1\%$ of the original matrix. We show the relative Hankel singular values $\sigma_{j}/\sigma_{1}$ of the original and perturbed systems. In this case, the three curves are almost overlapping.

## Appendix H Adjusting the discretization step via NUFFT

Assume we have a discrete LTI system whose transfer function is $\overline{G}$. As explained in section 4, if we assume $\Delta t=1$, then the output of the system given an input sequence can be computed as

$$ $\mathbf{y}=\texttt{iFFT}(\texttt{FFT}(\mathbf{u})\circ\overline{G}(\bm{\omega} ^{(L)})),\qquad\overline{G}(\bm{\omega}^{(L)})=\begin{bmatrix}\overline{G}( \omega_{0}^{(L)})&\cdots&\overline{G}(\omega_{L-1}^{(L)})\end{bmatrix}^{\top}.$ $$

To understand why this relationship holds, assume there exists a continuous function $u$ on the unit circle $\partial\mathbb{D}$ where the discrete inputs $\mathbf{u}$ are sampled from. Then, FFT allows us to write $u$ into the Fourier expansion:

$$ $u(z)=\sum_{j=0}^{L-1}[\texttt{FFT}(\mathbf{u})]_{j}\;\text{exp}\left(-2\pi i \frac{j}{L}z\right).$ $$

By the property of the transfer function eq. 3, we know that the output function $y$ is equal to

$$ $y(z)=\sum_{j=0}^{L-1}\underbrace{\left([\texttt{FFT}(\mathbf{u})]_{j}\; \overline{G}(\omega_{j}^{(L)})\right)}_{\hat{\mathbf{y}}_{j}}\text{exp}\left(- 2\pi i\frac{j}{L}z\right).$ $$

To compute the discrete output $\mathbf{y}$, one samples $y$ at $z=\omega_{0}^{(L)},\ldots,\omega_{L-1}^{(L)}$, which is equivalent to an inverse FFT on $\texttt{FFT}(\mathbf{u})\circ\overline{G}(\bm{\omega}^{(L)})$.

Now, if we want to change $\Delta t$, one way to think of it is as if our LTI system is unchanged, but the time domain of $u(z)$ is scaled by a factor of $\Delta t$. That is, we now have(^5^55Note that we could alternatively scale the angular domain instead of the time domain, i.e., $z^{(\Delta t)}=z/\Delta t$. The difference is on the level of discretization. However, we find that discretizing the time domain gives us a better performance in general.)

$$ $u^{(\Delta t)}(z)=\sum_{j=0}^{L-1}[\texttt{FFT}(\mathbf{u})]_{j}\;\text{exp} \left(-2\pi i\frac{j}{L}z^{(\Delta t)}\right),\qquad z^{(\Delta t)}=\frac{1+s/ \Delta_{t}}{1-s/\Delta_{t}},\qquad s=\frac{z-1}{z+1}.$ $$

The output function $y^{(\Delta t)}$ is now equal to

$$ $y^{(\Delta t)}(z)=\sum_{j=0}^{L-1}\left([\texttt{FFT}(\mathbf{u})]_{j}\; \overline{G}(\omega_{j}^{(L,\Delta t)})\right)\text{exp}\left(-2\pi i\frac{j}{ L}z^{(\Delta t)}\right).$ $$

The only real difficulty is that when we sample $y^{(\Delta t)}$ at $z=\omega_{0}^{(L)},\ldots,\omega_{L-1}^{(L)}$ to obtain $\mathbf{y}^{(\Delta t)}$, we note that $z^{(\Delta t)}=\omega_{0}^{(L,\Delta t)},\ldots,\omega_{L-1}^{(L,\Delta t)}$ are not uniform on the unit circle. Hence, it cannot be achieved via inverse FFT. However, this sampling can be done via the so-called nonuniform FFT (NUFFT), which also takes $\mathcal{O}(L\log L)$. In general, one can interpret FFT as the procedure of evaluating a function

$$ $f(\omega)=\sum_{j=0}^{n-1}f_{j}\text{exp}\left(-j\omega\right),$ $$

at the degree-$L$ roots of unity $\bm{\omega}=(\omega_{0},\ldots,\omega_{L-1})$. The NUFFT is the procedure of evaluating exactly the same function $f$, but at potentially nonuniform nodes $\tilde{\bm{\omega}}\neq\bm{\omega}$. We remark that NUFFT is only used to simplify the representation of our derivation of Algorithm 1. It is not explicitly used in the algorithm, even though the numerical stability of the NUFFT procedure is studied in [, , ], and fast algorithms can be found in [, , ].

Using the NUFFT at the uneven samples $\bm{\omega}^{(L,\Delta t)}$, we obtain

$$ $\mathbf{y}^{(\Delta t)}=\texttt{NUFFT}(\texttt{FFT}(\mathbf{u})\circ\overline{ G}(\bm{\omega}^{(L,\Delta t)})).$ $$

Now, consider the following function

$$ $\tilde{y}^{(\Delta t)}(z)=y^{(\Delta t)}(z^{(1/\Delta t)})=\sum_{j=0}^{L-1} \left([\texttt{FFT}(\mathbf{u})]_{j}\;\overline{G}(\omega_{j}^{(L,\Delta t)}) \right)\text{exp}\left(-2\pi i\frac{j}{L}z\right).$ $$

This output function is a scaled version of $y^{(\Delta t)}$, where we scale the time domain by a factor of $1/\Delta t$. One can sample $\tilde{y}^{(\Delta t)}$ at $z=\omega_{0}^{(L)},\ldots,\omega_{L-1}^{(L)}$ using iFFT:

$$ $\tilde{\mathbf{y}}^{(\Delta t)}=\texttt{iFFT}(\texttt{FFT}(\mathbf{u})\circ \overline{G}(\bm{\omega}^{(L,\Delta t)})).$ $$

This leads us to eq. 7.

Figure: Figure 9: A visualization of appendix H.
Refer to caption: extracted/5896611/figures/NUFFT.png

## Appendix I Experimental details

In this section, we provide the details of the three experiments presented in the main text.

### I.1 Analyzing Hankel singular values using the sCIFAR-10 task

As mentioned in section 3, in the experiments presented in Figure 2 and Figure 5, we always train an SSM with $4$ layers and $128$ channels. Each channel in a layer is modeled by an LTI system with $n=64$ states. When we parameterize the LTI systems using $\mathbf{A},\mathbf{B},\mathbf{C}$, and $\mathbf{D}$, we assign a learning rate of $0.001$ to $\mathbf{A}$ and of $0.01$ to the rest. When we freeze the system matrices, then we set the learning rate of $\mathbf{A},\mathbf{B}$ and $\mathbf{C}$ to 0 0 while keeping that of $\mathbf{D}$ to be $0.01$. Note that the matrix $\mathbf{D}$ does not affect the Hankel singular values. All other model parameters are trained with a learning rate of $0.01$. For an LTI system parameterized by the Hankel matrix $\overline{\mathbf{H}}$, we adopt the same setting, except that $\overline{\mathbf{H}}$ is trained with a non-reduced learning rate of $0.01$. To compute the Hankel singular values of an LTI system $(\mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D})$, we use its balanced realization (see Appendix B). To compute the Hankel singular values of a system parameterized by $\overline{\mathbf{H}}$, we apply an SVD to the matrix $\overline{\mathbf{H}}$.

### I.2 Testing HOPE-SSMs long memory using noisy-sCIFAR

In this experiment (see Figure 6), we modify the sequential CIFAR-10 dataset by padding random sequences to the right. For each sequence of length $1024$ from the original dataset, we pad another sequence of length $1024$ to the end of it. The entries are sampled independently from the Gaussian distribution on the same magnitude as the entries in the original sequences. We adopt the same model architectures and learning rates as described in Section I.1 but make two exceptions. First, we fix the discretization size to be $\Delta t=0.1$, and therefore, it does not need a learning rate. In addition, a canonical SSM decodes the output sequence by first doing a pooling. Here, instead of pooling over all $2048$ output vectors, to make the problem more challenging and require longer memory, we only pool over the last $1024$ output vectors. These correspond to the output vectors when the noises are fed. We also test the models using different discretization sizes $\Delta t$. We see that when $\Delta t=1$, the S4D model fails to converge while our HOPE-SSM performs relatively well; on the other hand, when $\Delta t=0.01$, both models tend to have a relatively good performance. These observations align with our theory because a larger $\Delta t$ means that we put the discrete data on a continuous time domain with a longer span; hence, longer memory capacity is needed.

Figure: Figure 10: Performance of the HOPE-SSM and the S4D model on the noise-padded sCIFAR-10 task using different values of $\Delta t$.

### I.3 Hyperparameters of HOPE-SSMs in the Long-Range Arena

In this section, we present the table of hyperparameters used to train our HOPE-SSM on the LRA tasks [] (Apache License, Version 2.0). (See Table 2.) Our codes are adapted from the code associated with the original S4 and S4D papers [, ] (Apache License, Version 2.0). Note that compared to the hyperparameters used to train S4 and S4D, we use the same model hyperparameters and only slightly tune the training hyperparameters. All experiments are done on a NVIDIA A30 Tensor Core GPU with 24 GB of memory. The time efficiency of our model is roughly the same as that of the S4D model.

**Table 2: Configurations of the HOPE-SSM model, where DO, LR, BS, and WD stand for dropout rate, learning rate, batch size, and weight decay, respectively.**
| Task | Depth | #Features | Norm | Prenorm | DO | LR | BS | Epochs | WD | $\Delta$ Range |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| ListOps | 8 | 256 | BN | False | 0. | 0.01 | 20 | 100 | 0.03 | (0.001,0.1) |
| Text | 6 | 256 | BN | True | 0.01 | 0.01 | 16 | 150 | 0.05 | (0.001,0.1) |
| Retrieval | 6 | 128 | BN | True | 0. | 0.008 | 32 | 80 | 0.03 | (0.001,0.1) |
| Image | 6 | 128 | LN | False | 0.1 | 0.004 | 32 | 1500 | 0.01 | (0.001,10) |
| Pathfinder | 6 | 256 | BN | True | 0. | 0.001 | 16 | 250 | 0.03 | (0.0001,0.1) |
| Path-X | 6 | 128 | BN | True | 0. | 0.001 | 16 | 100 | 0.04 | (0.0001,1) |

## Peer Reviews

confidence: 3
contribution: 2
initial_rating: 5
presentation: 2
questions: 1. In section 3 before 3.1, how should I think about $\epsilon$? Under init3, if $\epsilon=$1e-4, then the frozen model has a higher numerical rank, since there is considerable mass below 1e-4 for the trained model, but the frozen model has worse performance?
2. Why is numerical rank the correct metric instead of the sum of the smallest eigenvalues? ROM quantifies the approximation error using the latter.
3. Theorem 2's upper and lower bounds are off by n, which is significant when the hidden dimension is large (regime of interest for high-order LTIs). In addition, Theorem 4 has a $\sqrt{n}$ on the right hand side, but it's not clear whether the worse $n$ scaling for default SSM parameterization is tight. Can authors comment on this?
soundness: 2
strengths: This paper proposes a theoretical framework for analyzing SSMs, and there are many supporting graphs and experiments. The plots are well-executed, and the experiments are LRA has a good set of baselines. I also appreciate the noise-padded sCIFAR-10 experiment as an ablation. The algorithm is clearly written.
summary: This work proposes an alternative parameterization for SSMs using the Hankel operator of the dynamical system, called HOPE. HOPE learns the Hankel operator of a linear dynamical system, instead of directly learning the dynamics matrices as in other SSMs. By using techniques and results in Hankel operator theory, authors show that HOPE is more robust to noise and has long memory. Finally, experiments on LRA show that HOPE is competitive with other SSMs and outperforms baselines on several tasks.
weaknesses: 1. Motivation: random init in models like Mamba also works, so the initialization issues seem to be limited to certain SSMs.
2. HOPE-SSM does not have non-decaying memory: these are asymptotically stable systems, instead of marginally stable systems, so they cannot have arbitrarily long memory. This point is not clarified upfront in the paper.
3. HOPE-SSM on Path-X, which is the task with the longest memory in LRA, does not outperform S5. This raises some concerns on the long-memory capability of HOPE-SSM in practice.
4. Though authors provide expressivity results on HOPE, there is no guarantee on learning with this parameterization (convergence, etc). 

Please also see questions.

---

confidence: 2
contribution: 3
initial_rating: 8
presentation: 3
questions: 1. If Mamba achieves SOTA performance, could you explain why there is no experimental comparison with it? If there are no technical limitations and the comparison is meaningful, including it could further enhance the quality of the paper.

2. Could you please explain how Theorem 4 guarantees the numerical stability of HOPE?

3. In Figure 6 (right), why does the memory decay of HOPE still appear significantly smaller than that of S4D after $t=n=64$? Shouldn't HOPE have no memory after $n$?

4. The Hankel matrix looks very similar to the SSM convolution kernel. Could you please explain the differences between the two?
soundness: 3
strengths: 1. I believe HOPE is very well motivated by Section 3. The analysis presented in that section clearly illustrates all the issues with current methods based on standard parametrization. In a nutshell, Figures 2 and 3 perfectly summarize the entire section.

2. The three advantages of HOPE are well explained and essential. It is impressive that the authors managed to achieve all three benefits with a simple new parametrization.

3. The entire framework seems very simple, yet powerful.
summary: The paper first examines the stability of SSM training and the recently used initialization method through the lens of Hankel operator theory. The authors argue that by decreasing the effective rank of the Hankel operator, the system loses its expressivity for modeling complex sequences. They then propose a new parametrization of LTI systems based on the Hankel operator, which improves robustness and performance due to the stable spectrum of the Hankel operator during training. Finally, they demonstrate that with this parametrization, memory remains practically constant within a fixed time window. They test this framework on the sCIFAR-10 dataset.
weaknesses: Some results lack sufficient motivation. For example, the setup introduced before Theorem 1 could be better explained—why is it reasonable to sample entries of $A$ from $F_a$​ and entries of $\overline{B}\circ\overline{C}^\top$ from a standard normal distribution? Similarly, the discussion preceding Algorithm 1, as well as the algorithm itself, appears more complicated than necessary, with essentially the same information repeated three times. I would prefer having only Algorithm 1 with the full expressions for variables, as in the current prelude to Algorithm 1.

---

confidence: 3
contribution: 3
initial_rating: 8
presentation: 4
questions: 1. How does the parameterisation behave in tasks where a recency bias can be beneficial? For example this is the case for generative Language Modelling.
2. Can the authors comment on the relationship between the learned discretization parameter $\Delta_t$ and the size of the window of non-discarded inputs of the HOPE SSM? Is the full input always considered? How does this change during training? Such discussion and potential experimental evidence could strengthen the paper.

Comments:
1. Information on the actual runtime could be moved the main paper since now it is quite hard to find (it is in the last page of the appendix). 
2. The paper could benefit from having a discussion on how HOPE differs from spectral SSMs (Agarwal et al.2023) and the stable parameterisation proposed by Wang & Li (2023), which are probably the closest works This could be placed in the appendix. Such discussion could potentially address pros and cons of each approach in addition to the brief discussion in the introduction.

Minor: 
- Line 159-161: Can the authors provide a definition of the function space $\ell^2(\mathbb{N})$   and $L^2([0,\infty])$?
- Line 208: imaginal axis -> imaginary axis
soundness: 4
strengths: 1. Strong empirical and theoretical evidence in favour of the proposed parameterisation.
2. Clearly written.
3. The algorithm is novel, albeit SSM working in the frequency domain have been proposed by (Agarwal et al. 2023).
summary: This paper proposes a new parametrisation, named HOPE, for Linear Time Invariant (LTI) State Space Models (SSMs), based on the Hankel operator, that improves the performance and makes training more stable. The proposed parametrisation reduces the SSM parameters to 1/3 of the original and works in the frequency domain, hence the Fast Fourier Transform (FFT) and its inverse need to be used during the forward pass. Fortunately, the time and memory complexity of the forward pass remains the same as that of standard LTI SSMs using the convolutional form. The authors first show how current parameterisations can perform poorly both at initialisation or during training and how this can be explained by their corresponding Henkel operator being low-rank (Theorem 1, Figure 3) and the parameterisation being numerically unstable (Theorem 2, Figure 4). Then, they show how their parameterisation effectively mitigates both issues (Theorems 3-4, Figure 5) and does not have an exponentially decaying memory (Advantage III). Each claim Is supported both theoretically and empirically. Empirically, the authors test different parameterisation, initialisation and training schemes on the (noisy-padded) sequential CIFAR dataset, and on the Long Range Arena, also showing how HOPE outperforms previous LTI parameterisations in terms of Test Accuracy.
weaknesses: 1. Some limitations are not experimentally investigated. In particular, the method seems to weigh the input time steps belonging to a fixed window equally and discard the others. This could be a limitation but all tasks considered in the experiments do not seem to benefit from a recency bias (see also the first question). 
2. Limited discussion of Related works (Minor).

---

confidence: 3
contribution: 3
initial_rating: 5
presentation: 3
questions: * Why do the authors choose to show the Hankel singular values of SSMs trained after 10 epochs instead of after the whole training process (e.g. Figure 3 & 5)?

* How do the authors validate the assumption on the distribution of the diagonal entry of $\bar{A}$ in Theorem 1? For linear recurrent unit [1], $\bar{A}$ is initialized such that $a_j$ is uniformly distributed in a disc $[r_{\min}, r_{\max}]$ where $r_{\min}, r_{\max}$ are both close to $1$.

* Also for Theorem 1, the theoretical bound does not indicate $\beta < 1$, but $\beta \leq 1$ instead. Is there a way to mitigate the logarithmic term in theory?





[1] Orvieto, Antonio, et al. "Resurrecting recurrent neural networks for long sequences." International Conference on Machine Learning. PMLR, 2023.
soundness: 2
strengths: * The paper is written clearly and the ideas have been presented in the proper order.

* The idea using Hankel operator theory to parameterize state space models is novel, which can provide more insights for other recurrent architectures on sequence modeling.
summary: This paper develops a new parameterization scheme to help improve training of state space models (SSMs) in sequence modeling. The new parameterization method is based on Hankel operator theory, which resolves the low-rank and numerical instability issues of an LTI system, and the parameterized SSMs could enjoy non-decaying memory. Numerical experiments are conducted on the long range arena benchmark.
weaknesses: * The connections of the Hankel rank/stability and the SSM performance seem not that positively correlated in this paper. For example, from Figure we see that the numerical rank of HOPE-SSM is even lower than init3 before and after training, but these two initializations produce similar numerical performance on the sCIFAR-10 task. Does it mean that even slightly low rank Hankel matrix can achieve good results? 

* From Algorithm1, the Hankel operator method looks very similar to the convolution method, it would be better if the authors could add more discussion on these two methods. For example, [1] can be a good reference. 

* The notations for $A$ and $\bar{A}$ are not consistently used in this paper. From my understanding, $A$ is the hidden matrix for the continuous system and $\bar{A}$ is the hidden matrix for the discrete system. The HiPPO-Legs initialization is for $A$ matrix, and $\bar{A}$ contains the information of the timescale $\Delta t$. It would be better to distinguish these two notations when making statements. For example, in Figure 4, the perturbation is added to the continuous system while the Hankel singular values are calculated  for the discrete system, and the discrete system contains the information for $\Delta t$. It seems that $\Delta t$ is set to be $1$ in this paper, but in practice, the timescale is scaled based on the sequence length. So what is the reason to take $\Delta t = 1$?
Also, it would be better if the authors could add more details on how the perturbation is added here.

[1] Li, Yuhong, et al. "What Makes Convolutional Models Great on Long Sequence Modeling?." The Eleventh International Conference on Learning Representations.

---

confidence: 4
contribution: 2
initial_rating: 5
presentation: 4
questions: 1. Can the authors clarify what is the fundamental difference, if any, of the HOPE-SSM and classical convolutional models that directly parameterises a truncated impulse response function? Looking at Alg 1, the only possible departure appears to be in the for loop in steps 6-7. 
2. Can the authors clarify how (or if) the proposed model can process long sequences without increasing the filter size, while still meaningfully extract memory structures? 
3. If different, can the authors comment and demonstrate both theoretically and numerically their advantage against convolutional-type sequence models (e.g. WaveNet/temporal CNNs, and directly paramterising the impulse response function)?
4. Is the conclusion of Theorem 1 verified by results in Fig 3 quantitatively?

Post rebuttal/revision: Appreciate the comparison with discrete convolutional models to assess the model's performance. I think while the novelty is still limited (it is at most a rescaled temporal convolution model with large filters), the performance appears good and the paper is well written. I have increased my score.
soundness: 2
strengths: While the paper's identified limitation (i.e. favoring exponential decay in memory)  of SSM is widely known and studied, its use of the Hankel numerical rank to study it is interesting. The theoretical results are carefully derived and the practical results show good promise. Finally, the paper is carefully written and the explanations of their approach and results are excellent.
summary: This paper proposes an alternative parameterisation of the state space models, where the SSM layer's dynamics is replaced by parameterising rows of the Hankel matrix. It is shown that classical SSM parameterisation tends to give rise to Hankel matrices with fast decaying Hankel singular values, and is posited to correlate with poor performance for long memory tasks. The proposed parameterisation does not involve repeated multiplication of the state space matrix, hence avoiding this shortcoming. These are supported by theoretical results, and the proposed method is shown to improve upon existing SSM variants on sCIFAR and long range arena.
weaknesses: The main weakness, to my understanding, is the premise of the proposed approach. In classical LTI and SSM literature on system identification (e.g. Tangirala, A. K. *Principles of System Identification: Theory and Practice*), the very reason why SSMs were introduced in favor of convolutional models is its ability to handle infinite sequences with a finite parameterisation. The present approach appears to go in the reverse direction:

1. In particular, I cannot see what is the difference between HOPE-SSM and the classical convolutional model (the most basic sequence to sequence model in LTI literature, where one directly parameterises the impulse response function), with the filter size $n$ is chosen as the same the original SSM's hidden state dimension. In fact,  the parameterisation (5) via $h$ seems to be exactly the convolution filter (impulse response), since each row of the Hankel matrix is just the impulse response (shifted by 1 for each subsequent row). If this were the case, the method has little to do with Hankel matrix, but is rather just a convolutional model based on parameterising the impulse response function.
2. The HOPE-SSM model's size scales with sequence length if it must process long sequences. In particular, all memory functions that can be represented are finitely supported, so it is not accurate to say that this improves upon SSM parameterisation, which is used in preference to convolutional models *precisely* because it allows the parameterisation of memory of infinite support with finite number of parameters (albeit introducing curse of memory due to forcing an exponential decay, as observed in previous literature). 
3. The observed strengths of HOPE-SSM (non-exponential decay) are also present in convolutional networks (e.g. WaveNet with small filters but multiple layers, or single layer direct parameterisation of the impulse response), thus the numerical comparison should be made with those models, instead of just SSM variants.

Meta-review: The paper introduces HOPE, a novel parameterization scheme for linear, time-invariant SSMs, leveraging Hankel operator theory to enhance initialization and training stability. The reviewers praised the robust theoretical formulation and strong empirical results demonstrating improved performance on Long-Range Arena asks and sequential CIFAR-10. Strengths highlighted include the insightful use of Hankel operators, non-decaying memory within a fixed time window, and practical improvements over existing HiPPO-based SSMs. The paper is well-executed and addresses critical challenges in training SSMs. Some reviewers raised minor concerns about presentation clarity and specific implementation details. The authors' rebuttal successfully addressed these points, and no major concerns remained post-discussion. Overall, the combination of strong theoretical grounding, novel insights, and empirical improvements makes this a valuable contribution. I recommend acceptance.

Decision: Accept (Poster)