Year: 2021

Abstract:
A central goal of sequence modeling is designing a single principled model that can address sequence data across a range of modalities and tasks, particularly on long-range dependencies. Although conventional models including RNNs, CNNs, and Transformers have specialized variants for capturing long dependencies, they still struggle to scale to very long sequences of 10000 or more steps. A promising recent approach proposed modeling sequences by simulating the fundamental state space model (SSM) x′(t) = Ax(t) +Bu(t), y(t) = Cx(t) +Du(t), and showed that for appropriate choices of the state matrix A, this system could handle longrange dependencies mathematically and empirically. However, this method has prohibitive computation and memory requirements, rendering it infeasible as a general sequence modeling solution. We propose the Structured State Space (S4) sequence model based on a new parameterization for the SSM, and show that it can be computed much more efficiently than prior approaches while preserving their theoretical strengths. Our technique involves conditioning A with a low-rank correction, allowing it to be diagonalized stably and reducing the SSM to the well-studied computation of a Cauchy kernel. S4 achieves strong empirical results across a diverse range of established benchmarks, including (i) 91% accuracy on sequential CIFAR-10 with no data augmentation or auxiliary losses, on par with a larger 2-D ResNet, (ii) substantially closing the gap to Transformers on image and language modeling tasks, while performing generation 60× faster (iii) SoTA on every task from the Long Range Arena benchmark, including solving the challenging Path-X task of length 16k that all prior work fails on, while being as efficient as all competitors.

heading: 1 INTRODUCTION
text: A central problem in sequence modeling is efficiently handling data that contains long-range dependencies (LRDs). Real-world time-series data often requires reasoning over tens of thousands of time steps, while few sequence models address even thousands of time steps. For instance, results from the long-range arena (LRA) benchmark (Tay et al., 2021) highlight that sequence models today perform poorly on LRD tasks, including one (Path-X) where no model performs better than random guessing.
Since LRDs are perhaps the foremost challenge for sequence models, all standard model families such as continuous-time models (CTMs), RNNs, CNNs, and Transformers include many specialized variants designed to address them. Modern examples include orthogonal and Lipschitz RNNs (Arjovsky et al., 2016; Erichson et al., 2021) to combat vanishing gradients, dilated convolutions to increase context size (Bai et al., 2018; Oord et al., 2016), and an increasingly vast family of efficient Transformers that reduce the quadratic dependence on sequence length (Katharopoulos et al., 2020; Choromanski et al., 2020). Despite being designed for LRDs, these solutions still perform poorly on challenging benchmarks such as LRA (Tay et al., 2021) or raw audio classification (Gu et al., 2021).
An alternative approach to LRDs was recently introduced based on the state space model (SSM) (Fig. 1). SSMs are a foundational scientific model used in fields such as control theory, computational neuroscience, and many more, but have not been applicable to deep learning for concrete theoretical reasons. In particular, Gu et al. (2021) showed that deep SSMs actually struggle even on simple tasks, but can perform exceptionally well when equipped with special state matrices A recently derived to solve a problem of continuous-time memorization (Voelker et al., 2019; Gu et al., 2020a). Their
Linear State Space Layer (LSSL) conceptually unifies the strengths of CTM, RNN and CNN models, and provides a proof of concept that deep SSMs can address LRDs in principle.
Unfortunately, the LSSL is infeasible to use in practice because of prohibitive computation and memory requirements induced by the state representation. For state dimension N and sequence length L, computing the latent state requires O(N2L) operations and O(NL) space – compared to a Ω(L + N) lower bound for both. Thus for reasonably sized models (e.g. N = 256 in Gu et al. (2021)), the LSSL uses orders of magnitude more memory than comparably-sized RNNs or CNNs. Although theoretically efficient algorithms for the LSSL were proposed, we show that these are numerically unstable. In particular, the special A matrix is highly non-normal in the linear algebraic sense, which prevents the application of conventional algorithmic techniques. Consequently, although the LSSL showed that SSMs have strong performance, they are currently computationally impractical as a general sequence modeling solution.
In this work, we introduce the Structured State Space (S4) sequence model based on the SSM that solves the critical computational bottleneck in previous work. Technically, S4 reparameterizes the structured state matrices A appearing in Voelker et al. (2019); Gu et al. (2020a) by decomposing them as the sum of a low-rank and skew-symmetric term. Additionally, instead of expanding the standard SSM in coefficient space, we compute its truncated generating function in frequency space, which can be simplified into a multipole-like evaluation. Combining these two ideas, we show that the low-rank term can be corrected by the Woodbury identity while the skew-symmetric term can be diagonalized stably, ultimately reducing to a well-studied and theoretically stable Cauchy kernel (Pan, 2001; 2017). This results in Õ(N + L) computation and O(N + L) memory usage, which is essentially tight for sequence models. Compared to the LSSL, S4 is up to 30× faster with 400× less memory usage, while exceeding the LSSL’s performance empirically.
Empirically, S4 significantly advances the state-of-the-art for LRD. On the LRA benchmark for efficient sequence models, S4 is as fast as all baselines while outperforming them by 20+ points on average. S4 is the first model to solve the difficult LRA Path-X task (length-16384), achieving 88% accuracy compared to 50% random guessing for all prior work. On speech classification with length-16000 sequences, S4 halves the test error (1.7%) of specialized Speech CNNs – by contrast, all RNN and Transformer baselines fail to learn (≥ 70% error).
Towards a general-purpose sequence model. Beyond LRD, a broad goal of machine learning is to develop a single model that can be used across a wide range of problems. Models today are typically specialized to solve problems from a particular domain (e.g. images, audio, text, time-series), and enable a narrow range of capabilities (e.g. efficient training, fast generation, handling irregularly sampled data). This specialization is typically expressed via domain-specific preprocessing, inductive biases, and architectures. Sequence models provide a general framework for solving many of these problems with reduced specialization – e.g. Vision Transformers for image classification with less
2D information (Dosovitskiy et al., 2020). However, most models such as Transformers generally still require substantial specialization per task to achieve high performance.
Deep SSMs in particular have conceptual strengths that suggest they may be promising as a general sequence modeling solution. These strengths include a principled approach to handling LRDs, as well as the ability to move between continuous-time, convolutional, and recurrent model representations, each with distinct capabilities (Fig. 1). Our technical contributions enable SSMs to be applied successfully to a varied set of benchmarks with minimal modification:
• Large-scale generative modeling. On CIFAR-10 density estimation, S4 is competitive with the best autoregressive models (2.85 bits per dim). On WikiText-103 language modeling, S4 substantially closes the gap to Transformers (within 0.8 perplexity), setting SoTA for attention-free models.
• Fast autoregressive generation. Like RNNs, S4 can use its latent state to perform 60× faster pixel/token generation than standard autoregressive models on CIFAR-10 and WikiText-103.
• Irregularly sampled time-series. Like specialized CTMs, S4 can adapt to changes in time series sampling frequency without retraining, e.g. at 0.5× frequency on speech classification.
• Learning with weaker inductive biases. With no architectural changes, S4 surpasses Speech CNNs on speech classification, outperforms the specialized Informer model on 40/50 time-series forecasting problems, and matches a 2-D ResNet on sequential CIFAR with over 90% accuracy.

heading: 2 BACKGROUND: STATE SPACES
text: Sections 2.1 to 2.4 describe the four properties of SSMs in Fig. 1: the classic continuous-time representation, addressing LRDs with the HiPPO framework, the discrete-time recurrent representation, and the parallelizable convolution representation. In particular, Section 2.4 introduces the SSM convolution kernel K, which is the focus of our theoretical contributions in Section 3.

heading: 2.1 STATE SPACE MODELS: A CONTINUOUS-TIME LATENT STATE MODEL
text: The state space model is defined by the simple equation (1). It maps a 1-D input signal u(t) to an N -D latent state x(t) before projecting to a 1-D output signal y(t).
x′(t) = Ax(t) + Bu(t)
y(t) = Cx(t) + Du(t) (1)
SSMs are broadly used in many scientific disciplines and related to latent state models such as Hidden Markov Models (HMM). Our goal is to simply use the SSM as a black-box representation in a deep sequence model, where A,B,C,D are parameters learned by gradient descent. For the remainder of this paper, we will omit the parameter D for exposition (or equivalently, assume D = 0) because the term Du can be viewed as a skip connection and is easy to compute.

