Efficiently Modeling Long Sequences with Structured State Spaces

A sequence model that reaches 16,000 steps back while carrying 64 numbers.

Underneath the layer is a linear differential equation borrowed from control theory. Making one fast enough to train came down to the shape of a single matrix.

The Long Range Arena benchmark shipped with a task nobody could do. Eleven Transformer variants were measured on it and every one scored exactly what you get by guessing.

The task is called Path-X. You get a 128×128 picture containing two small circles and a tangle of dashed curves, and you answer one question: are the two circles joined by a path? A person solves it by tracing. A convolutional network solves the 32×32 version easily. The catch is that the model is not handed a picture. It is handed the pixels one at a time, in reading order, all 16,384 of them, and it never gets to see them laid out in two dimensions.

That flattening is what makes the task brutal, and you can see why before any of the machinery. Two pixels sitting directly above and below each other in the picture arrive 128 tokens apart in the sequence. The two circles you have to connect can be ten thousand tokens apart. Drag the teal circle around and watch the gap on the strip below; a single step downward in the picture costs a full row of tokens.

S4 solved it, at 96.35%. The paper that did it is not really an architecture paper. It takes a model that control engineers have used since the 1960s, points out that nobody could run it fast enough to be a deep learning layer, and spends its length fixing that. The fix is a piece of numerical linear algebra, and it is the reason this line of work exists: Mamba, Mamba-2, and the rest of the state space family descend from the algorithm in this paper.

The argument runs in a few steps: what the model is, what step size means, which matrix to use, why that matrix cannot be handled the obvious way, and the trick that handles it anyway.

Why length is the difficulty

Before the model, be precise about the enemy. If you want to relate a token to something 16,000 steps earlier, you have exactly two options, and both had walls.

Attention compares every position with every other position. That gives any token direct access to any other, at any distance, which is why Transformers are good at range. It costs $L^2$ comparisons. At $L = 16{,}384$ that is 268 million pairwise scores per head per layer, and a score matrix of the same size to hold. (FlashAttention later removed the memory half of that problem by never materializing the matrix; the 268 million multiply-adds stay.) The eleven efficient-Transformer variants in the benchmark all attack the same quadratic term with sparsity or low-rank approximations, and all of them still failed Path-X, some because they ran out of memory and some because they simply did not learn it.

A recurrent network has the opposite shape. It carries a hidden state forward and pays a constant amount per step, so length costs nothing in memory. What it loses is the signal: to influence step 16,000, an input at step 1 has to survive 16,000 multiplications by the recurrent weight matrix, and the product of 16,000 matrices either shrinks to nothing or explodes. That is the vanishing gradient problem, and it is why LSTMs top out in the low thousands of steps.

So: one family with perfect range and quadratic cost, one family with linear cost and no range. S4 is a third thing that gets linear-ish cost and real range, and the way it does that is by being a continuous-time system underneath.

A linear system as a layer

The state space model is one equation, and it predates deep learning by decades:

Read it as a machine with a dial panel. $u(t)$ is a single incoming number that changes over time, a scalar signal. $x(t)$ is a vector of $N$ numbers, the state, which is the machine's entire memory of everything that has happened. $y(t)$ is a single outgoing number. The first equation says how fast each of the $N$ state numbers is changing right now: partly from the other state numbers (through $\bm{A}$), partly from whatever just arrived (through $\bm{B}$). The second says the output is a fixed weighted sum of the state (through $\bm{C}$), plus a direct copy of the input scaled by $\bm{D}$.

Pin down the word "state" before anything else. In control theory a state is a summary of the past that is sufficient: given $x(t)$ and everything that arrives after $t$, you never need to look at the input history again. So those $N$ numbers are not a window, a cache, or a sample. They are a compression of all of history into a fixed budget, and the model can only remember what survives that compression. S4 fixes $N = 64$ throughout the Long Range Arena experiments.

$\bm{D}$ is the least interesting symbol here. It is a direct wire from input to output that skips the state entirely, so the paper drops it from the exposition and the code adds it back as a per-channel skip connection.

