Tuning Frequency Bias of State Space Models

## Abstract

Abstract State space models (SSMs) leverage linear, time-invariant (LTI) systems to effectively learn sequences with long-range dependencies.
By analyzing the transfer functions of LTI systems, we find that SSMs exhibit an implicit bias toward capturing low-frequency components more effectively than high-frequency ones.
This behavior aligns with the broader notion of frequency bias in deep learning model training.
We show that the initialization of an SSM assigns it an innate frequency bias and that training the model in a conventional way does not alter this bias.
Based on our theory, we propose two mechanisms to tune frequency bias: either by scaling the initialization to tune the inborn frequency bias; or by applying a Sobolev-norm-based filter to adjust the sensitivity of the gradients to high-frequency inputs, which allows us to change the frequency bias via training.
Using an image-denoising task, we empirically show that we can strengthen, weaken, or even reverse the frequency bias using both mechanisms.
By tuning the frequency bias, we can also improve SSMs’ performance on learning long-range sequences, averaging an 88.26 % percent 88.26 88.26\% 88.26 % accuracy on the Long-Range Arena (LRA) benchmark tasks.

## 1 Introduction

Sequential data are ubiquitous in fields such as natural language processing, computer vision, generative modeling, and scientific machine learning.
Numerous specialized classes of sequential models have been developed, including recurrent neural networks
(RNNs) [, , , , , ], convolutional neural networks (CNNs) [, ], continuous-time models (CTMs) [, ], transformers [, , , , ], state space models (SSMs) [, , , ], and Mamba [, ]. Among these, SSMs stand out for their ability to learn sequences with long-range dependencies.

Using the continuous-time linear, time-invariant (LTI) systems,

$$ $\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) $$

where $\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}$, an SSM computes the output time-series $\mathbf{y}(t)$ from the input $\mathbf{u}(t)$ via a latent state vector $\mathbf{x}(t)$. Compared to an RNN, a major computational advantage of an SSM is that the LTI system can be trained both efficiently (i.e., the training can be parallelized for long sequences) and numerically robustly (i.e., it does not suffer from vanishing and exploding gradients). An LTI system can be computed in the time domain via convolution:

$$ $\mathbf{y}(t)=(\mathbf{h}*\mathbf{u}+\mathbf{D}\mathbf{u})(t)=\int_{-\infty}^{ \infty}\mathbf{h}(t-\tau)\mathbf{u}(\tau)d\tau+\mathbf{D}\mathbf{u}(t),\qquad \mathbf{h}(t)=\mathbf{C}\text{exp}(t\mathbf{A})\mathbf{B}.$ $$

Alternatively, it can be viewed as an action in the frequency domain:

$$ $\hat{\mathbf{y}}(s)=\mathbf{G}(is)\hat{\mathbf{u}}(s),\qquad\mathbf{G}(is):= \mathbf{C}(is\mathbf{I}-\mathbf{A})^{-1}\mathbf{B}+\mathbf{D},\quad s\in \mathbb{R},$ (2) $$

where $i$ is the imaginary unit and $\mathbf{I}$ is the identity matrix. The function $\mathbf{G}$ is called the transfer function of the LTI system.(^1^11Equation 2 often appears in the form of the Laplace transform instead of the Fourier transform. We restrict ourselves to the Fourier transform, due to its widespread familiarity, by assuming decay properties of $\mathbf{u}(t)$. All discussions in this paper nevertheless apply to the Laplace domain as well.) It is a rational function whose poles are at the spectrum of $\mathbf{A}$.

The frequency-domain characterization of the LTI systems in eq. 2 sets the stage for understanding the so-called frequency bias of an SSM.
The term “frequency bias” originated from the study of a general overparameterized multilayer perceptron (MLP) [], where it was observed that the low-frequency content was learned much faster than the high-frequency content.
It is a form of implicit regularization [].
Frequency bias is a double-edged sword: on one hand, it partially explains the good generalization capability of deep learning models, because most high-frequency noises are not learned until the low-frequency components are well-captured; on the other hand, it puts a curse on learning the useful high-frequency information in the target.

In this paper, we aim to understand the frequency bias of SSMs. In Figure 1, we observe that, similar to most deep learning models, SSMs are also better at learning the low frequencies than the high ones.
To understand that, we develop a theory that connects the spectrum of $\mathbf{A}$ to the SSM’s capability of processing high-frequency signals.
Then, based on the spectrum of $\mathbf{A}$, we analyze the frequency bias in two steps.
First, we show that the most popular initialization schemes [, , ] lead to SSMs that have an innate frequency bias.
More precisely, they place the spectrum of $\mathbf{A}$, $\Lambda(\mathbf{A})$, in the low-frequency region in the $s$-plane, preventing LTI systems from processing high-frequency input, regardless of the values of $\mathbf{B}$ and $\mathbf{C}$.
Second, we consider the training of the SSMs.
Using the decay properties of the transfer function, we show that if an eigenvalue $a_{j}\in\Lambda(\mathbf{A})$ is initialized in the low-frequency region, then its gradient is insensitive to the loss induced by the high-frequency input content.
Hence, if an SSM is not initialized with the capability of handling high-frequency inputs, then it will not be trained to do so by conventional training.

Figure: Figure 1: In a synthetic example to illustrate the frequency bias of SSMs, we form the inputs by superposing three waves of low, moderate, and high frequencies, respectively. We train an S4D model to regress the magnitudes of the three waves. We observe that the magnitudes of the low-frequency waves can be approximated much better compared to those of the high-frequency waves. In Figure 9, we show how to tune the frequency bias in this example.

The initialization of the LTI systems equip an SSM with a certain frequency bias, but this is not necessarily the appropriate implicit bias for a given task.
Depending on whether an SSM needs more expressiveness or generalizability, we may want less or more frequency bias, respectively (see Figure 9).
Motivated by our analysis, we propose two ways to tune the frequency bias:

- 1.
Instead of using the HiPPO initialization, we can scale the initialization of $\Lambda(\mathbf{A})$ to lower or higher-frequency regions. This serves as a “hard tuning strategy” that marks out the regions in the frequency domain that can be learned by our SSM.
- 2.
Motivated by the Sobolev norm, which applies weights to the Fourier domain, we can apply a multiplicative factor of $(1+|s|)^{\beta}$ to the transfer function $\mathbf{G}(is)$. This is a “soft tuning strategy” that reweighs each location in the frequency domain. By selecting a positive or negative $\beta$, we make the gradients more or even less sensitive to the high-frequency input content, respectively, which changes the frequency bias during training.

One can think of these two mechanisms as ways to tune frequency bias at initialization and during training, respectively.
After rigorously analyzing them, we present an experiment on image-denoising with different noise frequencies to demonstrate their effectiveness.
We also show that tuning the frequency bias enables better performance on tasks involving long-range sequences.
Equipped with our two tuning strategies, a simple S4D model can be trained to average an $88.26\%$ accuracy on the Long-Range Arena (LRA) benchmark tasks [].

Contribution.
Here are our main contributions:

- 1.
We formalize the notion of frequency bias for SSMs and quantify it using the spectrum of $\mathbf{A}$. We show that a diagonal SSM initialized by HiPPO has an innate frequency bias.
We are the first to study the training of the state matrix $\mathbf{A}$, and we show that training the SSM does not alter this frequency bias.
- 2.
We propose two ways to tune frequency bias, by scaling the initialization, and by applying a Sobolev-norm-based filter to the transfer function of the LTI systems.
We study the theory of both strategies and provide guidelines for using them in practice.
- 3.
We empirically demonstrate the effectiveness of our tuning strategies using an image-denoising task.
We also show that tuning the frequency bias helps an S4D model to achieve state-of-the-art performance on the Long-Range Arena tasks and provide ablation studies.

To make the presentation cleaner, throughout this paper, we focus on a single single-input/single-output (SISO) LTI system $\Gamma=(\mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D})$ in an SSM, i.e., $m=p=1$, although all discussions naturally extend to the multiple-input/multiple-output (MIMO) case, as in an S5 model []. Hence, the transfer function $\mathbf{G}:\mathbb{C}\rightarrow\mathbb{C}$ is complex-valued.
We emphasize that while we focus on a single system $\Gamma$, we do *not* isolate it from a large SSM; in fact, when we study the training of $\Gamma$ in section 4, we backpropagate through the entire SSM.

Related Work.
The frequency bias, also known as the spectral bias, of a general neural network (NN) was initially observed and studied in [, , ]. The name spectral bias stemmed from the spectral decomposition of the so-called neural tangent kernels (NTKs) [], which provides a means of approximating the training dynamics of an overparameterized NN [, , ]. By carefully analyzing the eigenfunctions of the NTKs, [, ] proved the frequency bias of an overparameterized two-layer NN for uniform input data. The case of nonuniform input data was later studied in [, ]. The idea of Sobolev-norm-based training of NNs has been considered in [, , , , , , , ].

The initialization of the LTI systems in SSMs plays a crucial role, which was first observed in []. Empirically successful initialization schemes called “HiPPO” were proposed in [, , ]. Other efforts in improving the initialization of an SSM were studied in [, ]. Later, [, ] attributed the success of HiPPO to the proximity of the spectrum of $\mathbf{A}$ to the imaginary axis (i.e., the real parts of the eigenvalues of $\mathbf{A}$ are close to zero). This paper considers the imaginary parts of the eigenvalues of $\mathbf{A}$, which was also discussed in the context of the approximation-estimation tradeoff in []. The training of SSMs has mainly been considered in [, ], where the matrix $\mathbf{A}$ is assumed to be fixed, making the optimization convex. To our knowledge, we are the first to consider the training of $\mathbf{A}$. While we consider the decay of the transfer functions of the LTI systems in the frequency domain, there is extensive literature on the decay of the convolutional kernels in the time domain (i.e., the memory) [, , , , , ].

## 2 What is the Frequency Bias of an SSM?

Figure: Figure 2: The frequency bias of an SSM says that the frequency response has more variation in the low-frequency area than the high-frequency one.