heading: 2.2 ADDRESSING LONG-RANGE DEPENDENCIES WITH HIPPO
text: Prior work found that the basic SSM (1) actually performs very poorly in practice. Intuitively, one explanation is that linear first-order ODEs solve to an exponential function, and thus may suffer from gradients scaling exponentially in the sequence length (i.e., the vanishing/exploding gradients problem (Pascanu et al., 2013)). To address this problem, the LSSL leveraged the HiPPO theory of continuous-time memorization (Gu et al., 2020a). HiPPO specifies a class of certain matrices A ∈ RN×N that when incorporated into (1), allows the state x(t) to memorize the history of the input u(t). The most important matrix in this class is defined by equation (2), which we will call the HiPPO matrix. For example, the LSSL found that simply modifying an SSM from a random matrix A to equation (2) improved its performance on the sequential MNIST benchmark from 50% to 98%.
(HiPPO Matrix) Ank =  (2n+ 1)1/2(2k + 1)1/2 if n > k n+ 1 if n = k 0 if n < k . (2)

heading: 2.3 DISCRETE-TIME SSM: THE RECURRENT REPRESENTATION
text: To be applied on a discrete input sequence (u0, u1, . . . ) instead of continuous function u(t), (1) must be discretized by a step size ∆ that represents the resolution of the input. Conceptually, the inputs uk can be viewed as sampling an implicit underlying continuous signal u(t), where uk = u(k∆).
To discretize the continuous-time SSM, we follow prior work in using the bilinear method (Tustin, 1947), which converts the state matrix A into an approximation A . The discrete SSM is
xk = Axk−1 + Buk A = (I −∆/2 ·A)−1(I + ∆/2 ·A) yk = Cxk B = (I −∆/2 ·A)−1∆B C = C.
(3)
Equation (3) is now a sequence-to-sequence map uk 7→ yk instead of function-to-function. Moreover the state equation is now a recurrence in xk, allowing the discrete SSM to be computed like an RNN. Concretely, xk ∈ RN can be viewed as a hidden state with transition matrix A. Notationally, throughout this paper we use A,B, . . . to denote discretized SSM matrices defined by (3). Note that these matrices are a function of both A as well as a step size ∆; we suppress this dependence for notational convenience when it is clear.

heading: 2.4 TRAINING SSMS: THE CONVOLUTIONAL REPRESENTATION
text: The recurrent SSM (3) is not practical for training on modern hardware due to its sequentiality. Instead, there is a well-known connection between linear time-invariant (LTI) SSMs such as (1) and continuous convolutions. Correspondingly, (3) can actually be written as a discrete convolution.
For simplicity let the initial state be x−1 = 0. Then unrolling (3) explicitly yields
x0 = Bu0 x1 = ABu0 + Bu1 x2 = A 2 Bu0 + ABu1 + Bu2 . . .
y0 = CBu0 y1 = CABu0 + CBu1 y2 = CA 2 Bu0 + CABu1 + CBu2 . . .
This can be vectorized into a convolution (4) with an explicit formula for the convolution kernel (5).
yk = CA k Bu0 + CA k−1 Bu1 + · · ·+ CABuk−1 + CBuk
y = K ∗ u. (4)
K ∈ RL := KL(A,B,C) := ( CA i B ) i∈[L] = (CB,CAB, . . . ,CA L−1 B). (5)
In other words, equation (4) is a single (non-circular) convolution and can be computed very efficiently with FFTs, provided that K is known. However, computing K in (5) is non-trivial and is the focus of our technical contributions in Section 3. We call K the SSM convolution kernel or filter.

heading: 3 METHOD: STRUCTURED STATE SPACES (S4)
text: Our technical results focus on developing the S4 parameterization and showing how to efficiently compute all views of the SSM (Section 2): the continuous representation (A,B,C) (1), the recurrent representation (A,B,C) (3), and the convolutional representation K (4).
Section 3.1 motivates our approach, which is based on the linear algebraic concepts of conjugation and diagonalization, and discusses why the naive application of this approach does not work. Section 3.2 gives an overview of the key technical components of our approach and formally defines the S4 parameterization. Section 3.3 sketches the main results, showing that S4 is asymptotically efficient (up to log factors) for sequence models. Proofs are in Appendices B and C.

heading: 3.1 MOTIVATION: DIAGONALIZATION
text: The fundamental bottleneck in computing the discrete-time SSM (3) is that it involves repeated matrix multiplication by A. For example, computing (5) naively as in the LSSL involves L successive multiplications by A, requiring O(N2L) operations and O(NL) space.
To overcome this bottleneck, we use a structural result that allows us to simplify SSMs.
Algorithm 1 S4 CONVOLUTION KERNEL (SKETCH)
Input: S4 parameters Λ,p, q,B,C ∈ CN and step size ∆ Output: SSM convolution kernel K = KL(A,B,C) for A = Λ− pq∗ (equation (5))
1: C ← C ( I −AL ) . Truncate SSM generating function (SSMGF) to length L
2: [ k00(ω) k01(ω) k10(ω) k11(ω) ] ← [ C q ]∗ ( 2 ∆ 1−ω 1+ω −Λ )−1 [B p] . Black-box Cauchy kernel
3: K̂(ω)← 21+ω [ k00(ω)− k01(ω)(1 + k11(ω))−1k10(ω) ] . Woodbury Identity 4: K̂ = {K̂(ω) : ω = exp(2πi kL )} . Evaluate SSMGF at all roots of unity ω ∈ ΩL 5: K ← iFFT(K̂) . Inverse Fourier Transform
Lemma 3.1. Conjugation is an equivalence relation on SSMs (A,B,C) ∼ (V −1AV ,V −1B,CV ).
Proof. Write out the two SSMs with state denoted by x and x̃ respectively:
x′ = Ax+ Bu x̃′ = V −1AV x̃+ V −1Bu
y = Cx y = CV x̃
After multiplying the right side SSM by V , the two SSMs become identical with x = V x̃. Therefore these compute the exact same operator u 7→ y, but with a change of basis by V in the state x.
Lemma 3.1 motivates putting A into a canonical form by conjugation1, which is ideally more structured and allows faster computation. For example, if A were diagonal, the resulting computations become much more tractable. In particular, the desired K (equation (4)) would be a Vandermonde product which theoretically only needs O((N +L) log2(N +L)) arithmetic operations (Pan, 2001). Unfortunately, the naive application of diagonalization does not work due to numerical issues. First, Vandermonde multiplication is itself a famously ill-conditioned problem (Pan, 2016). Furthermore, we derive the explicit diagonalization for the HiPPO matrix (2) and show it has entries exponentially large in the state sizeN , rendering the diagonalization numerically infeasible (e.g. CV in Lemma 3.1 would not be computable). We note that Gu et al. (2021) proposed a different (unimplemented) algorithm to compute K faster than the naive algorithm. In Appendix B, we prove that it is also numerically unstable for related reasons.
Lemma 3.2. The HiPPO matrix A in equation (2) is diagonalized by the matrix Vij = ( i+j i−j ) . In
particular, V3i,i = ( 4i 2i ) ≈ 24i. Therefore V has entries of magnitude up to 24N/3.