The equation has two properties that everything else in the paper hangs off, so they deserve names. It is linear: feed it two inputs added together and you get the two outputs added together, with no gate, no $\tanh$, nothing that bends. And it is time-invariant: $\bm{A},\bm{B},\bm{C}$ do not depend on $t$ or on the input, so a spike fed in at step 10 produces exactly the same response shape as a spike fed in at step 5,000, just shifted. Together those two properties are called LTI, and they force a strong conclusion: the output has to be a convolution.

The argument is short. Any input is a sum of spikes, one per time step, scaled by the input values. Linearity means the response to the sum is the sum of the responses. Time-invariance means the response to each spike is one fixed shape, just shifted to where the spike was. So the output is that one fixed shape, slid along and scaled, added up: a convolution. The fixed shape has a name in signal processing, the impulse response, and knowing it is knowing the whole system.

To get it in a form you can compute, first turn the ODE into a recurrence over discrete steps (the next section says how) and then unroll it from a zero initial state:

The bars mean "discretized". $\bm{\overline{K}}$ is the impulse response written out: feed in a single 1 followed by zeros and the output is $\bm{\overline{K}}$, term by term. Since $\bm{\overline{B}}$ is $N \times 1$ and $\bm{\overline{C}}$ is $1 \times N$, each entry $\bm{\overline{C}}\bm{\overline{A}}^{\,i}\bm{\overline{B}}$ is an ordinary number, so $\bm{\overline{K}}$ is a vector of $L$ numbers, a filter as long as the sequence.

Nothing was approximated in that last step. The recurrence and the convolution are the same function, not two nearby ones, provided the initial state is zero (the paper assumes it, and it is what the code does). That equality is the deal S4 lives on: train with the convolution, because every output can be computed at once, and generate with the recurrence, because it costs a fixed amount per token no matter how long the context has grown. Two schedules for one set of parameters. Switch to the convolutional view below and check that the output numbers are literally the same.

The numbers in that figure are worth doing by hand once, because the rest of the paper is about computing them for large $N$ and large $L$. Take the smallest real version: $N = 2$, $\Delta = 1$, with $\bm{A}$ and $\bm{B}$ from the HiPPO matrix two sections below and a readout that picks off the first state coordinate:

Discretizing at $\Delta = 1$ gives $\bm{\overline{A}} = \left[\begin{smallmatrix} 1/3 & 0 \\ -1/\sqrt{3} & 0\end{smallmatrix}\right]$ and $\bm{\overline{B}} = (2/3,\, 1/\sqrt{3})^{\top}$. Because $\bm{\overline{A}}$ is triangular with a $1/3$ in the corner, the kernel is a plain geometric decay, $\bm{\overline{K}}_i = \tfrac{2}{3}(\tfrac13)^i$, so

Run the pulse $u = (1,1,1,1,0,0,0,0)$ through it and you get $y = (0.667,\, 0.889,\, 0.963,\, 0.988,\, 0.329,\, 0.110,\, 0.037,\, 0.012)$, whether you unroll the recurrence, convolve directly, or use an FFT. The output charges toward 1 while the pulse is on and decays by a factor of three per step afterwards, which is a leaky accumulator, which is exactly what a running average of recent history should look like.

Choosing the step size

Equation (1) has no time steps in it. It describes a signal flowing continuously, and a sequence model has to work on $u_0, u_1, u_2, \dots$. Bridging those needs one number, the step size $\Delta$, which declares how much continuous time passes between two consecutive tokens: $u_k = u(k\Delta)$. Nothing in the data tells you what $\Delta$ should be. It is a parameter, and S4 learns it.

You could discretize with the crudest thing available, forward Euler, $x_k = x_{k-1} + \Delta(\bm{A}x_{k-1} + \bm{B}u_k)$, but forward Euler can turn a stable system unstable when the step is large, and stability is the one property this model must not lose. S4 uses the bilinear transform (Tustin, 1947), which is the trapezoidal rule applied to the ODE:

Read it as half a step forward followed by half a step backward. The reason to prefer it is a fact about where eigenvalues land. Since $\bm{\overline{A}}$ is a rational function of $\bm{A}$, each eigenvalue $\lambda$ of $\bm{A}$ becomes an eigenvalue of $\bm{\overline{A}}$ by the same formula:

That map sends the entire left half of the complex plane onto the inside of the unit circle, exactly, for every $\Delta > 0$. A continuous system is stable when its eigenvalues have negative real part; a discrete one is stable when its eigenvalues are inside the unit circle. So the bilinear transform preserves stability at any step size, while forward Euler does not. You can check the map by hand: $\lambda = -1$ with $\Delta = 0.1$ goes to $0.9048$, and $\lambda=-8$ goes to $0.4286$. Both inside.

What $\Delta$ buys, besides stability, is a horizon. A small $\Delta$ means each token advances the underlying system by very little, so the state decays slowly across tokens and the model looks far back; a large $\Delta$ means each token is a big jump and the past is forgotten quickly. The code initializes one $\Delta$ per channel, drawn log-uniformly from $[0.001,\ 0.1]$, which puts $1/\Delta$ anywhere between 10 and 1000 tokens of modelled context. A single layer is therefore a rack of filters running at different timescales at once.

It also explains one of the paper's odder results. Speech Commands is recorded at 16 kHz; S4 trained at that rate scores 96.30% on test audio resampled to half the rate, with no retraining, by doubling $\Delta$. Halving the sample rate doubles the real time between tokens, so doubling $\Delta$ describes the same underlying continuous system. Nothing about the learned $\bm{A},\bm{B},\bm{C}$ changes. That is a manual rescale rather than automatic invariance, and it works because the parameters describe a continuous object in the first place.

Which A? The HiPPO matrix

Everything so far works for any $\bm{A}$. It also fails for almost any $\bm{A}$. Take a random Gaussian matrix, scale it down until training does not produce NaNs, and the model learns badly: the LSSL, the paper's immediate predecessor, measured about 62% on permuted MNIST with a random $\bm{A}$ against 98.76% with the right one.

The right one comes from HiPPO (Gu et al., 2020), a paper about a different question: if you can only keep $N$ numbers, what is the best possible summary of a signal's history, and can you maintain it online? HiPPO's answer is a running least-squares polynomial fit. Keep the $N$ coefficients of the best degree-$N$ polynomial through the history so far, in an orthogonal basis. The result they proved is that those coefficients obey a linear ODE, so you never re-fit anything. You integrate. The matrix that does it, for the scaled Legendre basis, is:

It looks arbitrary and it is not. The square roots are the normalization of the Legendre polynomials, the lower-triangular structure says coefficient $n$ is fed by all the coarser coefficients below it, and the diagonal $-(n+1)$ is the rate at which the $n$-th coefficient forgets. Higher coefficients carry finer detail and decay faster, which is what you would design by hand if you had to.

S4 makes one change to HiPPO's system, and it changes what the state remembers. HiPPO's original system is time-varying: the exact scaled-Legendre operator is $\tfrac{d}{dt}x = \tfrac1t\bm{A}x + \tfrac1t\bm{B}u$, with a $1/t$ that widens the window as time passes. S4 keeps the matrix and drops the $1/t$, running it at a fixed $\Delta$ so the system stays time-invariant, because time-invariance is what makes the convolution exist. And that substitution changes the memory. Run equation (2) as a fixed linear system and its impulse response has a closed form,

which is the Legendre polynomial $P_n$ evaluated not at time but at $e^{-t}$. So S4 fits a polynomial to the past in a warped variable: recent history is stretched out and gets most of the resolution, distant history is squeezed toward zero. The memory decays smoothly with no hard cutoff, rather than averaging uniformly over everything since the start. The paper does not flag the swap; it matters because it is the difference between "remembers all of history equally" and "remembers recent history well and old history vaguely", and the second is what you get.

The figure below is that statement drawn. The horizontal axis is the warped variable, so "now" sits at the right edge and age runs left, compressed. Slide $N$ up from 1 and watch which parts of the history sharpen first.

This initialization holds up over long sequences partly because the system is contractive with no transient. The symmetric part of equation (2) works out to $\tfrac12(\bm{A}+\bm{A}^{\top}) = -\tfrac12\bm{I}-\tfrac12\bm{P}\bm{P}^{\top}$ with $\bm{P}_n=(2n+1)^{1/2}$, which is at most $-\tfrac12\bm{I}$, so $\lVert e^{t\bm{A}}\rVert_2 \le e^{-t/2}$ for every $t$. The state shrinks monotonically instead of swelling first and shrinking later, which is the failure mode that wrecks badly conditioned recurrences. Hold onto that $\bm{P}$; it is about to do the paper's real work.