In Figure 1, we see an example where an S4D model is better at predicting the magnitude of a low-frequency component in the input than a high-frequency one.
This coincides with our intuitive interpretation of frequency bias: the model is better at “handling” low frequencies than high frequencies.
To rigorously analyze this phenomenon for SSMs, however, we need to formalize the notion of frequency bias.
This is our goal in this section. One might imagine that an SSM has a frequency bias if, given a time-series input $\mathbf{u}(t)$ that has rich high-frequency information, its time-series output $\mathbf{y}(t)$ lacks high-frequency content.
Unfortunately, this is not the case: an SSM is capable of generating high-frequency outputs.
Indeed, the skip connection $\mathbf{D}$ of an LTI system is an “all-pass” filter, multiplying the whole input $\mathbf{u}(t)$ by a factor of $\mathbf{D}$ and adding it to the output $\mathbf{y}(t)$.
On the other hand, the secret of a successful SSM hides in $\mathbf{A}$, $\mathbf{B}$, and $\mathbf{C}$ [].
In an ablation study, when $\mathbf{D}$ is removed, an S4D model only loses less than $2\%$ of accuracy on the sCIFAR-10 task [, ], whereas the model completely fails when we remove $\mathbf{A}$, $\mathbf{B}$, and $\mathbf{C}$.
This can be ascribed to the LTI system’s power to model complicated behaviors in the frequency domain.
That is, each Fourier mode in the input has its own distinct “pass rate” (see [, , ] for why this is an important feature of LTI systems).
For example, the task in Figure 1 can be trivially solved if the LTI system can filter out a single mode $a_{1}\cos(x)$, $a_{2}\cos(16x)$, or $a_{3}\cos(256x)$ from the superposition of the three; the skip connection $\mathbf{D}$ alone is not capable of doing that.

Given that, we can formulate the frequency bias of an LTI system as follows (see Figure 2):
{adjustwidth}15pt15pt
Frequency bias of an SSM means that the frequency responses (i.e., the transfer functions $\mathbf{G}$) of LTI systems have more variation in the low-frequency area than the high-frequency area.
More precisely, given the transfer function $\mathbf{G}(is)$, we can study its total variation in a particular interval $[a,b]$ in the Fourier domain defined by

$$ $V_{a}^{b}(\mathbf{G}):=\sup_{\begin{subarray}{c}a=s_{0}<s_{1}<\cdots<s_{N}=b, \\ N\in\mathbb{N}\end{subarray}}\sum_{j=1}^{N}|\mathbf{G}(is_{j})-\mathbf{G}(is_{ j-1})|,\qquad-\infty\leq a<b\leq\infty.$ $$

Intuitively, $V_{a}^{b}(\mathbf{G})$ measures the total change of $\mathbf{G}(is)$ when $s$ moves from $a$ to $b$. The larger it is, the better an LTI system is at distinguishing the Fourier modes with frequencies between $a$ and $b$. Frequency bias thus says that for a fixed-length interval $[a,b]$, $V_{a}^{b}(\mathbf{G})$ is larger when $[a,b]$ is near the origin than when it lies in the high-frequency region, i.e., when it is far from the origin.

## 3 Frequency Bias of an SSM at Initialization

Our exploration of the frequency bias of an SSM starts with the initialization of a SISO LTI system $(\mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D})$, where $\mathbf{A}=\text{diag}(a_{1},\ldots,a_{n})\in\mathbb{C}^{n\times n}$ is diagonal.
The system is assumed to be stable, meaning that $a_{j}=x_{j}+iy_{j}$ for some $x_{j}<0$ for all $1\leq j\leq n$, where $x_{j}$ and $y_{j}$ are the real and the imaginary parts of $a_{j}$, respectively.
Note that the diagonal structure of $\mathbf{A}$ is indeed the most popular choice for SSMs [, ], in which case it suffices to consider the Hadamard (i.e., entrywise) product $\mathbf{B}\circ\mathbf{C}^{\top}=\begin{bmatrix}c_{1}&\cdots&c_{n}\end{bmatrix
}^{\top}\in\mathbb{C}^{n}$, where $c_{j}=\xi_{j}+i\zeta_{j}$ for all $1\leq j\leq n$. Then, the transfer function $\mathbf{G}$ is naturally represented in partial fractions:

$$ $\mathbf{G}(is)=\frac{c_{1}}{is-a_{1}}+\cdots+\frac{c_{n}}{is-a_{n}}+\mathbf{D} =\frac{\xi_{1}+i\zeta_{1}}{-x_{1}+i(s-y_{1})}+\cdots+\frac{\xi_{n}+i\zeta_{n}} {-x_{n}+i(s-y_{n})}+\mathbf{D}.$ $$

In most cases, the input $\mathbf{u}(t)$ is real-valued.
To ensure that output is real-valued, the standard practice is to take the real part of $\mathbf{G}$ as a real-valued transfer function before applying eq. 2:

$$ $\tilde{\mathbf{G}}(is):=\text{Re}(\mathbf{G}(is))=\sum_{j=1}^{n}\frac{\zeta_{j }(s-y_{j})-\xi_{j}x_{j}}{x_{j}^{2}+(s-y_{j})^{2}}+\mathbf{D}.$ (3) $$

We now derive a general statement for the total variation of $\tilde{\mathbf{G}}$ given the distribution of $y_{j}$.

###### Lemma 1 .

Let $\tilde{\mathbf{G}}$ be the transfer function defined in eq. 3. Given any $B>\max_{j}|y_{j}|$, we have

$$ $V_{-\infty}^{-B}(\tilde{\mathbf{G}})\leq\sum_{j=1}^{n}\frac{\left|c_{j}\right| }{\left|y_{j}+B\right|},\qquad V_{B}^{\infty}(\tilde{\mathbf{G}})\leq\sum_{j=1 }^{n}\frac{\left|c_{j}\right|}{\left|y_{j}-B\right|}.$ $$

Given the formula of the transfer function, the proof of 1 is almost immediate.
In this paper, we leave all proofs to Appendix A and C. 1 illustrates a clear and intuitive concept:
{adjustwidth}15pt15pt
If the imaginary parts of $a_{j}$ are distributed in the low-frequency region, i.e., $\left|y_{j}\right|$ are small, the transfer function has a small total variation in the high-frequency areas $(-\infty,-B]$ and $[B,\infty)$, inducing a frequency bias of the SSM.
We can now apply 1 to study the innate frequency bias of the HiPPO initialization [].
While there are many variants of HiPPO, we choose the one that is commonly used in practice [].
All other variants can be similarly analyzed.

###### Corollary 1 .

Assume that $a_{j}=-0.5+i(-1)^{j}\lfloor j/2\rfloor\pi$ and $\xi_{j},\zeta_{j}\sim\mathcal{N}(0,1)$ i.i.d., where $\mathcal{N}(0,1)$ is the standard normal distribution. Then, given $B>n\pi/2$ and $\delta>0$, we have

$$ $V_{-\infty}^{-B}(\tilde{\mathbf{G}}),V_{B}^{\infty}(\tilde{\mathbf{G}})\leq \frac{\sqrt{2n}(\sqrt{n}+\sqrt{\ln(1/\delta)})}{B-n/2}\qquad\text{with probability}\geq 1-\delta.$ $$

In particular, 1 tells us that the HiPPO initialization only captures the frequencies $s\in[-B,B]$ up to $B=\mathcal{O}(n)$, because when $B=\omega(n)$, we see that $V_{-\infty}^{-B}(\tilde{\mathbf{G}}),V_{B}^{\infty}(\tilde{\mathbf{G}})$ vanish as $n$ increases.
This means that no complicated high-frequency responses can be learned.

###### Lemma 1 .

###### Corollary 1 .

## 4 Frequency Bias of an SSM during Training

In section 3, we see that the initialization of the LTI systems equips an SSM with an innate frequency bias. A natural question to ask is whether an SSM can be trained to adopt high-frequency responses. Analyzing the training of an SSM (or many other deep learning models) is not an easy task, and we lack theoretical characterizations. Two notable exceptions are [, ], where the convergence of a trainable LTI system to a target LTI system is analyzed, assuming that the state matrix $\mathbf{A}$ is fixed to make the optimization problem convex. Unfortunately, this assumption is too strong to be applied for our purpose. Indeed, 1 characterizes the frequency bias using the distribution of $y_{j}=\text{Im}(a_{j})$, making the training dynamics of $\mathbf{A}$ a crucial element in our analysis. Even if we set aside the issue of $\mathbf{A}$, analyzing an isolated LTI system in an SSM remains unrealistic: when an SSM, consisting of hundreds of LTI systems, is trained for a single task, there is no clear notion of “ground truth” for each individual LTI system within the model.

To make our discussion truly generic, we assume that there is a loss function $\mathcal{L}(\bm{\Theta})$ that depends on all parameters $\bm{\Theta}$ of an SSM. In particular, $\bm{\Theta}$ contains $x_{j},y_{j},\xi_{j},\zeta_{j}$, and $\mathbf{D}$ from every LTI system within the SSM, as well as the encoder, decoder, and inter-layer connections. With mild assumptions on the regularity of the loss function $\mathcal{L}$, we provide a quantification of the gradient of $\mathcal{L}$ with respect to $y_{j}$ that leads to a qualitative statement about the frequency bias during training.

###### Theorem 1 .

Let $\mathcal{L}(\bm{\Theta})$ be a loss function and $(\mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D})$ be a diagonal LTI system in an SSM defined in section 3. Let $\tilde{\mathbf{G}}$ be its associated real-valued transfer function defined in eq. 3. Suppose the functional derivative of $\mathcal{L}(\bm{\Theta})$ with respect to $\tilde{\mathbf{G}}(is)$ exists and is denoted by $(\partial/\partial\tilde{\mathbf{G}}(is))\mathcal{L}$. Then, if $|{(\partial/\partial\tilde{\mathbf{G}}(is))\mathcal{L}}|=\mathcal{O}(|s|^{p})$ for some $p<1$, we have

$$ $\frac{\partial\mathcal{L}}{\partial y_{j}}=\int_{-\infty}^{\infty}\frac{ \partial\mathcal{L}}{\partial\tilde{\mathbf{G}}(is)}\cdot K_{j}(s)\;ds,\qquad K _{j}(s):=\frac{\zeta_{j}((s-y_{j})^{2}-x_{j}^{2})-2\xi_{j}x_{j}(s-y_{j})}{[x_{ j}^{2}+(s-y_{j})^{2}]^{2}},$ (4) $$

for every $1\leq j\leq n$. In particular, we have that $\left|K_{j}(s)\right|=\mathcal{O}\left(|\zeta_{j}s^{-2}|+|\xi_{j}s^{-3}|\right)$ as $|s|\rightarrow\infty$.