heading: 3.2 THE S4 PARAMETERIZATION: NORMAL PLUS LOW-RANK
text: The previous discussion implies that we should only conjugate by well-conditioned matrices V . The ideal scenario is when the matrix A is diagonalizable by a perfectly conditioned (i.e., unitary) matrix. By the Spectral Theorem of linear algebra, this is exactly the class of normal matrices. However, this class of matrices is restrictive; in particular, it does not contain the HiPPO matrix (2).
We make the observation that although the HiPPO matrix is not normal, it can be decomposed as the sum of a normal and low-rank matrix. However, this is still not useful by itself: unlike a diagonal matrix, powering up this sum (in (5)) is still slow and not easily optimized. We overcome this bottleneck by simultaneously applying three new techniques.
• Instead of computing K directly, we compute its spectrum by evaluating its truncated generating function ∑L−1 j=0 Kjζ j at the roots of unity ζ. K can then be found by applying an inverse FFT.
• This generating function is closely related to the matrix resolvent, and now involves a matrix inverse instead of power. The low-rank term can now be corrected by applying the Woodbury identity which reduces (A + pq∗)−1 in terms of A−1, truly reducing to the diagonal case.
• Finally, we show that the diagonal matrix case is equivalent to the computation of a Cauchy kernel 1
ωj−ζk , a well-studied problem with stable near-linear algorithms (Pan, 2015; 2017).
1Note that although we ultimately care about A, conjugation commutes with discretization.
Our techniques apply to any matrix that can be decomposed as Normal Plus Low-Rank (NPLR). Theorem 1. All HiPPO matrices from (Gu et al., 2020a) have a NPLR representation
A = V ΛV ∗ − pq> = V (Λ− V ∗p(V ∗q)∗)V ∗ (6) for unitary V ∈ CN×N , diagonal Λ, and low-rank factorization p, q ∈ RN×r. These matrices HiPPO- LegS, LegT, LagT all satisfy r = 1 or r = 2. In particular, equation (2) is NPLR with r = 1.

heading: 3.3 S4 ALGORITHMS AND COMPUTATIONAL COMPLEXITY
text: By equation (6), note that NPLR matrices can be conjugated into diagonal plus low-rank (DPLR) form (now overC instead ofR). Theorems 2 and 3 describe the complexities of SSMs where A is in DPLR form. S4 is optimal or near-optimal for both recurrent and convolutional representations. Theorem 2 (S4 Recurrence). Given any step size ∆, computing one step of the recurrence (3) can be done in O(N) operations where N is the state size.
Theorem 2 follows from the fact that the inverse of a DPLR matrix is also DPLR (e.g. also by the Woodbury identity). This implies that the discretized matrix A is the product of two DPLR matrices and thus has O(N) matrix-vector multiplication. Appendix C.2 computes A in closed DPLR form. Theorem 3 (S4 Convolution). Given any step size ∆, computing the SSM convolution filter K can be reduced to 4 Cauchy multiplies, requiring only Õ(N + L) operations and O(N + L) space.
Appendix C, Definition 3 formally defines Cauchy matrices, which are related to rational interpolation problems. Computing with Cauchy matrices is an extremely well-studied problem in numerical analysis, with both fast arithmetic and numerical algorithms based on the famous Fast Multipole Method (FMM) (Pan, 2001; 2015; 2017). The computational complexities of these algorithms under various settings are described in Appendix C, Proposition 5.
We reiterate that Theorem 3 is our core technical contribution, and its algorithm is the very motivation of the S4 parameterization. This algorithm is formally sketched in Algorithm 1.

heading: 3.4 ARCHITECTURE DETAILS OF THE DEEP S4 LAYER
text: Concretely, an S4 layer is parameterized as follows. First initialize a SSM with A set to the HiPPO matrix (2). By Lemma 3.1 and Theorem 1, this SSM is unitarily equivalent to some (Λ−pq∗,B,C) for some diagonal Λ and vectors p, q,B,C ∈ CN×1. These comprise S4’s 5N trainable parameters. The overall deep neural network (DNN) architecture of S4 is similar to prior work. As defined above, S4 defines a map from RL → RL, i.e. a 1-D sequence map. Typically, DNNs operate on feature maps of sizeH instead of 1. S4 handles multiple features by simply definingH independent copies of itself, and then mixing theH features with a position-wise linear layer for a total ofO(H2)+O(HN) parameters per layer. Nonlinear activation functions are also inserted between these layers. Overall, S4 defines a sequence-to-sequence map of shape (batch size, sequence length, hidden dimension), exactly the same as related sequence models such as Transformers, RNNs, and CNNs.

heading: 4 EXPERIMENTS
text: Section 4.1 benchmarks S4 against the LSSL and efficient Transformer models. Section 4.2 validates S4 on LRDs: the LRA benchmark and raw speech classification. Section 4.3 investigates whether S4 can be used as a general sequence model to perform effectively and efficiently in a wide variety of settings including image classification, image and text generation, and time series forecasting.

heading: 4.1 S4 EFFICIENCY BENCHMARKS
text: We benchmark that S4 can be trained quickly and efficiently, both compared to the LSSL, as well as efficient Transformer variants designed for long-range sequence modeling. As outlined in Section 3, S4 is theoretically much more efficient than the LSSL, and Table 1 confirms that the S4 is orders of magnitude more speed- and memory-efficient for practical layer sizes. In fact, S4’s speed and memory use is competitive with the most efficient Transformer variants benchmarked by Tay et al. (2021)— Linear Transformer (Katharopoulos et al., 2020) and Performer (Choromanski et al., 2020)—in a parameter-matched setting (Table 2, following the protocol of Tay et al. (2021)).
Table 1: Deep SSMs: The S4 parameterization with Algorithm 1 is asymptotically more efficient than the LSSL.
Table 2: Benchmarks vs. efficient Transformers

heading: 4.2 LEARNING LONG RANGE DEPENDENCIES
text: As described in Sections 2.2 and 3.1, S4 uses a principled approach to address LRDs based on the HiPPO theory of continuous-time memorization. Our goal in this section is to validate that S4 achieves high performance on difficult tasks that require long-range reasoning. We focus here on two problems: (i) the Long-Range Arena, a well-known benchmark designed to test efficient sequence models on LRDs, and (ii) a speech classification problem as a real-world test of LRDs.
Long Range Arena (LRA). LRA (Tay et al., 2021) contains 6 tasks with lengths 1K-16K steps, encompassing modalities and objectives that require similarity, structural, and visuospatial reasoning. Table 3 compares S4 against the 11 Transformer variants from Tay et al. (2021) as well as follow-up work. S4 substantially advances the SoTA, outperforming all baselines on all tasks and averaging 80.88% compared to less than 60% for every baseline. Notably, S4 solves the Path-X task, an extremely challenging task that involves reasoning about LRDs with sparse signal over sequences of length 128× 128 = 16384. All previous models have failed (i.e. random guessing) due to memory or computation bottlenecks, or simply being unable to learn such long dependencies.
We analyze S4’s performance on Path-X by visualizing its learned representations, in particular 1-D convolution kernels K which are the focus of our technical results in Section 3. Fig. 2 shows that S4
learns a variety of filters that display spatially consistent structure and demonstrate awareness of the 2-D nature of the data. In particular, the lower layers learn simple kernels that extract features from just a few rows of local context while ignoring the rest of the image. On the other hand, higher layers aggregate information globally across full columns of the image at varying spatial frequencies. Filters in these higher layers span the entire context (16384 pixels), confirming S4’s ability to learn LRDs.
Raw Speech Classification. Speech is a typical real-world time series domain, involving signals sampled from an underlying physical process at high frequency. We perform speech classification using the Speech Commands dataset (Warden, 2018). While most sequence models for speech rely on extensive preprocessing (e.g. to MFCC features), we classify raw speech (length-16000) following Romero et al. (2021). S4 achieves 98.3% accuracy, higher than all baselines that use the 100× shorter MFCC features, and validates that a powerful LRD model is able to extract more information from the raw data and outperform hand-crafted pre-processing. Additionally, we include a baseline CNN specifically designed for raw speech, the discriminator from the WaveGAN model (Donahue et al., 2019), which performs worse than S4 while having 90× more parameters and incorporating many more architectural heuristics (Appendix D.2).