Why you cannot just diagonalize it

Now the problem. Building the kernel $\bm{\overline{K}}$ from equation (5) means computing $\bm{\overline{C}}\bm{\overline{A}}^{\,i}\bm{\overline{B}}$ for $i$ up to $L-1$, which naively means $L$ successive multiplications by an $N \times N$ matrix: $O(N^2L)$ operations and $O(NL)$ space to hold all the intermediate states. At $N = 256$, the size the LSSL used, that is orders of magnitude more memory than a comparable RNN, which is why the LSSL could not be run at the size it wanted.

The textbook escape is diagonalization. If $\bm{A} = \bm{V}\bm{\Lambda}\bm{V}^{-1}$ with $\bm{\Lambda}$ diagonal, then $\bm{A}^i = \bm{V}\bm{\Lambda}^i\bm{V}^{-1}$, and powering a diagonal matrix is just raising $N$ numbers to the $i$-th power. And you are allowed to change basis: the paper's Lemma 3.1 notes that $(\bm{A},\bm{B},\bm{C})$ and $(\bm{V}^{-1}\bm{A}\bm{V},\, \bm{V}^{-1}\bm{B},\, \bm{C}\bm{V})$ compute the identical map from $u$ to $y$, because the change of basis cancels. In exact arithmetic, conjugating is free.

Equation (2) is even friendly to the idea. It is lower triangular, so its eigenvalues are sitting on the diagonal in plain sight: exactly $-1, -2, \dots, -N$, all real, all distinct. Distinct eigenvalues guarantee a full set of eigenvectors. The matrix is diagonalizable, and you can write down both the eigenvalues and the eigenvectors in closed form.

Then you try it in floating point and it disintegrates. The eigenvectors, though independent, are nearly parallel, so the matrix $\bm{V}$ that holds them is catastrophically ill-conditioned: expressing an ordinary vector in that basis takes enormous coordinates that then almost cancel, and a float64 carries only about 16 significant digits, so the cancellation eats the answer. The condition number of $\bm{V}$ is about $7.7\times10^4$ at $N=8$, $8.3\times10^{10}$ at $N=16$, and (computed in 260-digit arithmetic, because float64 cannot even measure it) $1.2\times10^{23}$ at $N=32$ and $3.2\times10^{47}$ at $N=64$, the size S4 uses. In bits: you would need 158 of them at $N=64$, against float64's 53.

The paper puts it as a bound on the size of the entries of $\bm{V}$, $2^{4N/3}$, which at $N=64$ is around $10^{25}$. Either way the conclusion is the same, and the failure deserves a precise name. Nothing overflows, and the matrix is not defective. The eigenvalues stay exact at every size. The eigenvectors are what cannot be represented, so it is a loss of relative precision, not a breakdown of the algebra.

The figure demonstrates that rather than asserting it. The left panel is equation (2); the right panel is the same matrix plus the rank-one correction the next section builds. The meter underneath takes the raw matrix apart into its exactly-known eigenvalues and eigenvectors, multiplies them back together in double precision, and reports how far the result is from where it started. Walk $N$ up and watch the answer go from correct to meaningless somewhere in the low twenties.

The class of matrices where diagonalization is safe has a name. A matrix is normal when it commutes with its own transpose, and the spectral theorem says exactly the normal matrices are the ones diagonalizable by a unitary $\bm{V}$, meaning a rigid rotation, meaning $\text{cond}(\bm{V}) = 1$ and nothing is lost. Symmetric matrices are normal. Skew-symmetric matrices are normal. Equation (2) is not, and that is the obstruction.

Normal plus rank one

Demanding that $\bm{A}$ be normal would be a real restriction, and it would throw away HiPPO. The paper instead asks for something weaker: let $\bm{A}$ be a normal matrix plus a low-rank correction. That class is strictly bigger, and it contains the matrix we want. The paper calls it NPLR, normal plus low-rank, and after the unitary change of basis that diagonalizes the normal part, DPLR, diagonal plus low-rank.