In 1, we use a technical tool called the functional derivative []. The assumption that $(\partial/\partial\tilde{\mathbf{G}}(is))\mathcal{L}$ exists is easily satisfied, and we leave a survey of functional derivatives to Appendix B. The assumption that $|{(\partial/\partial\tilde{\mathbf{G}}(is))\mathcal{L}}|$ grows at most sublinearly is to guarantee the convergence of the integral in eq. 4; it is also easily satisfiable. We will see that the growth/decay rate of $|{(\partial/\partial\tilde{\mathbf{G}}(is))\mathcal{L}}|$ plays a more important role when we start to tune the frequency bias using the Sobolev-norm-based method (see section 5.2). As usual, one can intuitively think of the functional derivative $(\partial/\partial\tilde{\mathbf{G}}(is))\mathcal{L}$ as a measurement of the “sensitivity” of the loss function $\mathcal{L}$ to an LTI system’s action on a particular frequency $s$ (i.e., $\tilde{\mathbf{G}}(is)$). The fact that it is multiplied by a factor of $K_{j}(s)$ in the computation of the gradient in eq. 4 conveys the following important message:
{adjustwidth}15pt15pt
The gradient of $\mathcal{L}$ with respect to $y_{j}$ highly depends on the part of the loss that has “local” frequencies near $s=y_{j}$. It is relatively unresponsive to the loss induced by high frequencies, with a decaying factor of $\mathcal{O}(|s|^{-2})$ as the frequency increases, i.e., as $|s|\rightarrow\infty$.
Hence, the loss landscape of the frequency domain contains many local minima, and an LTI system can rarely learn the high frequencies with the usual training. To verify this, we train an S4D model initialized by HiPPO to learn the sCIFAR-10 task for $100$ epochs. We measure the relative change of each parameter ${\theta}$: $\Delta(\bm{\theta})=(|{\theta}^{(0)}\!-\!{\theta}^{(100)}|)/(|{\theta}^{(0)}|)$, where the superscripts indicate the epoch number. As we will show in section 6, the HiPPO initialization is unable to capture the high frequencies in the CIFAR-10 pictures fully. From Table 1, however, we see that $\text{Im}(\text{diag}(\mathbf{A}))$ is trained very little: every $y_{j}$ is only shifted by $1.43\%$ on average. This can be explained by 1: $y_{j}$ is easily trapped by a low-frequency local minimum.

**Table 1: The average relative change of each LTI system matrix in an S4D model trained on the sCIFAR-10 task. We see that the imaginary parts of $\text{diag}(\mathbf{A})$ are almost unchanged during training.**
| Parameter | $\text{Re}(\text{diag}(\mathbf{A}))$ | $\text{Im}(\text{diag}(\mathbf{A}))$ | $\mathbf{B}\circ\mathbf{C}^{\top}$ | $\mathbf{D}$ |
| --- | --- | --- | --- | --- |
| $\Delta$ | $1002.705$ | $0.0143$ | $1.1801$ | $0.8913$ |

An Illustrative Example.
Our analysis of the training dynamics of ${y}_{j}$ in 1 is very generic, relying on the notion of the functional derivatives. To make the theorem more concrete, we consider a synthetic example (see Figure 3). We fall back to the case of approximating a target function

$$ $\tilde{\mathbf{F}}(is)=\text{Re}\left(\frac{5}{is-(-1-50i)}+\frac{0.2}{is-(-1+ 50i)}+0.01\cos\left(\frac{9}{4}s\right)\cdot\mathbbm{1}_{[-2\pi,2\pi]}\right), \qquad s\in\mathbb{R},$ $$

using a trainable $\tilde{\mathbf{G}}$, where $9/4$ is chosen to guarantee the continuity of $\tilde{\mathbf{F}}$. We set the number of states to be one, i.e., $n=1$. For illustration purposes, we fix $x=-1$ and $\zeta=0$; therefore, we have

$$ $\tilde{\mathbf{G}}(is)=\text{Re}\left(\frac{\xi}{is-(-1-yi)}\right),\qquad s \in\mathbb{R},$ $$

where our only trainable parameters are $y$ and $\xi$. Our target function $\tilde{\mathbf{F}}$ contains two modes and some small noises between $-2\pi$ and $2\pi$, whereas $\tilde{\mathbf{G}}$ is unimodal with a trainable position and height (see Figure 3 (left)). We apply gradient flow on $y$ and $\xi$ with respect to the $L^{2}$-loss in the Fourier domain, in which case the functional derivative $(\partial/\partial\tilde{\mathbf{G}}(is))\mathcal{L}$ simply reduces to the residual:

$$ $\frac{\partial\mathcal{L}}{\partial\tilde{\mathbf{G}}(is)}=-2(\tilde{\mathbf{F }}-\tilde{\mathbf{G}})(is),\qquad\mathcal{L}=\|\tilde{\mathbf{F}}(is)-\tilde{ \mathbf{G}}(is)\|_{L^{2}}.$ $$

In Figure 3 (middle), we show the training dynamics of $(y(\tau),\xi(\tau))$, initialized with different values $(y(0),\xi(0)=3)$, where $\tau$ is the time index of the gradient flow. We make two remarkable observations that corroborate our discussion of the frequency bias during training:

- 1.
Depending on the initialization of $y(0)$, it has two options of moving left or right. Since we fix $\zeta=0$, by 1, a mode $(\hat{y},\hat{\xi})=(-50,5)$ or $(50,0.2)$ in the residual $\tilde{\mathbf{F}}-\tilde{\mathbf{G}}$ impacts the gradient $(\partial/\partial y)\mathcal{L}$ inverse-proportionally to the cube of the distance between the mode $\hat{y}$ and the current $y(\tau)$. Since ${\left|24.5-50\right|^{3}}/{\left|24.5-(-50)\right|^{3}}\approx{0.2}/{5}$, we indeed observe that when $y(0)\leq 24.5$, it tends to move leftward, and rightward otherwise.
- 2.
Although the magnitude of the noises in $[-2\pi,2\pi]$ is only $5\%$ of the smaller mode at $\hat{y}=50$ and $0.2\%$ of the larger mode at $\hat{y}=-50$, once $y(\tau)$ of the trainable LTI system $\tilde{\mathbf{G}}$ enters the noisy region, it gets stuck in a local minimum and never converges to one of the two modes of $\tilde{\mathbf{F}}$ (see Region II in Figure 3). This corroborates our discussion that the training dynamics of $y$ is sensitive to local information and it rarely learns the high frequencies when initialized in the low-frequency region.

Figure: Figure 3: We train an LTI system to learn a noisy bimodal target transfer function. The convergence to a local minimum depends on the initial location of the pole. Left: the ground truth contains a large mode and a small mode, plus some small noises. We want to investigate which mode, if any, our trainable LTI system converges to. Middle: we train the LTI system with respect to the $L^{2}$-loss. We show the trajectories of $(y(\tau),\xi(\tau))$ given different initializations $(y(0),\xi(0)=3)$. The two local minima corresponding to the two modes of $\tilde{\mathbf{F}}$ are shown in red crosses. The green trajectories (initialized in Region I) converge to the mode at $y=-50$, the magenta trajectories (initialized in Region III) converge to the mode at $y=50$, and the black ones (initialized in Region II) converge to neither. Right: the experiment is repeated with the $H^{2}$-loss (see section 5.2).

###### Theorem 1 .

## 5 Tuning the Frequency Bias of an SSM

In section 3 and 4, we analyze the frequency bias of an SSM initialized by HiPPO and trained by a gradient-based algorithm. While we now have a theoretical understanding of the frequency bias, from a practical perspective, we want to be able to tune it. In this section, we design two strategies to enhance, reduce, counterbalance, or even reverse the bias of an SSM against the high frequencies. The two strategies are motivated by our discussion of the initialization (see section 3) and training (see section 4) of an SSM, respectively.

### 5.1 Tuning Frequency Bias by Scaling the Initialization

Since the initialization assigns an SSM some inborn frequency bias, a natural way to tune the frequency bias is to modify the initialization. Here, we introduce a hyperparameter $\alpha>0$ as a simple way to scale the HiPPO initialization defined in 1:

$$ $a_{j}=-0.5+i(-1)^{j}\lfloor j/2\rfloor\alpha\pi,\qquad 1\leq j\leq n.$ (5) $$

Compared to the original HiPPO initialization, we scale the imaginary parts of the eigenvalues of $\mathbf{A}$ by a factor of $\alpha$. By making the modification, we lose the “polynomial projection” interpretation of HiPPO that was originally proposed as a way of explaining the success of the HiPPO initialization; yet, as shown in [, ], this mechanism is no longer regarded as the key for a good initialization. By setting $\alpha<1$, the eigenvalues $a_{j}$ are clustered around the origin, enhancing the bias against the high-frequency modes; conversely, choosing $\alpha>1$ allows us to capture more variations in the high-frequency domain, reducing the frequency bias.

So far, our discussion in the paper is from a perspective of the continuous-time LTI systems acting on continuous time-series. For SSMs, however, the inputs come in a discrete sequence. Hence, we inevitably have to discretize our LTI systems. To study the scaling laws of $\alpha$, we assume in this work that an LTI system is discretized using the bilinear transform [, ] with a sampling interval $\Delta t>0$. Other discretization choices can be similarly studied. Then, given an input sequence $\mathbf{u}\in\mathbb{R}^{L}$ of length $L$, the output $\mathbf{y}$ can be computed by discretizing eq. 2:

$$ $\mathbf{y}=\texttt{iFFT}(\texttt{FFT}(\mathbf{u})\circ\tilde{\mathbf{G}}(\bm{s })),\qquad\bm{s}_{j}=\frac{2}{\Delta t}\frac{\text{exp}\left(i2\pi(j-1)/L \right)-1}{\text{exp}\left(i2\pi(j-1)/L\right)+1},$ (6) $$

with the same transfer function $\tilde{\mathbf{G}}$ in eq. 3. Our goal is to propose general guidelines for an upper bound of $\alpha$. We leave most technical details to Appendix C; but we intuitively explain why $\alpha$ cannot be arbitrarily large for discrete inputs. Given a fixed sampling interval $\Delta t$, there is an upper bound for the frequency, called the Nyquist frequency, above which a signal cannot be “seen” by sampling, causing the so-called aliasing errors [, , ]. As a straightforward example, one cannot distinguish between $\cos(5t)$ and $\cos(t)$ from their samples at $t=k\pi$, $k\in\mathbb{Z}$. Our next result tells us how to avoid aliasing by constraining the range of $\alpha$.

###### Proposition 1 .

Let $a=x+iy$ be given and define ${\mathbf{G}}(is)=1/(is-a)$. Let $\bm{g}={\mathbf{G}}(\bm{s})$ be the vector of length $L$, where $\bm{s}$ is defined in eq. 6. Then, there exist constants $C_{1},C_{2}>0$ such that

$$ $\displaystyle C_{1}\|\bm{g}\|_{2}\!\leq\!\left|y\right|\!\Delta t\!\leq\!C_{2} \|\bm{g}\|_{2},\quad C_{1}\|\bm{g}\|_{\infty}\!\leq\!\frac{1}{1\!+\!\left| \left|y\right|\!-\!2(\Delta t)^{-1}\tan\left(\left(1\!-\!(L\!-\!1)/L\right)\pi /2\right)\right|}\!\leq\!C_{2}\|\bm{g}\|_{\infty}.$ $$