heading: 4.3 S4 AS A GENERAL SEQUENCE MODEL
text: A key goal of sequence modeling research is to develop a single model that can be applied in many domains (e.g. images, audio, text, time-series) with a broad range of capabilities (e.g. efficient training, fast generation, handling irregularly sampled data). As a fundamental scientific model, SSMs are a promising candidate that come with a range of capabilities, and S4’s strong results on LRD benchmarks spanning images, text, and speech are evidence of S4’s potential as a general sequence model. In this section, we focus on understanding this question in more depth by highlighting key strengths of S4 in settings that usually require specialized models. The tasks we focus on (generative modeling, image classification, time-series forecasting) are considered as LRD tasks in the literature, and serve as additional validation that S4 handles LRDs efficiently.
Large-scale generative modeling. We investigate two well-studied image and text benchmarks to validate the scalability, flexibility, and efficiency of S4. These tasks require much larger models than our previous tasks – up to 250M parameters.
First, CIFAR density estimation is a popular benchmark for autoregressive models, where images are flattened into a sequence of 3072 RGB subpixels that are predicted one by one. Table 6 shows that with no 2D inductive bias, S4 is competitive with the best models designed for this task.
Second, WikiText-103 is an established benchmark for language modeling, an important task for large-scale sequence models where tokens are predicted sequentially based on past context. Although RNNs were the model of choice for many years, Transformers are now the dominant model in such applications that contain data that is inherently discrete. We show that alternative models
to Transformers can still be competitive in these settings. By simply taking a strong Transformer baseline (Baevski & Auli, 2018) and replacing the self-attention layers, S4 substantially closes the gap to Transformers (within 0.8 ppl), setting SoTA for attention-free models by over 2 ppl.
Fast autoregressive inference. A prominent limitation of autoregressive models is inference speed (e.g. generation), since they require a pass over the full context for every new sample. Several methods have been specifically crafted to overcome this limitation such as the Linear Transformer, a hybrid Transformer/RNN that switches to a stateful, recurrent view at inference time for speed.
As a stateful model, SSMs automatically have this ability (Fig. 1). By switching to its recurrent representation (Section 2.3), S4 requires constant memory and computation per time step – in contrast to standard autoregressive models which scale in the context length. On both CIFAR-10 and WikiText103, we report the throughput of various models at generation time, with S4 around 60× faster than a vanilla Transformer on both tasks (details in Appendix D.3.3).
Sampling resolution change. As a continuous-time model, S4 automatically adapts to data sampled at different rates, a challenging setting for time series with a dedicated line of work (Rubanova et al., 2019; De Brouwer et al., 2019; Romero et al., 2021). Without re-training, S4 achieves 96.3% accuracy at 0.5× the frequency on Speech Commands (Table 4), simply by changing its internal step size ∆ (Section 2.3).
Learning with weaker inductive bias. Beyond our results on speech (Section 4.2), we further validate that S4 can be applied with minimal modifications on two domains that typically require specialized domain-specific preprocessing and architectures. First, we compare S4 to the Informer (Zhou et al., 2021), a new Transformer architecture that uses a complex encoder-decoder designed for time-series forecasting problems. A simple application of S4 that treats forecasting as a generic sequence-to-sequence transformation (Appendix D.3.5) outperforms the Informer and other baselines on 40/50 settings across 5 forecasting tasks. Notably, S4 is better on the longest setting in each task, e.g. reducing MSE by 32% when forecasting 30 days of weather data (Appendix D.3.5).
Finally, we evaluate S4 on pixel-level sequential image classification tasks (Table 5), popular benchmarks which were originally LRD tests for RNNs (Arjovsky et al., 2016). Beyond LRDs, these benchmarks point to a recent effort of the ML community to solve vision problems with reduced domain knowledge, in the spirit of models such as Vision Transformers (Dosovitskiy et al., 2020) and MLP-Mixer (Tolstikhin et al., 2021). Sequential CIFAR is a particularly challenging dataset where outside of SSMs, all sequence models have a gap of over 25% to a simple 2-D CNN. By contrast, S4 is competitive with a larger ResNet18 (7.9M vs. 11.0M parameters), both with (93.16% vs. 95.62%) or without (91.12% vs. 89.46%) data augmentation. Moreover, it is much more robust to other architectural choices (e.g. 90.46% vs. 79.52% when swapping BatchNorm for LayerNorm).

heading: 5 CONCLUSION
text: We introduce S4, a sequence model that uses a new parameterization for the state space model’s continuous-time, recurrent, and convolutional views to efficiently model LRDs in a principled manner. Results across established benchmarks evaluating a diverse range of data modalities and model capabilities suggest that S4 has the potential to be an effective general sequence modeling solution.

heading: A DISCUSSION
text: Related Work. Our work is most closely related to a line of work originally motivated by a particular biologically-inspired SSM, which led to mathematical models for addressing LRDs. Voelker (2019); Voelker et al. (2019) derived a non-trainable SSM motivated from approximating a neuromorphic spiking model, and Chilkuri & Eliasmith (2021) showed that it could be sped up at train time with a convolutional view. Gu et al. (2020a) extended this special case to a general continuous-time function approximation framework with several more special cases of A matrices designed for long-range dependencies. However, instead of using a true SSM, all of these works fixed a choice of A and built RNNs around it. Most recently, Gu et al. (2021) used the full (1) explicitly as a deep SSM model, exploring new conceptual views of SSMs, as well as allowing A to be trained. As mentioned in Section 1, their method used a naive instantiation of SSMs that suffered from an additional factor of N in memory and N2 in computation.
Beyond this work, our technical contributions (Section 3) on the S4 parameterization and algorithms are applicable to a broader family of SSMs including these investigated in prior works, and our techniques for working with these models may be of independent interest.
Implementation. The computational core of S4’s training algorithm is the Cauchy kernel discussed in Sections 3.2 and 3.3 and Appendix C.3. As described in Appendix C.3 Proposition 5, there are many algorithms for it with differing computational complexities and sophistication. Our current implementation of S4 actually uses the naive O(NL) algorithm which is easily parallelized on GPUs and has more easily accessible libraries allowing it to be implemented; we leverage the pykeops library for memory-efficient kernel operations. However, this library is a much more general library
that may not be optimized for the Cauchy kernels used here, and we believe that a dedicated CUDA implementation can be more efficient. Additionally, as discussed in this work, there are asymptotically faster and numerically stable algorithms for the Cauchy kernel (Proposition 5). However, these algorithms are currently not implemented for GPUs due to a lack of previous applications that require them. We believe that more efficient implementations of these self-contained computational kernels are possible, and that S4 (and SSMs at large) may have significant room for further improvements in efficiency.
Limitations and Future Directions. In this work, we show that S4 can address a wide variety of data effectively. However, it may not necessarily be the most suitable model for all types of data. For example, Table 7 still found a gap compared to Transformers for language modeling. An interesting future direction is exploring combinations of S4 with other sequence models to complement their strengths. We are excited about other directions, including continuing to explore the benefits of S4 on audio data (e.g. pre-training or generation settings), and generalizing HiPPO and S4 to higher-dimensional data for image and video applications.

heading: B NUMERICAL INSTABILITY OF LSSL
text: This section proves the claims made in Section 3.1 about prior work. We first derive the explicit diagonalization of the HiPPO matrix, confirming its instability because of exponentially large entries. We then discuss the proposed theoretically fast algorithm from (Gu et al., 2021) (Theorem 2) and show that it also involves exponentially large terms and thus cannot be implemented.