For HiPPO the correction is rank one, and it is not an approximation or a fitted decomposition. It is an identity you can verify with algebra on a general entry. Take $\bm{P}_n = (2n+1)^{1/2}$ and add $\tfrac12\bm{P}\bm{P}^{\top}$ to equation (2). Above the diagonal, where $\bm{A}$ is zero, you get $+\tfrac12(2n+1)^{1/2}(2k+1)^{1/2}$. Below it you get $-(2n+1)^{1/2}(2k+1)^{1/2} + \tfrac12(2n+1)^{1/2}(2k+1)^{1/2}$, which is the same number with the opposite sign. On the diagonal you get $-(n+1) + \tfrac12(2n+1) = -\tfrac12$, the same value in every row. So:

A skew-symmetric matrix plus a multiple of the identity is normal, so the corrected matrix is unitarily diagonalizable, and its spectrum is pinned: every eigenvalue is $-\tfrac12 + i\omega$ for some real $\omega$. At $N=8$ the eigenvalues are $-0.5 \pm 19.857i$, $-0.5 \pm 5.354i$, $-0.5 \pm 1.958i$, $-0.5 \pm 0.428i$. Every mode decays at the same rate and they differ only in how fast they oscillate, which is the Legendre structure showing up as frequency. This is the object S4 stores and trains: a diagonal $\bm{\Lambda}$ plus one outer product.

So the parameters of an S4 layer, per channel, are $\bm{\Lambda}, \bm{P}, \bm{Q}, \bm{B}, \bm{C}$, five complex vectors of length $N$, which is the paper's $5N$. In the released code it is a little different, and the differences matter if you go read it: $\bm{Q}$ is tied to $\bm{P}$, only half of each conjugate pair is stored, and $\Delta$ and $\bm{D}$ are trained too, which the $5N$ count does not include. The kernel is never a parameter. At Path-X scale, roughly 66,000 numbers generate a $256 \times 16{,}384$ filter bank, 4.2 million taps, on every forward pass.

Powers out, inverses in

Writing $\bm{A} = \bm{\Lambda} - \bm{P}\bm{Q}^{*}$ does not by itself make anything fast, and the reason explains every move that follows. Powering a diagonal-plus-rank-one matrix does not stay diagonal-plus-rank-one. Multiply it out: $(\bm{\Lambda}-\bm{P}\bm{Q}^{*})^2$ has a $\bm{\Lambda}^2$ term and three cross terms, and the correction now needs two rank-one pieces to describe. Cube it and you need three. The rank of the correction climbs by one with every multiplication, so after about $N$ steps the matrix is structurally dense and you are back to $O(N^2L)$.

Inverting is the opposite. The Woodbury identity says a rank-$p$ update of a matrix produces a rank-$p$ update of its inverse, and the only new inverse you need is $p \times p$:

At $p=1$ that inner inverse is a single number, so a rank-one-corrected inverse costs one division. Structure survives inversion and dies under powering. So the move is to find a way to write the kernel using inverses instead of powers.

That way is a generating function. Pack the $L$ kernel values into the coefficients of a polynomial and evaluate it, rather than building the coefficients directly:

This is the finite geometric series $1 + r + \dots + r^{L-1} = (1-r^L)/(1-r)$ with a matrix in place of $r$, and it is a polynomial identity, true for any $z$ where the inverse exists. Look at what happened: the powers went away and one inverse arrived.

Two more steps make it computable. First, evaluate at the $L$-th roots of unity, $z_k = e^{-2\pi i k/L}$. That choice does two jobs at once. It makes the values of $\hat{\mathcal{K}}_L$ literally the discrete Fourier transform of $\bm{\overline{K}}$, so one inverse FFT recovers the kernel exactly and no information is lost. And it kills the awkward term: at a root of unity $z^L = 1$, so $(\bm{I}-\bm{\overline{A}}^{L}z^{L})$ stops depending on $z$ and can be folded into $\bm{\overline{C}}$ once and for all as $\bm{\tilde{C}} = (\bm{I}-\bm{\overline{A}}^{L})^{*}\bm{C}$. S4 simply makes $\bm{\tilde{C}}$ the learned parameter and never computes that factor at all. Away from the unit circle the shortcut is wrong, so the roots of unity are load-bearing rather than merely convenient.