You may have noticed that in 1, we study the norm of the complex $\mathbf{G}$ instead of its real part restriction $\tilde{\mathbf{G}}$. The reason is that in an LTI system parameterized by complex numbers, we multiply $\mathbf{G}$ by a complex number $\xi+i\zeta$ and then extract its real part. Hence, both $\text{Re}(\mathbf{G})$ and $\text{Im}(\mathbf{G})$ are important. By noting that $\max_{1\leq j\leq n}|y_{j}|\approx n\pi/2$ in our scaled initialization in eq. 5, 1 gives us two scaling laws of $\alpha$ that prevent $\|\bm{g}\|_{2}$ and $\|\bm{g}\|_{\infty}$ from vanishing, respectively.
First, the $2$-norm of $\bm{g}$ measures the average-case contribution of the partial fraction $1/(is-a)$ to the input-to-output mapping in eq. 6.
{adjustwidth}15pt15pt
Rule I: (Law of Non-vanishing Average Information) For a fixed task, as $n$ and $\Delta t$ vary, one should scale $\alpha=\mathcal{O}(1/(n\Delta t))$ to preserve the LTI system’s impact on an average input.
Next, the $\infty$-norm of $\bm{g}$ tells us the maximum extent of the system’s action on any inputs. Therefore, if $\|\bm{g}\|_{\infty}$ is too small, then $1/(is-a)$ can be dropped without seriously affecting the system at all.
{adjustwidth}15pt15pt
Rule II: (Law of Nonzero Information) Regardless of the task, one should never take

$$ $\alpha\gg 4\tan\left(\left(1\!-\!(L\!-\!1)/L\right)\pi/2\right)/n\pi\Delta t$ $$

to avoid a partial fraction that does not contribute to the evaluation of the model.
We reemphasize that our scaling laws provide upper bounds of $\alpha$. Of course, one can always choose $\alpha$ to be much smaller to capture the low frequencies better.

###### Proposition 1 .

### 5.2 Tuning Frequency Bias by a Sobolev Filter

In section 5.1, we see that we can scale the HiPPO initialization to redefine the region in the Fourier domain that can be learned by an LTI system. Here, we introduce another way to tune the frequency bias: by applying a Sobolev-norm-based filter. The two strategies both tune the frequency bias, but by different means: scaling the initialization identifies a new set of frequencies that can be learned by the SSM, whereas the filter in this section introduces weights to different frequencies. Our method is rooted in the Sobolev norm, which extends a general $L^{2}$ norm. Imagine that we approximate a ground-truth transfer function $\tilde{\mathbf{F}}(is)$ using $\tilde{\mathbf{G}}(is)$. We can define the loss to be

$$ $\|\tilde{\mathbf{F}}-\tilde{\mathbf{G}}\|_{H^{\beta}}^{2}:=\int_{-\infty}^{ \infty}(1+|s|)^{2\beta}|\tilde{\mathbf{F}}(is)-\tilde{\mathbf{G}}(is)|^{2}\;ds$ (7) $$

for some hyperparameter $\beta\in\mathbb{R}$. The scaling factor $(1+|s|)^{2\beta}$ naturally reweighs the Fourier domain. When $\beta=0$, eq. 7 reduces to the standard $L^{2}$ loss. The high frequencies become less important when $\beta<0$ and more important when $\beta>0$. Unfortunately, as discussed in section 4, there lacks a notion of the “ground-truth” $\tilde{\mathbf{F}}$ for every single LTI system within an SSM, making eq. 7 uncomputable. To address this issue, instead of using a Sobolev loss function, we apply a Sobolev-norm-based filter to the transfer function $\tilde{\mathbf{G}}$ to redefine the dynamical system:

$$ $\hat{\mathbf{y}}(s)=\tilde{\mathbf{G}}^{(\beta)}(is)\hat{\mathbf{u}}(s),\qquad \tilde{\mathbf{G}}^{(\beta)}(is):=(1+|s|)^{\beta}\tilde{\mathbf{G}}(is).$ (8) $$

This equation can be discretized using the same formula in eq. 6 by replacing $\tilde{\mathbf{G}}$ with $\tilde{\mathbf{G}}^{(\beta)}$.

Equation 8 can be alternatively viewed as applying the filter $(1+|s|)^{\beta}$ to the FFT of the input $\mathbf{u}$, which clearly allows us to reweigh the frequency components. Surprisingly, there is even more beyond this intuition: applying the filter allows us to modify the training dynamics of $y_{j}$!

###### Theorem 2 .

Let $\mathcal{L}(\bm{\Theta})$ be a loss function and $\Gamma=(\mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D})$ be a diagonal LTI system in an SSM defined in section 3. For any $\beta\in\mathbb{R}$, we apply the filter in eq. 8 to $\Gamma$ and let $\tilde{\mathbf{G}}^{(\beta)}$ be the new transfer function. Suppose the functional derivative of $\mathcal{L}(\bm{\Theta})$ with respect to $\tilde{\mathbf{G}}^{(\beta)}(is)$ exists and is denoted by $(\partial/\partial\tilde{\mathbf{G}}^{(\beta)}(is))\mathcal{L}$. Then, if $|{(\partial/\partial\tilde{\mathbf{G}}^{(\beta)}(is))\mathcal{L}}|=\mathcal{O}
(|s|^{p})$ for some $p<1-\beta$, we have

$$ $\frac{\partial\mathcal{L}}{\partial y_{j}}\!=\!\int_{-\infty}^{\infty}\!\frac{ \partial\mathcal{L}}{\partial\tilde{\mathbf{G}}^{(\beta)}(is)}\cdot K^{(\beta) }_{j}(s)ds,\quad K^{(\beta)}_{j}(s)\!:=\!(1\!+\!|s|)^{\beta}\frac{\zeta_{j}(\! (s\!-\!y_{j})^{2}\!-\!x_{j}^{2})\!-\!2\xi_{j}x(s\!-\!y_{j})}{[x_{j}^{2}+(s-y_{ j})^{2}]^{2}},$ (9) $$

for every $1\leq j\leq n$. In particular, we have that $\left|K^{(\beta)}_{j}(s)\right|\!=\!\mathcal{O}\!\left(|\zeta_{j}s^{-2+\beta}|
\!+\!|\xi_{j}s^{-3+\beta}|\right)$ as $|s|\rightarrow\infty$.

Compared to 1, the gradient of $\mathcal{L}$ with respect to $y_{j}$ now depends on the loss at frequency $s$ by a factor of $\mathcal{O}(|s|^{-2+\beta})$. Thus, the effect of our Sobolev-norm-based filter is not only a rescaling of the inputs in the frequency domain, but it also allows better learning the high frequencies:
{adjustwidth}15pt15pt
The higher the $\beta$ is, the more sensitive $y_{j}$ is to the high-frequency loss. Hence, $y_{j}$ is no longer constrained by the “local-frequency” loss and will activately learn the high frequencies.
The decay constraint that $|{(\partial/\partial\tilde{\mathbf{G}}^{(\beta)}(is))\mathcal{L}}|=\mathcal{O}
(|s|^{p})$ for some $p<1-\beta$ is needed to guarantee the convergence of the integral in eq. 9. When it is violated, the theoretical statement breaks, but we could still implement the filter in practice, which is similar to []. In Figure 3 (right), we reproduce the illustrative example introduced in Section 4 using the Sobolev-norm-based filter in eq. 8 with $\beta=2$. This is equivalent to training an ordinary LTI system with respect to the $H^{2}$-loss function defined in eq. 7. We find that in this case, the trajectories of $(y(\tau),\xi(\tau))$ always converge to one of the two modes in $\tilde{\mathbf{F}}$ regardless of the initialization, with more of them converging to the high-frequency global minimum on the left. This verifies our theory, because by setting $\beta=2$, we amplify the contribution of the high-frequency residuals in the computation of $(\partial/\partial y)\mathcal{L}$, pushing a $y$ out of the noisy region between $-2\pi$ and $2\pi$. We leave more illustrative experiments to Appendix D, which show the effect of our tuning filter also when $\beta<0$.

###### Theorem 2 .

## 6 Experiments and Discussions

(I) SSMs as Denoising Sequential Autoencoders. We now provide an example of how our two mechanisms allow us to tune frequency bias. In this example, we train an SSM to denoise an image in the CelebA dataset []. We flatten an image into a sequence of pixels in the row-major order and feed it into an S4D model. We collect the corresponding output sequence and reshape it into an image. Similar to the setting of an autoencoder, our objective is to learn the identity map. To make the task non-trivial, we remove the skip connection $\mathbf{D}$ from the LTI systems. During inference, we add two different types of noises to the input images: horizontal or vertical stripes (see Figure 4). While the two types of noises may be visually similar to each other, since we flatten the images using the row-major order, the horizontal stripes turn into low-frequency noises while the vertical stripes become high-frequency ones (see Figure 8). In Figure 4, we show the outputs of the models trained with different values of $\alpha$ and $\beta$ as defined in section 5.1 and 5.2, respectively.

Figure: Figure 4: The outputs of image-denoising S4D models trained with different configurations.

From Figure 4, we can see that as $\alpha$ and $\beta$ increase, our model learns the high frequencies in the input better; consequently, the high-frequency noises get preserved in the outputs and the low-frequency noises are dampened. This corroborates our intuitions from section 5.1 and 5.2. We can further quantify the “pass rates” of the low and high-frequency noises. That is, we compute the percentage of the low and high-frequency noises that are preserved in the output. We show in Table 2 the ratio between the low-pass rate and the high-pass rate, which decreases as $\alpha$ and $\beta$ increase.

**Table 2: We compute the percentage of the low and high-frequency noises that are preserved in the output of a model trained with a pair of configurations $(\alpha,\beta)$. The table shows the ratio between the low-frequency pass rate and the high-frequency pass rate. The more bluish a cell is, the better our model learns the high frequencies. Circled in red is the S4D default.**
|  | $\beta$ |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
|  | $-1.0$ | $-0.5$ | $0.0$ | $0.5$ | $1.0$ |  |
| $\alpha$ | $0.1$ | 4.463e+07 | 2.409e+06 | 1.198e+05 | 4.613e+03 | 1.738e+02 |
| $1$ | 4.912e+05 | 2.124e+05 | 1.758e+04 | 9.595e+02 | 5.730e+01 |  |
| $10$ | 9.654e+04 | 7.465e+03 | 6.073e+02 | 5.699e+01 | 6.394e+00 |  |
| $100$ | 3.243e+00 | 3.745e-02 | 3.801e-03 | 7.299e-05 | 5.963e-06 |  |