heading: B.1 HIPPO DIAGONALIZATION
text: Proof of Lemma 3.2. The HiPPO matrix (2) is equal, up to sign and conjugation by a diagonal matrix, to
A =  1 −1 2 1 −3 3 −1 3 −5 4 1 −3 5 −7 5 −1 3 −5 7 −9 6 1 −3 5 −7 9 −11 7 −1 3 −5 7 −9 11 −13 8
... . . .  Ank =  (−1)n−k(2k + 1) n > k k + 1 n = k
0 n < k
.
Our goal is to show that this A is diagonalized by the matrix
V =
( i+ j
i− j ) ij =  1 1 1 1 3 1 1 6 5 1 1 10 15 7 1 1 15 35 28 9 1 ... . . .  ,
or in other words that columns of this matrix are eigenvectors of A.
Concretely, we will show that the j-th column of this matrix v(j) with elements
v (j) i = { 0 i < j( i+j i−j ) = ( i+j 2j ) i ≥ j
is an eigenvector with eigenvalue j + 1. In other words we must show that for all indices k ∈ [N ], (Av(j))k = ∑ i Akiv (j) i = (j + 1)v (j) k . (7)
If k < j, then for all i inside the sum, either k < i or i < j. In the first case Aki = 0 and in the second case v(j)i = 0, so both sides of equation (7) are equal to 0.
It remains to show the case k ≥ j, which proceeds by induction on k. Expanding equation (7) using the formula for A yields
(Av) (j) k = ∑ i Akiv (j) i = k−1∑ i=j (−1)k−i(2i+ 1) ( i+ j 2j ) + (k + 1) ( k + j 2j ) .
In the base case k = j, the sum disappears and we are left with (Av(j))j = (j+1) ( 2j 2j ) = (j+1)v (j) j , as desired.
Otherwise, the sum for (Av)(j)k is the same as the sum for (Av) (j) k−1 but with sign reversed and a few edge terms. The result follows from applying the inductive hypothesis and algebraic simplification:
(Av) (j) k = −(Av) (j) k−1 − (2k − 1)
( k − 1 + j
2j
) + k ( k − 1 + j
2j
) + (k + 1) ( k + j
2j ) = −(j + 1) ( k − 1 + j
2j
) − (k − 1) ( k − 1 + j
2j
) + (k + 1) ( k + j
2j ) = −(j + k) ( k − 1 + j
2j
) + (k + 1) ( k + j
2j ) = −(j + k) (k − 1 + j)!
(k − 1− j)!(2j)! + (k + 1)
( k + j
2j ) = − (k + j)!
(k − 1− j)!(2j)! + (k + 1)
( k + j
2j ) = −(k − j) (k + j)!
(k − j)!(2j)! + (k + 1)
( k + j
2j ) = (j − k)(k + 1) ( k + j
2j
) + (k + 1) ( k + j
2j ) = (j + 1)v
(j) k .

heading: B.2 FAST BUT UNSTABLE LSSL ALGORITHM
text: Instead of diagonalization, Gu et al. (2021, Theorem 2) proposed a sophisticated fast algorithm to compute
KL(A,B,C) = (CB,CAB, . . . ,CA L−1 B).
This algorithm runs in O(N log2N + L logL) operations and O(N + L) space. However, we now show that this algorithm is also numerically unstable.
There are several reasons for the instability of this algorithm, but most directly we can pinpoint a particular intermediate quantity that they use. Definition 1. The fast LSSL algorithm computes coefficients of p(x), the characteristic polynomial of A, as an intermediate computation. Additionally, it computes the coefficients of its inverse, p(x)−1 (mod xL).
We now claim that this quantity is numerically unfeasible. We narrow down to the case when A = I is the identity matrix. Note that this case is actually in some sense the most typical case: when
discretizing the continuous-time SSM to discrete-time by a step-size ∆, the discretized transition matrix A is brought closer to the identity. For example, with the Euler discretization A = I + ∆A, we have A→ I as the step size ∆→ 0. Lemma B.1. When A = I , the fast LSSL algorithm requires computing terms exponentially large in N .
Proof. The characteristic polynomial of I is
p(x) = det |I − xI| = (1− x)N .
These coefficients have size up to ( N N 2 ) ≈ 2 N√ πN/2 .
The inverse of p(x) has even larger coefficients. It can be calculated in closed form by the generalized binomial formula:
(1− x)−N = ∞∑ k=0 ( N + k − 1 k ) xk.
Taking this (mod xL), the largest coefficient is( N + L− 2 L− 1 ) = ( N + L− 2 N − 1 ) = (L− 1)(L− 2) . . . (L−N + 1) (N − 1)! .
When L = N − 1 this is ( 2(N − 1) N − 1 ) ≈ 2 2N √ πN
already larger than the coefficients of (1− x)N , and only increases as L grows.

heading: C S4 ALGORITHM DETAILS
text: This section proves the results of Section 3.3, providing complete details of our efficient algorithms for S4.
Appendices C.1 to C.3 prove Theorems 1 to 3 respectively.

heading: C.1 NPLR REPRESENTATIONS OF HIPPO MATRICES
text: We first prove Theorem 1, showing that all HiPPO matrices for continuous-time memory fall under the S4 normal plus low-rank (NPLR) representation.
Proof of Theorem 1. We consider each of the three cases HiPPO-LagT, HiPPO-LegT, and HiPPOLegS separately. Note that the primary HiPPO matrix defined in this work (equation (2)) is the HiPPO-LegT matrix.
HiPPO-LagT. The HiPPO-LagT matrix is simply
Ank =  0 n < k
− 12 n = k −1 n > k
A = −  1 2 . . . 1 12 1 1 12 1 1 1 12 ... . . .  .
Adding the matrix of all 12 , which is rank 1, yields
−  − 12 − 1 2 − 1 2 1 2 − 1 2 − 1 2
1 2
1 2 − 1 2
1 2
1 2
1 2  . This matrix is now skew-symmetric. Skew-symmetric matrices are a particular case of normal matrices with pure-imaginary eigenvalues.
Gu et al. (2020a) also consider a case of HiPPO corresponding to the generalized Laguerre polynomials that generalizes the above HiPPO-LagT case. In this case, the matrix A (up to conjugation by a diagonal matrix) ends up being close to the above matrix, but with a different element on the diagonal. After adding the rank-1 correction, it becomes the above skew-symmetric matrix plus a multiple of the identity. Thus after diagonalization by the same matrix as in the LagT case, it is still reduced to diagonal plus low-rank (DPLR) form, where the diagonal is now pure imaginary plus a real constant.
HiPPO-LegS. We restate the formula from equation (2) for convenience.
Ank = −  (2n+ 1)1/2(2k + 1)1/2 if n > k n+ 1 if n = k 0 if n < k .
Adding 12 (2n+ 1) 1/2(2k + 1)1/2 to the whole matrix gives
−  1 2 (2n+ 1)
1/2(2k + 1)1/2 if n > k 1 2 if n = k − 12 (2n+ 1) 1/2(2k + 1)1/2 if n < k
Note that this matrix is not skew-symmetric, but is 12I + S where S is a skew-symmetric matrix. This is diagonalizable by the same unitary matrix that diagonalizes S.

heading: HiPPO-LegT.
text: Up to the diagonal scaling, the LegT matrix is
A = −  1 −1 1 −1 . . . 1 1 −1 1 1 1 1 −1 1 1 1 1 ... . . .  . By adding −1 to this matrix and then the matrix2 2
2 2  the matrix becomes  −2 −22 −2
2 2  which is skew-symmetric. In fact, this matrix is the inverse of the Chebyshev Jacobi.
An alternative way to see this is as follows. The LegT matrix is the inverse of the matrix−1 1 0−1 1−1 1 −1 −1 
This can obviously be converted to a skew-symmetric matrix by adding a rank 2 term. The inverses of these matrices are also rank-2 differences from each other by the Woodbury identity.

heading: A final form is −1 1 −1 1−1 −1 1 −1−1 −1 −1 1
text: −1 −1 −1 −1 + 1 0 1 00 1 0 11 0 1 0 0 1 0 1  =  0 1 0 1−1 0 1 00 −1 0 1 −1 0 −1 0  This has the advantage that the rank-2 correction is symmetric (like the others), but the normal skew-symmetric matrix is now 2-quasiseparable instead of 1-quasiseparable.

heading: C.2 COMPUTING THE S4 RECURRENT VIEW
text: We prove Theorem 2 showing the efficiency of the S4 parameterization for computing one step of the recurrent representation (Section 2.3).
Recall that without loss of generality, we can assume that the state matrix A = Λ− pq∗ is diagonal plus low-rank (DPLR), potentially over C. Our goal in this section is to explicitly write out a closed form for the discretized matrix A.
Recall from equation (3) that
A = (I −∆/2 ·A)−1(I + ∆/2 ·A) B = (I −∆/2 ·A)−1∆B.
We first simplify both terms in the definition of A independently.
Forward discretization. The first term is essentially the Euler discretization motivated in Section 2.3.
I + ∆
2 A = I +
∆ 2 (Λ− pq∗)
= ∆
2
[ 2
∆ I + (Λ− pq∗) ] = ∆
2 A0
where A0 is defined as the term in the final brackets.
Backward discretization. The second term is known as the Backward Euler’s method. Although this inverse term is normally difficult to deal with, in the DPLR case we can simplify it using Woodbury’s Identity (Proposition 4).(
I − ∆ 2 A
)−1 = ( I − ∆
2 (Λ− pq∗) )−1 = 2
∆
[ 2
∆ −Λ + pq∗ ]−1 = 2
∆
[ D −Dp (1 + q∗Dp)−1 q∗D ] = 2
∆ A1
where D = (
2 ∆ −Λ )−1 and A1 is defined as the term in the final brackets. Note that (1 + q∗Dp)
is actually a scalar, and A1 is actually itself a DPLR matrix.
S4 Recurrence. Finally, the full bilinear discretization can be rewritten in terms of these matrices as
A = A1A0
B = 2
∆ A1∆B = 2A1B.
The discrete-time SSM (3) becomes
xk = Axk−1 + Buk = A1A0xk−1 + 2A1Buk
yk = Cxk.
Note that A0,A1 are accessed only through matrix-vector multiplications. Since they are both DPLR, they have O(N) matrix-vector multiplication, showing Theorem 2.

heading: C.3 COMPUTING THE CONVOLUTIONAL VIEW
text: The most involved part of using SSMs efficiently is computing K. This algorithm was sketched in Section 3.2 and is the main motivation for the S4 parameterization. In this section, we define the necessary intermediate quantities and prove the main technical result.
The algorithm for Theorem 3 falls in roughly three stages. Assuming A has been conjugated into diagonal plus low-rank form, we successively simplify the problem of computing K by applying the techniques outlined in Section 3.2.
Reduction 0: Diagonalization By Lemma 3.1, we can switch the representation by conjugating with any unitary matrix. For the remainder of this section, we can assume that A is (complex) diagonal plus low-rank (DPLR).
Note that unlike diagonal matrices, a DPLR matrix does not lend itself to efficient computation of K. The reason is that K computes terms CA i B which involve powers of the matrix A. These are
trivially computable when A is diagonal, but is no longer possible for even simple modifications to diagonal matrices such as DPLR.
Reduction 1: SSM Generating Function To address the problem of computing powers of A, we introduce another technique. Instead of computing the SSM convolution filter K directly, we introduce a generating function on its coefficients and compute evaluations of it.
Definition 2 (SSM Generating Function). We define the following quantities:
• The SSM convolution function is K(A,B,C) = (CB,CAB, . . . ) and the (truncated) SSM filter of length L
KL(A,B,C) = (CB,CAB, . . . ,CA L−1 B) ∈ RL (8)
• The SSM generating function at node z is
K̂(z;A,B,C) ∈ C := ∞∑ i=0 CA i Bzi = C(I −Az)−1B (9)
and the truncated SSM generating function at node z is
K̂L(z;A,B,C) ∈ C := L−1∑ i=0 CA i Bzi = C(I −ALzL)(I −Az)−1B (10)
• The truncated SSM generating function at nodes Ω ∈ CM is K̂L(Ω;A,B,C) ∈ CM := ( K̂L(ωk;A,B,C) ) k∈[M ]
(11)
Intuitively, the generating function essentially converts the SSM convolution filter from the time domain to frequency domain. Importantly, it preserves the same information, and the desired SSM convolution filter can be recovered from evaluations of its generating function.
Lemma C.1. The SSM function KL(A,B,C) can be computed from the SSM generating function KL(Ω;A,B,C at the roots of unity Ω = {exp(2π kL : k ∈ [L]} stably in O(L logL) operations.
Proof. For convenience define
K = KL(A,B,C) K̂ = K̂L(Ω;A,B,C).
Note that
K̂j = L−1∑ k=0 Kk exp(2πi jk L ).
Note that this is exactly the same as the Discrete Fourier Transform (DFT):
K̂ = FLK.
Therefore K can be recovered from K̂ with a single inverse DFT, which requires O(L logL) operations with the Fast Fourier Transform (FFT) algorithm.
Reduction 2: Woodbury Correction The primary motivation of Definition 2 is that it turns powers of A into a single inverse of A (equation (9)). While DPLR cannot be powered efficiently due to the low-rank term, they can be inverted efficiently by the well-known Woodbury identity. Proposition 4 (Binomial Inverse Theorem or Woodbury matrix identity Woodbury (1950); Golub & Van Loan (2013)). Over a commutative ring R, let A ∈ RN×N and U ,V ∈ RN×p. Suppose A and A + UV ∗ are invertible. Then Ip + V ∗A−1U ∈ Rp×p is invertible and
(A + UV ∗)−1 = A−1 −A−1U(Ip + V ∗A−1U)−1V ∗A−1
With this identity, we can convert the SSM generating function on a DPLR matrix A into one on just its diagonal component. Lemma C.2. Let A = Λ−pq∗ be a diagonal plus low-rank representation. Then for any evaluation node z ∈ C, the truncated generating function satisfies
KL(z) = 2
1 + z
[ C̃∗R(z)B − C̃∗R(z)p (1 + q∗R(z)p)−1 q∗R(z)B ] C̃ = C(I −AL)
R(z; Λ) =
( 2
∆ 1− z 1 + z
−Λ )−1 .
Proof. Directly expanding Definition 2 yields
KL(z;A,B,C) = C ∗ B + C ∗ ABz + · · ·+ C∗AL−1BzL−1 = C ∗ ( I −AL ) ( I −Az )−1 B
= C̃∗ ( I −Az )−1 B
where C̃∗ = C∗ ( I −AL ) .
We can now explicitly expand the discretized SSM matrices A and B back in terms of the original SSM parameters A and B. Lemma C.3 provides an explicit formula, which allows further simplifying
C̃∗ ( I −Az )−1 B = 2 1 + z C̃∗ ( 2 ∆ 1− z 1 + z −A )−1 B
= 2 1 + z C̃∗ ( 2 ∆ 1− z 1 + z − Ã + pq∗ )−1 B
= 2
1 + z
[ C̃∗R(z)B − C̃∗R(z)p (1 + q∗R(z)p)−1 q∗R(z)B ] .
The last line applies the Woodbury Identity (Proposition 4) where R(z) = (
2 ∆ 1−z 1+z − Ã
)−1 .
The previous proof used the following self-contained result to back out the original SSM matrices from the discretization.
Lemma C.3. Let A,B be the SSM matrices A,B discretized by the bilinear discretization with step size ∆. Then
C∗ ( I −A )−1 B = 2∆ 1 + z C∗ [ 2 1− z 1 + z −∆A ]−1 B
Proof. Recall that the bilinear discretization that we use (equation (3)) is
A = ( I − ∆
2 A
)−1( I + ∆
2 A ) B = ( I − ∆
2 A
)−1 ∆B
The result is proved algebraic manipulations.
C∗ ( I −Az )−1 B = C∗ [( I − ∆
2 A
)−1( I − ∆
2 A
) − ( I − ∆
2 A
)−1( I + ∆
2 A
) z ]−1 B
= C∗ [(
I − ∆ 2 A
) − ( I + ∆
2 A
) z ]−1( I − ∆
2 A
) B∆
= C∗ [ I − z − ∆
2 A(1 + z)
]−1 ∆B
= ∆
1− z C∗
[ I − ∆A
2 1−z1+z
]−1 B
= 2∆ 1 + z C∗ [ 2 1− z 1 + z −∆A ]−1 B
Note that in the S4 parameterization, instead of constantly computing C̃ = C ( I −AL ) , we can
simply reparameterize our parameters to learn C̃ directly instead of C, saving a minor computation cost and simplifying the algorithm.
Reduction 3: Core Cauchy Kernel We have reduced the original problem of computing K to the problem of computing the SSM generating function K̂(Ω;A,B,C) in the case that A is a diagonal matrix. We show that this is exactly the same as a Cauchy kernel, which is a well-studied problem with fast and stable numerical algorithms.
Definition 3. A Cauchy matrix or kernel on nodes Ω = (ωi) ∈ CM and Λ = (λj) ∈ CN is
M ∈ CM×N = M(Ω,Λ) = (Mij)i∈[M ],j∈[N ] Mij = 1
ωi − λj .
The computation time of a Cauchy matrix-vector product of size M ×N is denoted by C(M,N).
Computing with Cauchy matrices is an extremely well-studied problem in numerical analysis, with both fast arithmetic algorithms and fast numerical algorithms based on the famous Fast Multiple Method (FMM) (Pan, 2001; 2015; 2017).
Proposition 5 (Cauchy). A Cauchy kernel requires O(M +N) space, and operation count
C(M,N) =  O (MN) naively O ( (M +N) log2(M +N) ) in exact arithmetic
O ( (M +N) log(M +N) log 1ε ) numerically to precision ε.
Corollary C.4. Evaluating q∗R(Ω; Λ)p (defined in Lemma C.2) for any set of nodes Ω ∈ CL, diagonal matrix Λ, and vectors p, q can be computed in C(L,N) operations and O(L+N) space, where C(L,N) = Õ(L+N) is the cost of a Cauchy matrix-vector multiplication. Proof. For any fixed ω ∈ Ω, we want to compute ∑ j q∗j pj ω−λj . Computing this over all ωi is therefore exactly a Cauchy matrix-vector multiplication.
This completes the proof of Theorem 3. In Algorithm 1, note that the work is dominated by Step 2, which has a constant number of calls to a black-box Cauchy kernel, with complexity given by Proposition 5.

heading: D EXPERIMENT DETAILS AND FULL RESULTS
text: This section contains full experimental procedures and extended results and citations for our experimental evaluation in Section 4. Appendix D.1 corresponds to benchmarking results in Section 4.1, Appendix D.2 corresponds to LRD experiments (LRA and Speech Commands) in Section 4.2, and Appendix D.3 corresponds to the general sequence modeling experiments (generation, image classification, forecasting) in Section 4.3.

heading: D.1 BENCHMARKING
text: Benchmarking results from Table 1 and Table 2 were tested on a single A100 GPU.
Benchmarks against LSSL For a given dimension H , a single LSSL or S4 layer was constructed with H hidden features. For LSSL, the state size N was set to H as done in (Gu et al., 2021). For S4, the state size N was set to parameter-match the LSSL, which was a state size of N4 due to differences in the parameterization. Table 1 benchmarks a single forward+backward pass of a single layer.
Benchmarks against Efficient Transformers Following (Tay et al., 2021), the Transformer models had 4 layers, hidden dimension 256 with 4 heads, query/key/value projection dimension 128, and batch size 32, for a total of roughly 600k parameters. The S4 model was parameter tied while keeping the depth and hidden dimension constant (leading to a state size of N = 256).
We note that the relative orderings of these methods can vary depending on the exact hyperparameter settings. Other settings found S4 faster and more memory efficient than all Transformer variants, for example when the query/key/value dimension was set to 256 instead of 128, as is more standard, and when the batch size is increased. However, we chose to report a setting closer to the ones in (Tay et al., 2021).

heading: D.2 LONG-RANGE DEPENDENCIES
text: This section includes information for reproducing our experiments on the Long-Range Arena and Speech Commands long-range dependency tasks.
Long Range Arena Table 8 contains extended results table with all 11 methods considered in (Tay et al., 2021).
For the S4 model, hyperparameters for all datasets are reported in Table 9. For all datasets, we used the AdamW optimizer with a constant learning rate schedule. The S4 state size was always fixed to N = 64. As S4 is a sequence-to-sequence model with output shape (batch, length, dimension) and LRA tasks are classification, mean pooling along the length dimension was applied after the last layer.
We note that most of these results were trained for far longer than what was necessary to achieve SotA results (e.g., the Image task reaches SotA in 1 epoch). Results often keep improving with longer training times.
Hardware. All models were run on single GPU. Some tasks used an A100 GPU (notably, the Path-X experiments), which has a larger max memory of 40Gb. To reproduce these on smaller GPUs, the batch size can be reduced or gradients can be accumulated for two batches.
Path-X. We remark that an earlier version of this paper reported a higher score for Path-X. This earlier version used a different variant of the dataset, due to a misunderstanding of the properties of the dataset. More specifically, we found that S4 scored 93.68% accuracy on a version that involved taking the 256× 256 resolution version of the dataset and averaging every 2× 2 square; we erroneously thought that this version of the dataset was equivalent to the original Path-X.
After discussions with the LRA authors, we discovered that this is not equivalent to the 128× 128 resolution dataset (the correct Path-X), which is in fact much harder. In fact, Path-X is so difficult that a 2D ResNet without global receptive field (e.g. ResNet-18 or ResNet-34) cannot achieve above chance. This fact led to the original misunderstanding, as we could not solve this image classification task even with a ResNet and thought the data might have errors.
Speech Commands We provide details of sweeps run for baseline methods run by us—numbers for all others method are taken from Gu et al. (2021). The best hyperparameters used for S4 are included in Table 9.
Transformer (Vaswani et al., 2017) For MFCC, we swept the number of model layers {2, 4}, dropout {0, 0.1} and learning rates {0.001, 0.0005}. We used 8 attention heads, model dimension 128, prenorm, positional encodings, and trained for 150 epochs with a batch size of 100. For Raw, the Transformer model’s memory usage made training impossible.
Performer (Choromanski et al., 2020) For MFCC, we swept the number of model layers {2, 4}, dropout {0, 0.1} and learning rates {0.001, 0.0005}. We used 8 attention heads, model dimension 128, prenorm, positional encodings, and trained for 150 epochs with a batch size of 100. For Raw, we used a model dimension of 128, 4 attention heads, prenorm, and a batch size of 16. We reduced the
number of model layers to 4, so the model would fit on the single GPU. We trained for 100 epochs with a learning rate of 0.001 and no dropout.
ExpRNN (Lezcano-Casado & Martı́nez-Rubio, 2019) For MFCC, we swept hidden sizes {256, 512} and learning rates {0.001, 0.002, 0.0005}. Training was run for 200 epochs, with a single layer model using a batch size of 100. For Raw, we swept hidden sizes {32, 64} and learning rates {0.001, 0.0005} (however, ExpRNN failed to learn). Lipschitz RNN (Erichson et al., 2021) For MFCC, we swept hidden sizes {256, 512} and learning rates {0.001, 0.002, 0.0005}. Training was run for 150 epochs, with a single layer model using a batch size of 100. For Raw, we found that LipschitzRNN was too slow to train on a single GPU (requiring a full day for 1 epoch of training alone).
WaveGAN Discriminator (Donahue et al., 2019) The WaveGAN-D in Table 4 is actually our improved version of the discriminator network from the recent WaveGAN model for speech (Donahue et al., 2019). This CNN actually did not work well out-of-the-box, and we added several features to help it perform better. The final model is highly specialized compared to our model, and includes:
• Downsampling or pooling between layers, induced by strided convolutions, that decrease the sequence length between layers.
• A global fully-connected output layer; thus the model only works for one input sequence length and does not work on MFCC features or the frequency-shift setting in Table 4.
• It only works with Batch Normalization turned on, whereas S4 works equally well with either Batch Normalization or Layer Normalization.
• Almost 90× as many parameters as the S4 model (26.3M vs. 0.3M)

heading: D.3 GENERAL SEQUENCE MODELING
text: This subsection corresponds to the experiments in Section 4.3. Because of the number of experiments in this section, we use subsubsection dividers for different tasks to make it easier to follow: CIFAR-10 density estimation Appendix D.3.1, WikiText-103 language modeling Appendix D.3.2, autoregressive generation Appendix D.3.3, sequential image classification Appendix D.3.4, and time-series forecasting Appendix D.3.5.

heading: D.3.1 CIFAR DENSITY ESTIMATION
text: This task used a different backbone than the rest of our experiments. We used blocks of alternating S4 layers and position-wise feed-forward layers (in the style of Transformer blocks). Each feedforward intermediate dimension was set to 2× the hidden size of the incoming S4 layer. Similar to Salimans et al. (2017), we used a UNet-style backbone consisting of B identical blocks followed by a downsampling layer. The downsampling rates were 3, 4, 4 (the 3 chosen because the sequence consists of RGB pixels). The base model had B = 8 with starting hidden dimension 128, while the large model had B = 16 with starting hidden dimension 192.
We experimented with both the mixture of logistics from (Salimans et al., 2017) as well as a simpler 256-way categorical loss. We found they were pretty close and ended up using the simpler softmax loss along with using input embeddings.
We used the LAMB optimizer with learning rate 0.005. The base model had no dropout, while the large model had dropout 0.1 before the linear layers inside the S4 and FF blocks.

heading: D.3.2 WIKITEXT-103 LANGUAGE MODELING
text: The RNN baselines included in Table 7 are the AWD-QRNN (Merity et al., 2018), an efficient linear gated RNN, and the LSTM + Cache + Hebbian + MbPA (Rae et al., 2018), the best performing pure RNN in the literature. The CNN baselines are the CNN with GLU activations (Dauphin et al., 2017), the TrellisNet (Bai et al., 2019), Dynamic Convolutions (Wu et al., 2019), and TaLK Convolutions (Lioutas & Guo, 2020).
The Transformer baseline is (Baevski & Auli, 2018), which uses Adaptive Inputs with a tied Adaptive Softmax. This model is a standard high-performing Transformer baseline on this benchmark, used for example by Lioutas & Guo (2020) and many more.
Our S4 model uses the same Transformer backbone as in (Baevski & Auli, 2018). The model consists of 16 blocks of S4 layers alternated with position-wise feedforward layers, with a feature dimension of 1024. Because our S4 layer has around 1/4 the number of parameters as a self-attention layer with the same dimension, we made two modifications to match the parameter count better: (i) we used a GLU activation after the S4 linear layer (Section 3.4) (ii) we used two S4 layers per block. Blocks use Layer Normalization in the pre-norm position. The embedding and softmax layers were the Adaptive Embedding from (Baevski & Auli, 2018) with standard cutoffs 20000, 40000, 200000.
Evaluation was performed similarly to the basic setting in (Baevski & Auli, 2018), Table 5, which involves sliding non-overlapping windows of width 1024 tokens. Other settings are reported in (Baevski & Auli, 2018) that include more context at training and evaluation time and improves the score. Because such evaluation protocols are orthogonal to the basic model, we do not consider them and report the base score from (Baevski & Auli, 2018) Table 5.
Instead of SGD+Momentum with multiple cosine learning rate annealing cycles, our S4 model was trained with the simpler AdamW optimizer with a single cosine learning rate cycle with a maximum of 800000 steps. The initial learning rate was set to 0.0005. We used 8 A100 GPUs with a batch size of 8 per gpu and context size 1024. We used no gradient clipping and a weight decay of 0.1. Unlike (Baevski & Auli, 2018) which specified different dropout rates for different parameters, we used a constant dropout rate of 0.25 throughout the network, including before every linear layer and on the residual branches.

heading: D.3.3 AUTOREGRESSIVE GENERATION SPEED
text: Protocol. To account for different model sizes and memory requirements for each method, we benchmark generation speed by throughput. measured in images per second (Table 6) or tokens per second (Table 7). Each model generates images on a single A100 GPU, maximizing batch size to fit in memory. (For CIFAR-10 generation we limited memory to 16Gb, to be more comparable to the Transformer and Linear Transformer results reported from (Katharopoulos et al., 2020).)
Baselines. The Transformer and Linear Transformer baselines reported in Table 6 are the results reported directly from Katharopoulos et al. (2020). Note that the Transformer number is the one in their Appendix, which implements the optimized cached implementation of self-attention.
For all other baseline models, we used open source implementations of the models to benchmark generation speed. For the PixelCNN++, we used the fast cached version by Ramachandran et al. (2017), which sped up generation by orders of magnitude from the naive implementation. This code was only available in TensorFlow, which may have slight differences compared to the rest of the baselines which were implemented in PyTorch.
We were unable to run the Sparse Transformer (Child et al., 2019) model due to issues with their custom CUDA implementation of the sparse attention kernel, which we were unable to resolve.
The Transformer baseline from Table 7 was run using a modified GPT-2 baackbone from the HuggingFace repository, configured to recreate the architecture reported in (Baevski & Auli, 2018). These numbers are actually slightly favorable to the baseline, as we did not include the timing of the embedding or softmax layers, whereas the number reported for S4 is the full model.

heading: D.3.4 PIXEL-LEVEL SEQUENTIAL IMAGE CLASSIFICATION
text: Our models were trained with the AdamW optimizer for up to 200 epochs. Hyperparameters for the CIFAR-10 model is reported in Table 9.
For our comparisons against ResNet-18, the main differences between the base models are that S4 uses LayerNorm by default while ResNet uses BatchNorm. The last ablation in Section 4.3 swaps the normalization type, using BatchNorm for S4 and LayerNorm for ResNet, to ablate this architectural difference. The experiments with augmentation take the base model and train with mild data augmentation: horizontal flips and random crops (with symmetric padding).

heading: D.3.5 TIME SERIES FORECASTING COMPARED TO INFORMER
text: We include a simple figure (Appendix D.3.5) contrasting the architecture of S4 against that of the Informer (Zhou et al., 2021).
In Appendix D.3.5, the goal is to forecast a contiguous range of future predictions (Green, length F ) given a range of past context (Blue, length C ). We simply concatenate the entire context with a sequence of masks set to the length of the forecast window. This input is a single sequence of length C + F that is run through the same simple deep S4 model used throughout this work, which maps to an output of length C + F . We then use just the last F features as the forecasted predictions.
Tables 11 and 12 contain full results on all 50 settings considered by Zhou et al. (2021). S4 sets the best results on 40 out of 50 of these settings.
D.4 VISUALIZATIONS
We visualize the convolutional filter K̄ learned by S4 for the Pathfinder and CIFAR-10 tasks in Appendix D.4 and Appendix D.4.

Review summaries:
This paper presents a novel parameterization for the established state-space models (SSM), which tackles the scalability problem in linear state-space layers (LSSL) in modeling long-range dependencies for very long sequences. The proposed technique reparameterizes the structured state matrix in LSSL that allows the state to memorize the past, a key component following the continuous-time memorization theory. A complete theoretical description of this reparameterization and empirical evaluations on a diverse set of benchmarks are presented. The proposed model, S3, achieves astonishing results both in terms of performance and computational efficiency. 

The paper extends previous work on linear state space models (SSM), where the state transition matrix is fixed to be a highly structured _HiPPO_ matrix, which has provably beneficial properties for memorizing long-term information from continuous-time signals. The main contribution of the paper is regarding the computational aspect of the model. Namely, a novel approach is proposed for computing the convolutional kernel associated to the discretized SSM unrolled over time. This is done by first showing that the state transition matrix can be decomposed as the sum of a normal (i.e. diagonalizable with orthonormal eigenvectors) and a low-rank (rank 1 or 2) matrix, and this representation is then combined with techniques from numerical linear algebra to reduce the problem to computing the Cauchy kernel (i.e. a well-studied problem). The resulting structured state space model (S3) is then placed into a deep neural network setting, and extensive experiments are carried out on various tasks, such as: 1) Long Range Arena, a benchmark collection for scalable transformers, 2) raw speech classification (length-16k audio signals), 3) generative modelling on CIFAR-10 and WikiText-103, 4) sequential image classification on sMNIST, pMNIST and sCIFAR. Overall, the model seems to perform very well on each task either performing close to SoTA or setting a new high score.

The authors propose a sequence modeling approach called the structured state space model (S3) which parameterized the SSM in a more computationally efficient manner.  This is done through decomposing the structured state matrices A into a low-rank and skew-symmetric term and expanding the SSM in frequency space and using a multipole-like evaluation.  This approach maintains much of the gains and efficiencies of past SSM approaches while being more computationally stable and efficient as demonstrated on several tasks across broad domains (speech, images, text).  The paper first steps through the theoretical motivations and justifications of this approach and then performs a number of experiments to demonstrate the competitive or superior performance on a wide range of LRD tasks.