Second, unfold the remaining inverse with Woodbury. Writing $s = \tfrac{2}{\Delta}\tfrac{1-z}{1+z}$ and the resolvent $\bm{R}(z) = (s - \bm{\Lambda})^{-1}$, the evaluation collapses to:

Every one of those four bracketed quantities has the same shape. Since $\bm{\Lambda}$ is diagonal, $\bm{v}^{*}\bm{R}\bm{w}$ is nothing but $\sum_j \overline{v_j}w_j/(s-\lambda_j)$, a sum of $N$ divisions. A matrix of entries $1/(s_i - \lambda_j)$ is called a Cauchy matrix, and multiplying by one is a classical problem in numerical analysis with fast and numerically stable algorithms behind it. That is Theorem 3: the kernel reduces to four Cauchy multiplies, $\tilde{O}(N+L)$ operations and $O(N+L)$ space, down from $O(N^2L)$ and $O(NL)$.

Two details make that sum numerically safe as well as compact. The node $s$ from above is the bilinear image of a point on the unit circle, and for $z = e^{-i\theta}$ it works out to $s = \tfrac{2i}{\Delta}\tan(\theta/2)$, which is purely imaginary. The poles $\lambda_j$ all sit at real part $-\tfrac12$. So $\lvert s - \lambda_j\rvert \ge \tfrac12$ always, and no division in the sum is ever near a blow-up. The rank-one correction bought numerical safety twice over.

Here is the algorithm, in the shape you would actually write it:

# Algorithm 1: the length-L kernel, straight from N poles.
# Lam, P, B, Ct are complex vectors of length N. Ct is the learned C~.

z = np.exp(-2j * np.pi * np.arange(L) / L)     # the L-th roots of unity
s = (2 / dt) * (1 - z) / (1 + z)               # bilinear image: imaginary

def cauchy(v):                                 # one divide per pole per node
    return (v[None, :] / (s[:, None] - Lam[None, :])).sum(-1)

k00 = cauchy(Ct.conj() * B)                    # four sums, that is all
k01 = cauchy(Ct.conj() * P)
k10 = cauchy(P.conj() * B)                     # released code ties Q = P
k11 = cauchy(P.conj() * P)

Khat = 2 / (1 + z) * (k00 - k01 / (1 + k11) * k10)   # Woodbury: a divide
K = np.fft.ifft(Khat).real                     # back to the time domain

Two footguns in that listing will catch a reader implementing from the paper. The roots of unity carry a minus sign, matching Lemma C.2 and the released code; Algorithm 1 prints a plus, and pairing that with its own inverse FFT returns the kernel time-reversed, an anti-causal filter. And when $L$ is even, one node lands on $z=-1$, where the $2/(1+z)$ factor blows up. That is a removable singularity: the resolvent vanishes at exactly the same rate and the limit is $\bm{\tilde{C}}^{*}(\bm{I}+\bm{\overline{A}})^{-1}\bm{\overline{B}}$. Keep the node.

The headline complexity comes with a caveat. The $\tilde{O}(N+L)$ in Theorem 3 needs a fast multipole implementation of the Cauchy multiply, and the paper's own Appendix A says that was never written: the release computes the sum the naive way, $O(NL)$ time, in $O(N+L)$ space because the $N\times L$ matrix is never materialized. It costs little because the kernel is built once per batch rather than once per sample, but the shipped constant is not the theorem's.

The figure below runs exactly that listing, on the real HiPPO initialization at $N=64$, in your browser. The left panel is the kernel it produces; the right panel is where the discretized poles sit inside the unit circle. Drag $\Delta$ across its full initialization range and past both ends. At the bottom the filter stretches out over more than a thousand taps and the poles crowd into a narrow fan near $z=1$; at the top the fan opens all the way round, the fastest modes alias, and the filter turns into a short jagged burst. The poles never leave the circle, at any $\Delta$, which is the bilinear transform's promise from two sections ago made visible.

The layer, and what it costs