(II) Tuning Frequency Bias in the Long-Range Arena. The two tuning strategies section 5.1 and 5.2 are not only good when one needs to deal with a particular high or low frequency, but they also improve the performance of an SSM on general long-range tasks. In Table 3, we show that equipped with the two tuning strategies, our SSM achieves state-of-the-art performance on the Long-Range Arena (LRA) tasks [].

**Table 3: Test accuracies in the Long-Range Arena of different variants of SSMs. 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. The row labeled “Ours” stands for the S4D model equipped with our two tuning strategies. Experiments were run with $5$ random seeds and the medians and the standard deviations are reported. The S4 and S4D results are from the original papers [, ]. The sizes of our models are the same or smaller than the corresponding S4D models.**
| 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$ | $\underline{89.60}$ | $90.00$ | - | $\underline{95.60}$ | $90.10$ | - |
| Liquid S4 [] | $\bm{62.75}$ | $89.02$ | $91.20$ | $\underline{89.50}$ | $94.80$ | $96.66$ | $87.32$ |
| S5 [] | $62.15$ | $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$ |
| Ours | $\bm{62.75}$ | $\bm{89.76}$ | $\bm{92.45}$ | $\bm{90.89}$ | $\bm{95.89}$ | $\underline{97.84}$ | $\bm{88.26}$ |
|  | $\qquad\pm 0.78$ | $\qquad\pm 0.22$ | $\qquad\pm 0.16$ | $\qquad\pm 0.35$ | $\qquad\pm 0.13$ | $\qquad\pm 0.21$ |  |

(III) Ablation Studies. To show the effectiveness of our two tuning mechanisms, we perform ablation studies by training a smaller S4D model to learn the grayscale sCIFAR-10 task. From Figure 5, we see that we obtain better performance when we slightly increase $\alpha$ or decrease $\beta$. One might feel it as a contradiction because increasing $\alpha$ helps the high frequencies to be better learned while decreasing $\beta$ downplays their role. It is not: $\alpha$ and $\beta$ control two different notions of frequency bias: scaling $\alpha$ affects “which frequencies we can learn;” scaling $\beta$ affects “how much we want to learn a certain frequency.” They can be used collaboratively and interactively to attain the optimal extent of frequency bias that we need for a problem.

Figure: Figure 5: Two ablation studies of the tuning strategies proposed in this paper. We train an S4D model with varying parameters of $\alpha$ and $\beta$, respectively. On the left, we see that holding $\beta=0$ (the default value), the model achieves its best performance when $\alpha=4$; on the right, when we fix $\alpha=1$ (the default value), the model performs the best when $\beta=-0.5$.

(IV) More Experiments. One can find more supplementary experiments in Appendix F on wave prediction (see Figure 1) and video generation.

## 7 Conclusion

We formulated the frequency bias of an SSM and showed its existence by analyzing both the initialization and training. We proposed two different tuning mechanisms based on scaling the initialization and on applying a Sobolev-norm-based filter to the transfer function. As a future direction, one could develop ways to analyze the spectral information of the inputs of a problem and use it to guide the selection of the hyperparameters in our tuning mechanisms.

#### Acknowledgments

AY was supported by the SciAI Center, funded by the Office of Naval Research under Grant Number N00014-23-1-2729. SHL would like to acknowledge support from the Wallenberg Initiative on Networks and Quantum Information (WINQ) and the Swedish Research Council (VR/2021-03648). MWM would like to acknowledge LBL’s LDRD initiative for providing partial support of this work. 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.

## Appendix A Proofs

In this section, we provide the proofs of all theoretical statements in the manuscript.

First, we prove the statements about the initialization of the LTI systems. The total variation can be bounded straightforwardly using the decay of the transfer functions.

###### Proof of 1 .

Since $\tilde{\mathbf{G}}$ is the real part of $\mathbf{G}$, its total variation is always no larger than the total variation of $\mathbf{G}$. Hence, we have

$$ $V_{-\infty}^{-B}(\tilde{\mathbf{G}})\leq V_{-\infty}^{-B}(\mathbf{G})\leq\sum_ {j=1}^{n}\left|\frac{c_{j}}{-iB-a_{j}}\right|\leq\sum_{j=1}^{n}\left|\frac{c_{ j}}{B+y_{j}}\right|.$ $$

The other bound is similarly obtained.
∎

###### Proof of 1 .

By 1 and the Hölder’s inequality, we have

$$ $V_{-\infty}^{-B}(\tilde{\mathbf{G}})\leq\sum_{j=1}^{n}\frac{\left|c_{j}\right| }{\left|y_{j}+B\right|}\leq\|\mathbf{B}\circ\mathbf{C}^{\top}\|_{2}\|\mathbf{y }\|_{2},$ $$

where $\mathbf{y}_{j}=1/(y_{j}+B)$. Since $\xi_{j},\zeta_{j}\sim\mathcal{N}(0,1)$ i.i.d., we have that $\|\mathbf{B}\circ\mathbf{C}^{\top}\|_{2}^{2}$ follows the $\chi^{2}$-distribution with degree $2n$. By [], we have with probability at least $1-\delta$ that

$$ $\|\mathbf{B}\circ\mathbf{C}^{\top}\|_{2}^{2}\leq 2n+\sqrt{2n\ln(1/\delta)}+2 \ln(1/\delta)\leq 2(\sqrt{n}+\sqrt{\ln(1/\delta)})^{2}.$ $$

By the definition of $B$, we have that

$$ $\|\mathbf{y}\|_{2}^{2}\leq\frac{n}{(B-n/2)^{2}}.$ $$

The result for $V_{-\infty}^{-B}(\tilde{\mathbf{G}})$ follows from the last two inequalities. The bound on $V_{B}^{\infty}(\tilde{\mathbf{G}})$ can be derived similarly.
∎

We skip the proof of 1 and defer it to Appendix C when we present a detailed derivation of the scaling laws. We next prove the statement about the training dynamics of the imaginary parts of $\text{diag}(\mathbf{A})$.

###### Proof of 1 .

Fix some $1\leq j\leq n$, we first view the transfer function $\tilde{\mathbf{G}}(s,y_{j})$ as a function of two variables. We compute the derivative of $\tilde{\mathbf{G}}(s,y_{j})$ with respect to $y_{j}$:

$$ $\displaystyle\frac{\partial\tilde{\mathbf{G}}(s,y_{j})}{\partial y_{j}}$ $\displaystyle=\frac{\partial}{\partial y_{j}}\left(\sum_{j=1}^{n}\frac{\zeta_{ j}(s-y_{j})-\xi_{j}x_{j}}{x_{j}^{2}+(s-y_{j})^{2}}\right)+\mathbf{D}=\frac{ \partial}{\partial y_{j}}\left(\frac{\zeta_{j}(s-y_{j})-\xi_{j}x_{j}}{x_{j}^{2 }+(s-y_{j})^{2}}\right)$ $\displaystyle=\frac{(x_{j}^{2}+(s-y_{j})^{2})\cdot(\partial/\partial y_{j})( \zeta_{j}(s-y_{j})-\xi_{j}x_{j})}{(x_{j}^{2}+(s-y_{j})^{2})^{2}}$ $\displaystyle\qquad-\frac{(\zeta_{j}(s-y_{j})-\xi_{j}x_{j})\cdot(\partial/ \partial y_{j})(x_{j}^{2}+(s-y_{j})^{2})}{(x_{j}^{2}+(s-y_{j})^{2})^{2}}$ $\displaystyle=\frac{-(x_{j}^{2}+(s-y_{j})^{2})\zeta_{j}}{(x_{j}^{2}+(s-y_{j})^ {2})^{2}}+\frac{2(\zeta_{j}(s-y_{j})-\xi_{j}x_{j})(s-y_{j})}{(x_{j}^{2}+(s-y_{ j})^{2})^{2}}$ $\displaystyle=\frac{\zeta_{j}(-x_{j}^{2}-(s-y_{j})^{2}+2(s-y_{j})^{2})-\xi_{j} (2x_{j}(s-y_{j}))}{(x_{j}^{2}+(s-y_{j})^{2})^{2}}=K_{j}(s).$ $$

Since we assume that $|(\partial/\partial\tilde{\mathbf{G}}(is))\mathcal{L}|=\mathcal{O}(|s|^{p})$ for some $p<1$, the integral in eq. 4 converges. By [, Appendix A], eq. 4 holds.
∎

The statement about the training dynamics of $y_{j}$ given a Sobolev filter follows immediately from 1.

###### Proof of 2 .

Since we assume that $|(\partial/\partial\tilde{\mathbf{G}}(is))\mathcal{L}|=\mathcal{O}(|s|^{p})$ for some $p<1-\beta$, the integral in eq. 9 converges. The result follows by noting that

$$ $\frac{\partial\tilde{\mathbf{G}}^{(\beta)}(s,y_{j})}{\partial y_{j}}=(1+|s|)^{ \beta}\frac{\partial\tilde{\mathbf{G}}(s,y_{j})}{\partial y_{j}}=(1+|s|)^{ \beta}K_{j}(s)=K_{j}^{(\beta)}(s)$ $$

and applying the equation in [, Appendix A].
∎

###### Proof of 1 .

###### Proof of 1 .

###### Proof of 1 .

###### Proof of 2 .

## Appendix B Functional Derivatives

In this section, we briefly introduce the notion of functional derivatives (see Appendix A.1 in [] for a more technical overview). To make our discussion concrete, we do it in the context of 1. Consider the transfer function $\tilde{\mathbf{G}}(is)$ defined in eq. 3. It depends on the model parameters $x_{j},y_{j},\xi_{j}$, and $\zeta_{j}$. In this section, we separate out a single $y_{j}$ for a fixed $1\leq j\leq n$, leaving the remaining parameters unchanged. Then, for every $y\in\mathbb{R}$, we can define $f^{(y)}(s)$ to be the transfer function $\tilde{\mathbf{G}}(is)$ when $y_{j}=y$. Under this setting, the set of all possible transfer functions indexed by $y$, i.e., $\mathcal{F}=\{f^{(y)}|y\in\mathbb{R}\}$ is a subset of a Banach space, say $L^{2}(\mathbb{R})$. To avoid potential confusions, we shall remark that $f^{(y)}$ is not linear in its index $y$, i.e., $f^{(y_{1}+y_{2})}\neq f^{(y_{1})}+f^{(y_{2})}$ in general, neither is $\mathcal{F}$ a subspace of $L^{2}(\mathbb{R})$. This does not impact our following discussion.

Now, consider the loss function $\mathcal{L}$. Given a choice of $y_{j}=y$ and a corresponding transfer function $f^{(y)}$, the loss function maps $f^{(y)}$ to a real number that corresponds on the current loss. Hence, $\mathcal{L}$ can be viewed as a (not necessarily linear) functional of $f^{(y)}$. We would like to ask: how does $\mathcal{L}$ respond to a small change of $f^{(y)}(s)$ at some $s\in\mathbb{R}$? Ideally, this can be measured as

$$ $\lim_{\epsilon\rightarrow 0}\frac{\mathcal{L}(f^{(y)}+\epsilon\delta_{s})-f^{( y)}}{\epsilon},$ (10) $$

where $\delta_{s}$ is the Dirac delta function at $s$. However, the loss function $\mathcal{L}$ is not defined for distributions, making eq. 10 not directly well-defined. To fix this issue, we have to go through the functional derivatives. The idea, as usual in functional analysis, is to pass the difficulty of handling a distribution to smooth functions that approximate it. If there exists a function $(\partial/\partial f^{(y)})\mathcal{L}$ such that the equation

$$ $\int_{-\infty}^{\infty}\frac{\partial\mathcal{L}}{\partial f^{(y)}}(s)\phi(s) \;ds=\lim_{\epsilon\rightarrow 0}\frac{\mathcal{L}(f^{(y)}+\epsilon\phi)- \mathcal{L}(f^{(y)})}{\epsilon}$ $$

holds for all smooth $C^{\infty}_{0}$ functions $\phi$ that are infinitely differentiable and vanish at infinity, then $(\partial/\partial f^{(y)})\mathcal{L}$ is defined to be the functional derivative of $\mathcal{L}$ at $f^{(y)}$. Taking $\{\phi_{j}\}_{j=1}^{\infty}$ to be an approximate identity centered at $s$, we recover eq. 10 using $(\partial/\partial f^{(y)})\mathcal{L}(s)$.

One nice thing about the functional derivatives is that they allow us to write down a continuous analog of the chain rule, which is the meat of 1. To get some intuition, let us first consider a function $\tilde{\mathcal{L}}(y)=\tilde{\mathcal{L}}(f_{1}(y),\ldots,f_{k}(y))$ that depends on $y$ via $k$ intermediate variables $f_{1},\ldots,f_{k}$. Assuming sufficient smoothness conditions, the derivative of $\tilde{\mathcal{L}}$ with respect to $y$ can be calculated using the standard chain rule:

$$ $\frac{\partial\tilde{\mathcal{L}}}{\partial y}=\sum_{j=1}^{k}\frac{\partial \tilde{\mathcal{L}}}{\partial f_{j}}\frac{\partial f_{j}}{\partial y}=\begin{ bmatrix}\frac{\partial\tilde{\mathcal{L}}}{\partial f_{1}}&\cdots&\frac{ \partial\tilde{\mathcal{L}}}{\partial f_{k}}\end{bmatrix}\begin{bmatrix}\frac{ \partial f_{1}}{\partial y}\\ \vdots\\ \frac{\partial f_{k}}{\partial y}\end{bmatrix}.$ (11) $$

The only difference in the case of $\mathcal{L}$ is that instead of depending on $k$ discrete intermediate variables $f_{1},\ldots,f_{k}$, our $\mathcal{L}$ depends on a continuous family of intermediate variables $f^{(y)}(s)$ indexed by $s\in\mathbb{R}$. In this case, one would naturally expect that in eq. 11, the sum becomes an integral, or equivalently, the row and the column vectors become the row and the column functions. This is indeed the case given our functional derivative:

$$ $\frac{\partial{\mathcal{L}}}{\partial y}=\int_{-\infty}^{\infty}\frac{\partial {\mathcal{L}}}{\partial f^{(y)}}(s)\frac{\partial f^{(y)}(s)}{\partial y}\;ds.$ $$

This formula can be found in [, (A.24)] and [, sect. 2.3].

## Appendix C Scaling Laws of the Initialization

In this section, we expand our discussions in section 5.1 and give the proof of 1. Throughout this section, we assume that we use the bilinear transform to discretize our continuous-time LTI system. The length of our sequence is $L$ and the sampling interval is $\Delta t$. The bilinear transform is essentially a Mobius transform between the closed left half-plane of the $s$-domain and the closed unit disk in the $z$-domain. Hence, it gives us two ways to study this filter — by either transplanting the transfer function $\tilde{\mathbf{G}}$ onto the unit circle and analyzing in the discrete domain or by transplanting the FFT nodes from the $z$-domain to the imaginary axis in the $s$-domain. The two ways are equivalent, but we choose the second method for simplicity.

The output of an LTI system can be computed by

$$ $\mathbf{y}=\texttt{iFFT}(\texttt{FFT}(\mathbf{u})\circ\overline{\mathbf{G}}( \bm{\omega})),$ $$

where $\overline{\mathbf{G}}$ is the transfer function of the discrete system and where

$$ $\bm{\omega}=\begin{bmatrix}\text{exp}\left(2\pi i\frac{0}{L}\right)&\cdots& \text{exp}\left(2\pi i\frac{L-1}{L}\right)\end{bmatrix}^{\top}$ $$

is the length-$L$ vector consisting of $L$th roots of unity. We do not have direct access to $\overline{\mathbf{G}}$, but we do know $\tilde{\mathbf{G}}$, its continuous analog, in the partial fractions format. They are related by the following equation:

$$ $\overline{\mathbf{G}}(z)=\tilde{\mathbf{G}}(s),\qquad s=\frac{2}{\Delta t} \frac{z-1}{z+1}.$ $$

In that case, the vector $\overline{G}(\bm{\omega})$ can be equivalently written as

$$ $\overline{\mathbf{G}}(\bm{\omega}_{j})=\tilde{\mathbf{G}}\left(\frac{2}{\Delta t }\frac{\text{exp}\left(2\pi i\frac{j}{L}\right)-1}{\text{exp}\left(2\pi i\frac {j}{L}\right)+1}\right)=\tilde{\mathbf{G}}\left(i\frac{2}{\Delta t}\tan\left( \pi\frac{j}{L}\right)\right).$ $$

This is how we obtained eq. 6. The locations of the new samplers on the imaginary axis are shown in Figure 6, with $L=101$ and $\Delta t=0.01$.

Figure: Figure 6: The poles of the FFT samplers in the $z$-domain and the samplers in the $s$-domain for $L=101$ and $\Delta t=0.01$. Due to the essential difficulty of plotting the entire real line on a logarithmic scale, we only show the semilog plot of the upper half-plane of the $s$-domain (right). Hence, on the right figure, we omit the samplers corresponding to the lower arc of the unit circle on the left.
Refer to caption: x2.png

Note that the right figure (the $s$-domain) is on a logarithmic scale and only the upper half-plane is shown due to the scale. We also choose $L$ odd; when $L$ is even, a pole is placed “at infinity” in the $s$-domain, at which any partial fraction vanishes. Why do we go through all the pains to study this bilinear transformation? The reason is that it gives us a guideline for scaling the poles. For instance, for $L=101$ and $\Delta t=0.01$, if a pole has a much larger imaginary part than $10^{4}$, then the discrete sequence will hardly see the effect of this partial fraction even though the underlying continuous system will. This corresponds to the intuition behind the aliasing error that we discussed in the main text.

As in 1, it suffices to study a single partial fraction instead of all. Hence, instead of studying the entire transfer function together, we focus on one component of it:

$$ $\mathbf{G}(is)=\frac{1}{is-a},\qquad\text{Re}(s)=0,\quad\text{Re}(a)<0.$ $$

This is a partial fraction in the $s$-domain. For a fixed $L$ and $\Delta t$, this partial fraction corresponds to a bounded linear operator $\mathcal{G}:\ell^{2}([L])\rightarrow\ell^{2}([L])$ that maps an input sequence to an output sequence, where $[L]=\{1,\ldots,L\}$ is the set of the first $L$ natural numbers. We consider the norm of this operator, where we will find that as $|\text{Im}(s)|\rightarrow\infty$, the norm of the operator vanishes. The rate of vanishing will guide us in selecting an appropriate range for the pole. So, how can we tell the norm of this operator? By definition, the norm of $\mathcal{G}$ is defined by

$$ $\|\mathcal{G}\|_{\ell^{2}\rightarrow\ell^{2}}=\sup_{\|\mathbf{u}\|_{\ell^{2}}= 1}\|\mathbf{y}\|_{\ell^{2}}=\sup_{\|\hat{\mathbf{u}}\|_{\ell^{2}}=1}\|\hat{ \mathbf{y}}\|_{\ell^{2}},$ $$

where $\mathbf{y}$ is the output of the operator $\mathcal{G}$ given input $\mathbf{u}$, i.e., $\mathbf{y}=\mathcal{G}\mathbf{u}$, and the second step follows from the Parseval’s identity. By Hölder’s inequality, we further have

$$ $\|\mathcal{G}\|_{\ell^{2}\rightarrow\ell^{2}}=\sup_{\|\hat{\mathbf{u}}\|_{\ell ^{2}}=1}\|\hat{\mathbf{u}}\circ\overline{\mathbf{G}}(\bm{\omega})\|_{\ell^{2}} \leq\|\bm{g}\|_{\ell^{\infty}},$ $$

where $\overline{\mathbf{G}}(\bm{\omega})=\bm{g}$ is the sample vector of the bilinearly transformed transfer function
of $\mathbf{G}$ in the $z$-domain, i.e.,

$$ $\bm{g}=\overline{\mathbf{G}}(\bm{\omega}),\qquad\overline{\mathbf{G}}(\bm{ \omega})_{j}=\overline{\mathbf{G}}(\bm{\omega}_{j})=\mathbf{G}\left(i\frac{2}{ \Delta t}\tan\left(\pi\frac{j}{L}\right)\right)=\frac{1}{i2\tan\left(\pi j/L \right)/\Delta t-a}.$ $$

Hence, we have

$$ $\left|\bm{g}_{j}\right|^{2}=\frac{1}{\text{Re}(a)^{2}+\left(\text{Im}(a)-2\tan \left(\pi j/L\right)/\Delta t\right)^{2}}.$ $$

When $\text{Im}(a)>2\tan\left(\pi j/L\right)/\Delta t$ for all $j$, $\left|\bm{g}_{j}\right|^{2}$ is maximized when $j=\lfloor L/2\rfloor-1$, in which case we have

$$ $\|\mathcal{G}\|_{\ell^{2}\rightarrow\ell^{2}}^{2}\leq\frac{1}{\text{Re}(a)^{2} +\left(\text{Im}(a)-(2/\Delta t)\tan\left((\pi/2)(1-(L-1)/L)\right)\right)^{2}}.$ (12) $$