One SSM maps a scalar sequence to a scalar sequence. A network needs to move $H$ features. S4 does the obvious thing: run $H$ independent copies, one per channel, with no interaction, then mix the channels with a position-wise linear layer, then apply a nonlinearity. That structure has a name already. It is a depthwise-separable convolution, one filter per channel followed by a $1\times1$ mix, and the code literally builds the mixer as a kernel-size-1 convolution.

The analogy is exact in shape and breaks in three useful places. The kernels are as long as the entire sequence rather than a small window. They are generated from about $4N$ to $5N$ numbers per channel rather than stored, so the parameter count does not grow with $L$ and you can re-materialize the same filter at a different length or a different $\Delta$ at test time. And they have a recurrent twin, which no stored CNN kernel has.

Note also what is not in the layer. There is no gate anywhere in the state update. The recurrence $x_k = \bm{\overline{A}}x_{k-1} + \bm{\overline{B}}u_k$ is strictly linear and its parameters do not depend on the input, which is the property that makes the convolutional view exist at all. All the nonlinearity in a deep S4 model lives between layers. (Modern S4 blocks do put a GLU after the channel mixer, but that is downstream of the SSM, not inside it.) Hold onto that; it is the exact property Mamba gives up on purpose.

The convolution is done with FFTs, and an FFT computes circular convolution: the tail of the answer wraps around and lands on top of the beginning. Take $K = (1,2)$ and $u = (3,4)$. The linear convolution is $(3,10,8)$; the circular one of length 2 is $(11,10)$, because the 8 came back around and landed on the 3. The cure is zero-padding to at least $2L-1$ before transforming, and truncating afterwards. The code pads to $2L$ for tidiness. Skip the pad and the output is not slightly noisy, it is wrong by an amount the size of the signal.

# training: build the kernel once, then one FFT for all L outputs
K = kernel(Lam, P, B, Ct, dt, L)          # (N + L) numbers touched
y = irfft(rfft(u, 2*L) * rfft(K, 2*L), 2*L)[:L] + D * u

# generation: the same parameters, run as a recurrence
x = zeros(N)                              # 64 numbers, and that is the
for u_k in stream:                        # entire memory of the past
    x = Abar @ x + Bbar * u_k             # O(N), no matter how long
    y_k = C.conj() @ x                    # C, not Ct: solve (I - Abar^L)^H C = Ct once

Now the cost. Attention grows with the square of the length. The naive SSM already grew only linearly, so linearity in $L$ was never S4's contribution; what S4 fixed was the $N^2$ in front and the $NL$ of state the naive version had to keep. Slide the length across the figure and watch the three curves separate, then switch to memory and watch the naive SSM's problem appear.

Two qualifications about that chart. S4 is $O(L\log L)$ in training compute, not $O(L)$, because the FFT is in there; only memory and per-step generation are truly linear and constant. And at realistic widths the position-wise mixer, at $\Theta(BLH^2)$, is the larger of S4's two terms and is missing from the paper's own complexity table. At $B{=}32$, $H{=}256$, $L{=}16{,}384$ the FFT convolution is about $1.2\times10^{10}$ multiply-adds and the mixer about $3.4\times10^{10}$, against $2.2\times10^{12}$ for attention. Both terms together are about 2% of attention's count.

The measured version matches. Against the LSSL, at layer width 512, S4 is 29.6× faster per training step and uses 392× less memory. Against the efficient Transformers, parameter-matched at length 4096, S4 runs 5.19× faster than a vanilla Transformer on 9.1% of the memory, which puts it beside the Performer and the Linear Transformer rather than paying a premium for its extra range.

What it bought

On Long Range Arena, S4 leads every baseline on every one of the six tasks and averages 86.09% against under 60% for all of them. The telling number is the Path-X column, where the baselines are not merely behind; they are at 50.00, the value you get by flipping a coin. S4 is at 96.35.

A note on that 86.09, because the paper is inconsistent with itself. The prose in Section 4.2 still quotes 80.48%, which is the superseded row from the first version; the table beside it reports 86.09, and the appendix prints both rows. Appendix D.5 explains that the results were refreshed from the S4D and "How to Train Your HiPPO" follow-ups, same model, minor hyperparameter changes. The one Path-X-specific change behind the move from 88.10 to 96.35 is the one you would guess from the step-size section: the $\Delta$ initialization range was dropped from $[0.001,0.1]$ to $[0.0001,0.01]$, a ten-fold longer modelled horizon for a sixteen-thousand-token input. If you cite the table, cite 86.09.

The rest of the results are a spread rather than a single benchmark, which is the point the paper is making about generality. On raw speech, classifying 16,000-sample waveforms with no feature extraction, S4 gets 98.32%, beating every baseline that was handed the 100× shorter MFCC features instead, and beating a CNN discriminator purpose-built for audio that has 90× more parameters (26.3M against 0.3M). On sequential CIFAR-10, one pixel at a time with no 2-D structure available, S4 reaches 91.13%, which is competitive with a 2-D ResNet-18 half again as big: 7.9M parameters for S4 against 11.0M. On CIFAR-10 density estimation it hits 2.85 bits per dimension with no 2-D inductive bias, matching PixelSNAIL, which has 2-D convolutions and attention.

On WikiText-103 it does not win, and the paper says so. Swapping S4 for the attention layers in a standard Transformer language model gives 20.95 test perplexity at 249M parameters against 20.51 for the 247M Transformer. (The prose says "within 0.8 ppl", which was calibrated to the original 21.28 run; with the retrained number the gap is 0.44.) It is the best attention-free result by more than 2 points, and it is still behind.

What it does win on there is generation speed, by 60×: 48,000 tokens per second against 800. That comes straight from the recurrent view. A Transformer generating token $k$ attends over all $k$ previous tokens, so each new token costs more than the last; S4 carries 64 numbers and pays the same for every token forever. That 60× is throughput at maximum batch size on one A100 rather than single-sequence latency, and the speedup ranges from 65× down to 10.5× for the larger CIFAR model, but the mechanism behind it is the constant-cost step.

The ablation, and what came after

The most interesting experiment in the paper is the one that partly deflates it. Section 4.4 asks which piece is doing the work, and the answer is not the piece the title is named after. Random matrices put into S4's NPLR parameterization still perform badly. Training $\bm{A}$ helps every initialization, and with $\bm{A}$ trainable all of them reach perfect training accuracy, but their validation accuracies stay more than 15 points apart. Same capacity, same optimizer, same parameter count, different starting matrix, and the gap is a generalization gap.

So the HiPPO initialization is the payload and the NPLR machinery is the vehicle that makes it affordable. In that ablation grid the full method tops out at 84.27% on sequential CIFAR-10 at 100K parameters. Publishing the ablation at all was unusually candid, and it set the agenda for everything after: if the initialization is what matters, the parameterization can be simplified.

It was, almost immediately. DSS (Gupta et al., 2022) and then S4D showed that you can throw away the low-rank term entirely and keep only a diagonal $\bm{\Lambda}$, initialized in the spirit of HiPPO, and land within a point or so on average. That deletes the Woodbury correction, the Cauchy structure, and most of this paper's Appendix C, replacing them with a Vandermonde product. The tail is not free: S4D-LegS averages 84.89 against S4's 86.09, and the popular linear-spaced diagonal initialization fails Path-X outright. Both gaps land on Path-X.

The bigger departure went the other way. Everything in this paper rests on the model being time-invariant: $\bm{A},\bm{B},\bm{C}$ are the same at every position, which is what lets one fixed kernel represent the entire layer and one FFT compute it. That is also a real limitation. A filter that treats every position identically cannot decide to ignore a token or hold onto one, which is why S4 trails on language, where content matters more than distance. Mamba makes $\bm{B}, \bm{C}$ and $\Delta$ functions of the input, which buys exactly that selectivity and costs exactly the convolution: with parameters that change per token there is no single kernel left to FFT, so Mamba replaces it with a hardware-aware parallel scan. Mamba-2 then restricts the state matrix further so the scan becomes a matrix multiplication a GPU likes.

Read backwards from there, S4's contribution is sharper than "a fast sequence model". It established that a linear, time-invariant system with a carefully chosen state matrix can hold sixteen thousand steps of context in sixty-four numbers, and it supplied the numerical machinery that made such a system trainable at all. The Path-X column is the proof that the range is real. Everything since has been an argument about what to give up in exchange for selectivity.