This gives us the second rule (Law of Zero Information) when scaling the diagonal of $\mathbf{A}$. This is a worst-case analysis, where we essentially assume that the Fourier coefficient of $\mathbf{u}$ is one-hot at the highest frequency. In practice, of course, this assumption is a bit unrealistic; in fact, the Fourier coefficients usually decay as the frequency gets higher. Therefore, we should derive another rule for the average-case scenario. We consider the operator $\hat{\mathcal{G}}$ that maps the Fourier coefficients of the inputs to those of the outputs, then the norm $\|\hat{\mathcal{G}}\|_{\ell^{\infty}\rightarrow\ell^{2}}$ is a good average-case estimate, because

$$ $\text{arg\,max}_{\|\hat{\mathbf{u}}\|_{\ell^{\infty}}=1}\|\hat{\mathcal{G}} \hat{\mathbf{u}}\|_{\ell^{2}}$ $$

is necessarily at a vertex of the simplex defined by $\|\hat{\mathbf{u}}\|_{\ell^{\infty}}\leq 1$ (^2^22Note that we can use max instead of sup because the domain $\{\|\hat{\mathbf{u}}\|_{\ell^{\infty}}=1\}$ is compact.). That is, $\hat{\mathbf{u}}_{j}=\pm 1$ for all $1\leq j\leq L$. Now, using the Hölder’s inequality again, we have

$$ $\|\hat{\mathcal{G}}\|_{\ell^{\infty}\rightarrow\ell^{2}}=\sup_{\|\hat{\mathbf{ u}}\|_{\ell^{\infty}}=1}\|\hat{\mathbf{u}}\circ\bm{g}\|_{\ell^{2}}\leq\|\bm{g} \|_{\ell^{2}}.$ $$

Hence, instead of studying the $\ell^{\infty}$-norm of $\bm{g}$, we consider the $\ell^{2}$-norm for the average-case estimate. The precise computation of the $\ell^{2}$-norm can be hard, but let us write out the full expression:

$$ $\|\bm{g}\|_{\ell^{2}}^{2}=\sum_{j=1}^{L}|\bm{g}|^{2}\leq\sum_{j=1}^{L}\frac{1} {\text{Re}(a)^{2}+\left(\text{Im}(a)-(2/\Delta t)\tan\left(\pi j/L\right) \right)^{2}}.$ $$

Given the imaginary part of $a>0$, we grab all Fourier nodes on the $j\omega$ axis that are below $\text{Im}(a)/2$ and lower-bound them; we also grab all above $\text{Im}(a)/2$ and assume that they collapse to $\text{Im}(a)$. This gives us an estimate of the $\ell^{2}$ norm:

$$ $\displaystyle\|\bm{g}\|_{\ell^{2}}^{2}$ $\displaystyle\leq\left(\frac{\left|\left\{j\middle|(2/\Delta t)\tan\left(\pi j /L\right)\leq\text{Im}(a)/2\right\}\right|}{\text{Re}(a)^{2}+\text{Im}(a)^{2}/ 4}+\frac{\left|\left\{j\middle|(2/\Delta t)\tan\left(\pi j/L\right)>\text{Im}( a)/2\right\}\right|}{\text{Re}(a)^{2}}\right)$ (13) $\displaystyle\leq\left(\frac{L(2\arctan(\text{Im}(a)\Delta t/4)/\pi)+1}{\text{ Re}(a)^{2}+\text{Im}(a)^{2}/4}+\frac{L(1-2\arctan(\text{Im}(a)\Delta t/4)/\pi) +1}{\text{Re}(a)^{2}}\right)$ $\displaystyle=L\left(\underbrace{\frac{(2\arctan(\text{Im}(a)\Delta t/4)/\pi)+ 1/L}{\text{Re}(a)^{2}+\text{Im}(a)^{2}/4}}_{N_{1}}+\underbrace{\frac{(1-2 \arctan(\text{Im}(a)\Delta t/4)/\pi)+1/L}{\text{Re}(a)^{2}}}_{N_{2}}\right)$ $$

###### Proof of 1 .

Given eq. 12 and (13), 1 follows immediately.
∎

Let us take a closer look at this expression. Ideally, the $\ell^{2}$-norm of $\bm{g}$ should be independent of $L$ and $\Delta t$; that is, as $L\rightarrow\infty$ and $\Delta t\rightarrow 0$, we do not want $\|\bm{g}\|_{\ell^{2}}^{2}/L$ to diminish. First, we note that $N_{1}$ and $N_{2}$ are independent of $L$ as $L\rightarrow\infty$. As $\Delta t\rightarrow 0$, $N_{1}$ inevitably vanish, regardless of the location of $\text{Im}(a)$. In order to maintain $N_{2}$ a constant, we would need $\text{Im}(a)\Delta t/4$ to not blow up. This gives us the first rule (Law of Zero Information) for scaling the poles. We can further work out some constants in $\mathcal{O}$ to be used in practice. For example, to guarantee that $\text{Im}(a)$ is smaller than the top $5\%$ Fourier nodes, we would need that

$$ $\frac{2}{\pi}\arctan\left(\frac{\text{Im}(a)\Delta t}{4}\right)\leq 0.95 \Rightarrow\text{Im}(a)\leq\frac{50.82}{\Delta t}.$ $$

In particular, eq. 12 and eq. 13 together give us the proof of the lower bounds in 1. The upper bounds are proved by noting all all derivations in this section are asymptotically tight.

###### Proof of 1 .

## Appendix D More Numerical Experiments on the Illustrative Example

In section 4 and 5.2, we see that using a Sobolev-norm-based filter with $\beta>0$, one is able to escape from the local minima caused by small local noises. In this section, we present a similar set of experiments to show the effect of our filter, even when $\beta<0$. We choose our new objective function to be

$$ $\tilde{\mathbf{F}}(is)=\text{Re}\left(\frac{5}{is-(-1-75i)}+\frac{0.2}{is-(-1+ 25i)}\right),\qquad s\in\mathbb{R}.$ $$

Compared to the objective in section 4, we see two differences. First, we remove the sinusoidal noises around the origin. Second, we shift the locations of the two modes in the target: instead of locating at $s=-50$ and $s=50$, we shift them to $s=-75$ and $s=25$, respectively. This allows us to have a large high-frequency mode and a small low-frequency mode in the ground truth.

Figure: Figure 7: An SSM is trained with a filter based on the Sobolev-norm (see section 5.2). The plots are read in the same way as those in Figure 1. The transfer function converges to one of the two local minima. As $\beta$ ranges from $-2$ to $2$, the transfer function becomes more sensitive to the high-frequency global information rather than the local information. Hence, the edge between the two different convergences shifts rightward.

We show the results in Figure 7 when we train an LTI system using a Sobolev-norm-based filter with different values of $\beta$. Note that when $\beta=0$, the picture only differs from Figure 1 (middle) in the frequency labels because in that case, the gradient $(\partial/\partial y)\mathcal{L}$ only cares about the relative difference $y-s$ but not the absolute values of $s$. From Figure 7, we see that as $\beta$ increases, more trajectories converge to the local minimum near $(\xi=5,y=-75)$. The reason is that a larger $\beta$ favors a higher frequency (see 2).

## Appendix E Details of the Experiments

### E.1 Denoising Sequential Autoencoder

For every image in the CelebA dataset, we reshaped it to have a resolution of $1024\times 256$ pixels to allow for higher-frequency noises. We trained a single-layer S4D model with $n=128$ and $\texttt{d_model}=3$. We dropped the skip connection $\mathbf{D}$ from the model. The model was trained using the MSE loss. That is, for every predicted sequence of pixels, we compared the model against the true image and computed the $2$-norm of the difference vector.

Figure: Figure 8: When a noisy image is flattened into pixels using the row-major order, the noises induced by horizontal stripes become low-frequency noises while those induced by vertical stripes become high frequencies.

We trained the model on the original images, i.e., those without any noises. When we inferred from the model, we added noises to the inputs, introducing about $10$ cycles of horizontal or vertical stripes, respectively. Our noises were large, almost shielding the underlying images. When the value of a pixel was out of range, then we ignored such as issue during training; we clipped its value to the appropriate range when rendering the image in Figure 4.

To obtain the numbers in Table 2, we computed with our trained models, where we set the inputs to be pure horizontal or vertical noises with no underlying images. Then, we evaluated the size of the output image and took the ratio of the outputs over the inputs. We call this value the “pass rate” of a particular noise. Table 2 shows the ratio between the pass rate of the low-frequency noises over the high-frequency ones. Our model did not have a nonlinear activation function, which made the model linear. Hence, it does not matter what the magnitude of the inputs was.

### E.2 Long-Range Arena

In this section, we present the hyperparameters of our models trained on the Long-Range Arena tasks. Our model architecture and hyperparameters are almost identical to those of the S4D models reported in [], with only two exceptions: for the ListOps experiment, we set $n=2$ instead of $n=64$, which aligns with [] instead; for the PathX experiment, we set $\texttt{d_model}=128$ to reduce the computational burden. We do not report the dropout rates since they are set to be the same as those in []. Also, we made $\beta$ a trainable parameter.

**Table 4: Configurations of our S4D model, where LR, BS, and WD stand for learning rate, batch size, and weight decay, respectively. The hyperparameter $\alpha$ is the scaling factor introduced in section 5.1. We set the parameter in section 5.2 as a trainable parameter to reduce the need for hyperparameter tuning.**
| Task | Depth | #Features | Norm | Prenorm | $\alpha$ | LR | BS | Epochs | WD | $\Delta$ Range |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| ListOps | 8 | 256 | BN | False | 3 | 0.002 | 50 | 80 | 0.05 | (1e-3,1e0) |
| Text | 6 | 256 | BN | True | 5 | 0.01 | 32 | 300 | 0.05 | (1e-3,1e-1) |
| Retrieval | 6 | 128 | BN | True | 3 | 0.004 | 64 | 40 | 0.03 | (1e-3,1e-1) |
| Image | 6 | 512 | LN | False | 3 | 0.01 | 50 | 1000 | 0.01 | (1e-3,1e-1) |
| Pathfinder | 6 | 256 | BN | True | 3 | 0.004 | 64 | 300 | 0.03 | (1e-3,1e-1) |
| Path-X | 6 | 128 | BN | True | 5 | 0.001 | 20 | 80 | 0.03 | (1e-4,1e-1) |

## Appendix F Supplementary Experiments

### F.1 Predict the Magnitudes of Waves

In Figure 1, we see an example of the frequency bias of SSMs, where the model is better at extracting the wave information of a low-frequency wave than a high-frequency one. In this section, we produce more examples on the same task to show that one is able to tune frequency bias by playing with $\alpha$ and $\beta$ we introduced in section 5.1 and 5.2, respectively.

Figure: Figure 9: Reproduction of the experiment shown in Figure 1, with difference choices of $(\alpha,\beta)$. From the left to right, our choices of $(\alpha,\beta)$ are $(0.01,-1)$, $(1,0)$ (the default), $(10,0.5)$, and $(100,1)$, respectively. The legend applies to all pictures. We can see that these choices allow us to enhance, counterbalance, or even reverse the frequency bias.

We see from Figure 9 that by tuning hyperparameters, we can change the frequency bias of an SSM. In particular, when $\alpha=100$ and $\beta=1$, we reversed the frequency bias so that the magnitude of the low-frequency wave $\cos(t)$ cannot be well-predicted, while a high-frequency wave is captured relatively well.

### F.2 Tuning Frequency Bias in Moving MNIST Video Prediction

In this section, we present an experiment to tune frequency bias in a video prediction task. We show that our frequency bias analysis and the tuning strategies not only work for vanilla SSMs but also their variants. We examine a model architecture called ConvS5 that combines SSMs and spatial convolution []. We apply the model to predict movies from the Moving MNIST dataset []. In this dataset, two (or more) digits taken from the MNIST dataset [] move on a larger canvas and bounce when touching the border. This forms a video over time. In our experiment, we slightly modify the movies by coloring the two digits. In particular, every movie contains two moving digits — a fast-moving red one and a slow-moving blue one. The speed of the red digit is ten times that of the blue digit; consequently, the red digit can be considered as a “high-frequency” component, whereas the blue digit is a “low-frequency” component. Our goal in this experiment is to use a ConvS5 model to generate up to $100$ frames, conditioned on $500$ frames. The ConvS5 model applies LTI systems to the time domain (i.e., the axis of the frames), but in the meantime incorporates spatial convolutions in the LTI systems, where the LTI systems are still initialized by the HiPPO initialization.

In this experiment, we train two models using two different initializations. The first initialization we use is the default HiPPO initialization. Then, we try another initialization, where for every $y_{j}$ that is the imaginary part of an eigenvalue of $\mathbf{A}$, we transform $y_{j}$ by

$$ $y_{j}\mapsto\text{sign}(y_{j})(\left|y_{j}\right|+200).$ (14) $$

That is, we shift every $y_{j}$ away from the origin by $200$. This does not correspond to any $\alpha>0$ that we introduced in section 5.1, but our intuition is still based on our discussions in section 3 and 4. That is, when we move away every $y_{j}$ that is contained in $[-200,200]$, our model is incapable of handling the low frequencies. This is indeed observed in Figure 10: when we use the original HiPPO initialization, the high-frequency red digit cannot be predicted, whereas when we modify the initialization based on eq. 14, we well-predicted the red digit but the low-frequency blue digit is completely distorted.

Figure: Figure 10: A ConvS5 model is trained to predict the Moving MNIST videos. In “Model1”, we use the original HiPPO initialization; in “Model2”, we modify the initialization based on eq. 14. We see that when we use the HiPPO initialization, only the slow-moving blue digit can be generated; on the other hand, pushing all eigenvalues of $\mathbf{A}$ to the high-frequency regions (see eq. 14) allows us to predict the fast-moving red digit.

## Peer Reviews

confidence: 3
contribution: 4
initial_rating: 8
presentation: 4
questions: 1. The authors may need to define HiPPO in line 92 Page 2 before mentioning it. 

2. How will these discoveries help for general neural network settings? 

3. Would the Sobolev filter difficult to compute in practice (when doing the SSM training)? 

4. How to interpret Figure 4 when having low-freq. noise and high-freq. noise. What is the difference? Additionally, why low-freq. noise can be denoised (with higher-frequency-pass filter) but the images in the bottom rows will be reconstructed?
soundness: 4
strengths: This paper presents a remarkable analysis of the frequency bias in state-space models (SSMs). It has been commonly observed that neural networks tend to fit low-frequency information first, often filtering out high-frequency data in many cases. This paper critically emphasizes the importance of initialization and introduces a strategy to mitigate this problem using a Sobolev filter. The figures presented in the experiments are highly inspiring. As a theoretical work, this paper provides sufficient experimental results, demonstrating strong performance.
summary: This paper studies the frequency bias for State-space-models (SSMs). With strategies to select better initialization and Sobolev filtering (during training), this problem could be mitigated. In the experimental study, this manuscript also provides extensive results to show the effectiveness of their proposed strategies.
weaknesses: The weaknesses of this paper mostly come from two parts. 

1. The analysis is based on SSM which is potentially hard to be generalized on general neural networks. 

2. Even though the reason is getting clearer, these strategies of this paper could still be suboptimal for practitioners.

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: For the numerical experiments (Table 3) on the long range arena benchmark, how is the $A$ matrix initialized? Is $A$ initialized by the method in Corollary 1? There are some ablation studies on the scaling of the imaginary part of $A$ in the S4D paper [1]. It is shown in [1] that scaling all imaginary parts by a constant factor substantially reduce the performance of S4D, so I am curious that if we only scaling the imaginary part of $a$, does it help to improve the performance on the long range arena benchmark? It would be more convincing if the authors can provide more experiment results of only scale $a$, only train $\beta$, the results in Table 3 is the case for scale $a$ + train $\beta$.






[1] Gu, Albert, et al. "On the parameterization and initialization of diagonal state space models." Advances in Neural Information Processing Systems 35 (2022): 35971-35983.
soundness: 3
strengths: 1. The organization and presentation of this paper is smooth and clear, and it provides a better understanding on the training mechanisms of state space models on sequence modeling.

2. I find the paper to be well-written and easy to follow. The overall topic of state space models is important.
summary: In this paper, the authors find that when training state space models with canonical initialization methods, there is an implicit bias such that the models are implicitly toward to capture low-frequency components instead of high-frequency ones.
To encourage the models to capture a more broader frequency pattern, the authors propose two mechanisms: one is scaling the initialization and another method is to apply a Sobolev-norm-based filter on the state space models.
By tuning the frequency bias, the state space models are shown to be able to have better performance on the long range arena benchmark.
weaknesses: 1. In the end of Section 2, the authors state that the larger the total variation of the transfer function is, the better an LTI system is at distinguishing the Fourier modes. This claim is not intuitive for me, and the total variation cannot distinguish the following two different cases: the first one is large amplitude with low frequency and the second one is small amplitude with high frequency. For example, $G$ is an impulse response vs $G$ is a sinusoidal wave with a small amplitude. From my understanding, these two LTI systems have different ability on distinguishing the Fourier modes.

2. The statement after Lemma 1 that small $|y_i|$ induces small total variation seems to be wrong. From the upper bounds in Lemma 1, if $|y_i|$ decreases, the the upper bounds increase, which means that the total variation will be large.

3. The initialization method in Corollary1 is not a commonly used method. For S4D, the initialization methods for $a$ are mainly S4D-Legs and S4D-Lin. Why not choose these two initialization methods instead?

---

confidence: 4
contribution: 3
initial_rating: 8
presentation: 4
questions: 1. Can the authors explain more clearly why the rate of change of the transfer function is a proper measure of frequency bias? Traditionally (e.g. low pass filters) would consider the support of the transfer function, or its magnitude over frequency intervals, as a measure of its response to frequency components. Some proper definition of "frequency sensitivity" would make the definition of frequency bias in terms of TV more transparent.
2. After proposition 1, Rule II: what happens in the limit $L\to\infty$?
3. Eq 7: I find it strange that the weighting factor is $(1+|s|)^{2\beta}$ instead of $(1+|s|^2)^{\beta}$, as in the usual Bessel potential spaces. What's the motivation for this departure? In any case, the name "the Sobolev norm" is unclear since there are many Sobolev spaces. 
4. For the mitigation algorithms, what are some general tuning strategies for $\alpha$ and $\beta$? Fig 5 seems to say $\alpha$ is more sensitive than $\beta$. Some comments on this would be helpful.
soundness: 3
strengths: The paper is very clearly written, and the motivation and implications of the theoretical results are clearly explained. The identification of the frequency bias is useful for problems where high frequency components need to be extracted. The proposed methods to mitigate the frequency bias appear simple to implement, and are principled.
summary: This paper studies the frequency bias of state space models, where the authors identify a tendency, under generic initialisations of the state space matrix, to favor the response to low frequency signals. The authors give a definition of the frequency bias through the total variation of the transfer function, in that a low TV over an interval in frequency space is understood as a bias to not capture this region effectively. The authors then show that generic initialisations indeed gives lower TV for high frequency regions, and moreover that this cannot be fixed by optimisation (in the sense that the gradient on the frequency component is small if there is little initial overlap). 

Two solutions are proposed to mitigate this issue: better initialisation where a hyper-parameter is introduced to weigh higher frequencies differently, or changing the loss function by adding a weighting function in frequency space to promote training larger frequencies (or vice versa). Improvements are observed on both a synthetic denoising task and benchmarks from long range arena.
weaknesses: 1. The paper can be improved by clarifying some notational and definition issues, and these are detailed in the questions below.
2. Notational suggestion: the notation of using $y_j$ to denote imaginary part of the diagonal of $A$ can be changed to avoid confusion with the output of the network, also called $y$.
3. Since the SSM transfer function is differentiable, it may be simpler to define the total variation as the integral of its derivatives to avoid unncessary complications.

---

confidence: 2
contribution: 3
initial_rating: 8
presentation: 3
questions: Could you provide some guidelines for tuning the α and β parameters? How much tuning was necessary to achieve the results in Table 3?

How could the analysis be extended to non-diagonal systems?
soundness: 4
strengths: Well-written with effective illustrative explanations.

Provides a theoretically sound analysis of a complex problem.

Practical relevance is high.
summary: This paper addresses the frequency bias of State Space Models (SSMs) and provides methods to adjust this bias. The authors introduce a formal framework to characterize the frequency behavior of SSMs and reveal that standard SSMs, as commonly described in recent literature, tend to learn low frequencies. While this low-frequency bias is often advantageous, the authors argue that it is possible to improve the model's frequency response.

To achieve this, they propose two different approaches: the first modifies initialization to control which frequencies the model can learn and the second alters the transfer function to adjust how much the model emphasizes particular frequencies.

Although I am not a expert in SSMs, I found the paper well-structured and effective in explaining complex issues with clear, illustrative examples.
weaknesses: The analysis is limited to diagonal systems.

Tuning the additional hyperparameters, α and β, may pose challenges in practice.

Meta-review: The authors are concerned with the observation that state space models contain an implicit bias for ignoring high-frequency components. After conducting a theoretical analysis to shed light into this phenomenon, they formulate a theory (related to initialization) which, in turn, gives rise to their two proposed solutions. This is, therefore, an overall well-conducted piece of research, from motivation to solution proposal, which is a great addition to the literature.

Decision: Accept (Spotlight)