Before we begin

HiPPO emphasized the importance of the structure of matrix $A$ for establishing long-term range in an SSM, and LSSL generalized the construction across all continuously defined $A$. In the previous post, I explained LSSL as one of the foundations of Mamba modeling and included brief derivations of the relevant equations. To be honest, I still do not think I fully understand every one of these fundamentals, and trying to understand Mamba through this route can feel excessively roundabout. That may sound unconfident, but ultimately I do not believe that understanding Mamba modeling requires internalizing state-space modeling all the way down to the bone. My personal view is that top-down, rather than bottom-up, is the right direction.

Many blog posts also explain Mamba clearly in one pass with the help of visualizations. Even so, the reason I chose to write a long review is the same one I gave in the previous post: the attempt to understand Mamba’s foundations has value of its own. If we understand Mamba only at a surface level, we cannot explain precisely why it can address problems in the Transformer architecture—and we should not assume that an architecture is unconditionally good merely because it is new.

What comes easily is lost just as easily. A sound understanding of Mamba’s characteristics should reveal its fundamental advantages and disadvantages, which can then provide a foundation for future research and model development. This post covers S4 modeling, which seeks to resolve LSSL’s computational cost and numerical instability.

Returning to sequence modeling

The central questions in sequence modeling are how efficiently a model can process sequences of a useful length, and whether the amount of accessible context grows with the sequence length. A Transformer sacrifices computation to obtain attention information over the full sequence in a single operation and can extend its reach by increasing the token-embedding length, but ultimately runs into computational limits. RNNs and CNNs have their own characteristics, yet no modeling approach is universal because each comes with trade-offs. Papers on LMU, HiPPO, and related methods proposed state-space modeling (SSM) in response.

To summarize SSM

In one sentence, SSM seeks to construct a latent space with a Linear Time-Invariant system. The matrices in an LTI system’s differential equation are time invariant. At every observation unit obtained by discretizing a continuous signal, the differential equation therefore retains the same form and can create the same latent space $x(t)$ regardless of sequence length. Unlike an RNN that depends on gates, if matrix $A$ has a particular structure, it can extract a long-range-dependent latent—the feature vector most important to a deep-learning model’s predictions. Rather than computing attention token by token with a Transformer, or forming an implicit latent from a stream of data through an RNN, we can predict a genuinely time-invariant latent space. No matter how much time passes—even when the observation range has moved far from its beginning—the hidden space is implemented by the same function, thereby supporting long-range dependency. An SSM can consequently process longer sequences than a Transformer with a fixed attention context and can theoretically handle inputs of unlimited length.

Let us use the paper’s figure to build intuition. In the “Continuous State Space” section on the left, the SSM input $u(t)$ produces an output $y(t)$ modeled through the state of a linear system composed of matrices $A, B, C,$ and $D$. That state is called the hidden state, or latent state $x(t)$. If $x(t)$ acts as memory that continuously stores the input-to-output transformation $\mathbf{u} \rightarrow \mathbf{y}$, then fixed, theoretically invariant matrices $A, B, C,$ and $D$ can continually retain preceding information. This is the “Long-Range Dependencies” shown in the center, and the paper that demonstrated the effectiveness of a structured matrix for this purpose was “HiPPO: Recurrent Memory with Optimal Polynomial Projections.”

An SSM fundamentally has a recurrent system. Its structure maps input to output, and the change in system state resulting from an earlier input and output affects subsequent inputs and outputs. This can be implemented inductively as recurrent operations, or simplified into a convolution. Yet the matrix operations that necessarily grow for high-dimensional data—for example, operations on $A, B, C, D,$ or $\bar{K}$—still create a problem of high computational cost.

What if we use a “well”-structured SSM? (S4)

This brings us to the objective of “Efficiently Modeling Long Sequences with Structured State Spaces” (S4). Methods had been proposed to reduce the high computational cost of conventional LSSL—a plain SSM—but all lacked numerical stability or rigor. To improve numerical stability while taking advantage of algorithms applicable to existing “well”-structured matrices, S4 restructures matrix $A$, the basis of the SSM. It exploits the fact that $A$ can be decomposed into a component whose rank is far lower than that of the original $A$ and a normal-matrix component. Rank is the number of independent columns; a low-rank matrix can be represented in fewer dimensions after dependent columns are removed. A normal matrix is a special square matrix whose properties—including commutation and diagonalizability—make it faster to compute with than a non-normal matrix.

Whereas the conventional SSM used a coefficient-space approach—the function representing the latent was expressed as a weighted sum of coefficients on predefined orthogonal functions—S4 lifts the computation into the frequency domain. Convolution in the time domain becomes multiplication in the frequency domain.

This lets the low-rank matrix be handled through the Woodbury identity and the normal matrix through a Cauchy kernel, reducing computation from $O(N^2L)$ to $\tilde{O}(N+L)$ and memory from $O(NL)$ to $O(N+L)$.

Conventional SSMs and their representation

We use matrix $A$ that is modeled to preserve long contexts, such as a HiPPO matrix. This constructs a system expressed by the following simultaneous differential equations:

$$\begin{aligned} \mathbf{x}^\prime = A\mathbf{x} + B\mathbf{u} \\ \mathbf{y} = C\mathbf{x}+D\mathbf{u} \end{aligned},\quad A^{\text{HiPPO}}_{n, k}= -\begin{cases} (2n+1)^{1/2}(2k+1)^{1/2},&\text{if }n > k \\ n+1,&\text{if }n = k \\ 0,&\text{if } n < k \end{cases}.$$

A conventional computing system must convert the continuous system into a discrete one with discretized inputs. Applying that conversion produces the discrete SSM used in practice:

$$\begin{aligned} x_k = \overline{A}x_{k-1} + \overline{B}u_k\quad &\overline{A} = (I-\Delta/2 \cdot A)^{-1} (I+\Delta/2\cdot A) \\ y_k = \overline{C}x_k,\quad &\overline{B} = (1-\Delta/2\cdot A)^{-1}\Delta B \\ &C = \overline{C} \end{aligned}$$

For the proof and more detail, see the previous post on LSSL. The important implication of the discrete SSM is that, because an SSM ultimately has a recurrent computational structure, it shares the computational characteristics of an RNN. The discretized matrix $\overline{A}$ serves as the transition matrix for hidden state $x$. If we expand the hidden states and outputs above while assuming $x_{-1} = 0$, however, the operation can also be expressed with a convolution kernel:

$$\begin{aligned} x_0 =& \overline{B}u_0,\quad x_1 = \overline{AB}u_0 + \overline{B}u_1,\quad x_2 = \overline{A}^2\overline{B}u_0+\overline{A}\overline{B}u_1 +\overline{B}u_2,~\cdots \\ y_0 =& \overline{CB}u_0,\quad y_1 = \overline{CAB}u_0 + \overline{CB}u_1,\quad y_2 = \overline{CA}^2\overline{B}u_0+\overline{CA}\overline{B}u_1 +\overline{CB}u_2,~\cdots \\ y_k =& \overline{K} \ast \mathbf{u},\quad \overline{K}\in\mathbb{R}^L := \mathcal{K}_L(\overline{A}, \overline{B}, \overline{C}) := (\overline{C}\overline{A}^i\overline{B})_{0\le i < L} \end{aligned}$$

If the convolutional filter $\overline{K}$ were known, an FFT—a fast convolution algorithm—could improve computation speed. But computing the filter itself requires matrix multiplication, and generalizing the operation to a non-normal matrix $\overline{A}$ remains difficult. The paper’s task is therefore to compute this filter efficiently.

Method: Diagonalization

Matrix diagonalization uses a matrix’s eigenvalues $\lambda$ and eigenvectors to transform it into a matrix whose diagonal entries are the eigenvalues. If a matrix $A \in \mathbb{R^{n \times n}}$ is diagonalizable, for example, its $n$ eigenvalues and corresponding $n$ eigenvectors let us write $\Lambda = V^{-1}AV$. Here $\Lambda$ and $V$ are

$$\text{For eigenvalues }\{ \lambda_i\}_{i=1}^n \text{ and eigenvectors } \{v_i\}_{i=1}^n,\quad\Lambda = \begin{bmatrix} \lambda_1 & 0 & 0 & \cdots & 0 \\ 0 & \lambda_2 & 0 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0& \lambda_{n-1} & 0 \\ 0 & 0 & 0 & 0 & \lambda_n \end{bmatrix},\quad V = \begin{bmatrix} \mid & \mid & \cdots & \mid \\ v_1 & v_2 & \cdots & v_n \\ \mid & \mid & \cdots & \mid \\ \end{bmatrix}$$

Under these conditions, we can rewrite the original equations as follows. SSMs generally simplify with $D = 0$, so we will write $\mathbf{y} = C\mathbf{x}$.

$$\begin{aligned} \mathbf{x}^\prime =& A\mathbf{x} + B\mathbf{u} \\ \mathbf{y} =& C\mathbf{x} \end{aligned}\quad\rightarrow\quad \begin{aligned} \tilde{\mathbf{x}}^\prime =& V^{-1}AV\tilde{\mathbf{x}} + V^{-1}B\mathbf{u} \\ \mathbf{y} =& CV\tilde{\mathbf{x}} \end{aligned}$$

The two systems may look different, but multiplying every expression on the left side of the transformed system by $V$ shows that it is equivalent to an SSM with $x = V\tilde{x}$. The input-to-output relationship $\mathbf{u} \rightarrow \mathbf{y}$ remains the same; only the system latent changes according to eigenvector matrix $V$ of matrix $A$. If we can construct the diagonal matrix $V^{-1}AV$ that plays the role of $A$ in the transformed system, we can efficiently reduce the cost of repeated multiplication by $A$ in the convolution above. This is called a Vandermonde product; we will simply take that as given.

Unfortunately, the HiPPO matrix cannot be diagonalized stably.

The reason is that when the HiPPO matrix is diagonalized, the entries of the eigenvector matrix $V$ become far too large. Put simply, numerically stable computation—agreement with the true mathematical value—requires matrix entries not to grow excessively during matrix operations. Diagonalizing the HiPPO matrix $A$ shown earlier produces entries $V_{ij} = {i+j \choose i-j}$. As the state size $N$ grows, entries as large as $2^{4N/3}$ become unmanageable and can, for example, break the output computation $CV\tilde{\mathbf{x}}$.

Method: Normal + Low-Rank Matrix

The preceding section attempted to make computation easier through diagonalization, but showed that applying it to a basic HiPPO matrix requires additional work. The ideal case is for the target matrix $A$ to be diagonalizable by a special matrix such as a unitary matrix. In linear algebra, matrices $A$ satisfying this condition are called normal matrices. As you may have guessed, a HiPPO matrix is not normal, which is why its eigenvector entries grow during diagonalization.

Fortunately, the authors discovered that although HiPPO matrix $A$ is not normal, it can be represented as the sum of a normal matrix and a low-rank matrix. A new problem appears, however: computing the convolution filter requires taking powers of this sum (Normal + Low Rank), which is itself slow and needs optimization.

$$\overline{K}\in\mathbb{R}^L := \mathcal{K}_L(A, B, C) := (CA^iB)_{0\le i < L}$$

To address the problem, the kernel filter is computed with three additional algorithms. The following figure summarizes how they relate. We have not yet covered enough detail to understand it, so let us proceed in order.

To follow the algorithm, we first need to accept that matrix $A$ can be represented as NPLR (Normal Plus Low-Rank) or DPLR (Diagonal Plus Low-Rank). Appendix C.1 of the paper demonstrates this empirically for every HiPPO matrix. There is no real need to understand the equations in this part; it can simply be accepted as a fact.

The method then departs from the conventional calculation of the convolution kernel. Instead of directly computing $\overline{K}$, it uses $\hat{K}$, the Discrete Fourier Transform (DFT) of $\overline{K}$. The DFT converts a signal on a discrete time axis—here, think of the filter entries $CB, CAB, CA^2B, \cdots$ as samples on a continuous time axis—into a spectrum on a discrete frequency axis. The algorithms for the DFT and its inverse, the IDFT, are known as the Fast Fourier Transform (FFT), which runs in $O(L\log L)$.

$$\hat{K}_j = \sum_{k=0}^{L-1} \overline{K}_k\exp(-2\pi j\frac{k}{L})$$

We will cover the details below, but constructing a truncated SSM and computing in the spectral domain replaces repeated multiplications by $A$ during filter computation with a single matrix operation. This introduces the troublesome term $(1-\overline{A}^L)$—which apparently still requires exponentiation—but reparameterization avoids repeated computation, saves memory, and improves speed. The structure assumed above then shows that computing the spectral kernel $\hat{K}$ is equivalent to a Cauchy kernel, allowing an efficient algorithm. In brief, the sequence is as follows.

1. Every HiPPO matrix $A$ can be represented as NPLR (or DPLR). This structure reduces the discretization $A \rightarrow \overline{A}$ to $O(N)$ computation.
2. The truncated-SSM generating function of $\overline{K}$ is equivalent to the DFT. We can therefore compute through transforms into and out of the frequency domain, replacing repeated powers of $\overline{A}$ with a single operation.
3. This operation has the same structure as a Cauchy-kernel computation, so an efficient algorithm is available.
4. The inverse can be simplified with the Woodbury identity.

A more detailed explanation of the method

Let us examine the efficiency improvements in detail. Personally, I find this harder than diffusion papers. Still, let us press on.

Every HiPPO matrix can be represented as NPLR, or Normal Plus Low-Rank. Four HiPPO matrices exist; considering only the most common case, LegS, gives

$$A^{\text{HiPPO}} _{n, k}= -\begin{cases} (2n+1)^{1/2}(2k+1)^{1/2},&\text{if }n > k \\ n+1,&\text{if }n = k \\ 0,&\text{if } n < k \end{cases}.$$

Adding $\frac{1}{2}(2n+1)^{1/2}(2k+1)^{1/2}$ to every entry produces

$$-\begin{cases} \frac{1}{2}(2n+1)^{1/2}(2k+1)^{1/2},&\text{if }n > k \\ \frac{1}{2},&\text{if }n = k \\ -\frac{1}{2}(2n+1)^{1/2}(2k+1)^{1/2},&\text{if } n < k \end{cases}.$$

Separating the diagonal entries gives $-\frac{1}{2}I + S$, where $S$ is a skew-symmetric matrix, which is normal. The matrix that adds the same value to every entry of the original matrix also has rank $1$. See the paper for proofs covering every HiPPO matrix. The implication is that every HiPPO state matrix can be separated into a diagonal part and a low-rank part. We will write this as follows. The paper uses the conjugate symbol ($\ast$) in its proof; because that symbol is less familiar to me, I will replace it with the transpose symbol ($\top$).

$$A = \Lambda - PQ^\top$$

Using this substitution to express the discrete system matrices gives

$$\begin{aligned} \overline{A} &= (I-\Delta / 2 \cdot A)^{-1} (I+\Delta / 2 \cdot A) \\ \overline{B} &= (1-\Delta / 2 \cdot A)^{-1} \Delta B \\ \\ I + \Delta/2\cdot A &= I + \Delta / 2 \cdot\left(\Lambda - PQ^\top\right) \\ &= \frac{\Delta}{2} \left(\frac{2}{\Delta}I + \Lambda - PQ^\top\right) \\ &= \frac{\Delta}{2}A_0 \\ \\ (I-\Delta/2 \cdot A)^{-1} &= \frac{2}{\Delta}\left(\frac{2}{\Delta}-\Lambda+PQ^\top\right)^{-1} \\ &= \frac{2}{\Delta}\left(D-DP(I+Q^\top DP)^{-1}Q^\top D \right) \\ &= \frac{2}{\Delta}A_1, \text{ where } D = \left(\frac{2}{\Delta}-\Lambda\right)^{-1} \text{ (Diagonal term)} \end{aligned}$$

Computing $A_0$ poses no serious problem, but $A_1$ contains a major obstacle: the matrix inverse. The computation grows dramatically as the matrix dimension $N$ increases. The inverse of a DPLR matrix can be simplified as above with the Woodbury identity. A diagonal matrix is easy to invert, and the attached low-rank term can be handled without inverting the entire expression. Because the Woodbury identity appears repeatedly in the derivation below, it is useful to keep it in mind: it can be applied efficiently whenever a matrix has DPLR structure.

Woodbury's identity applies as follows. For three matrices whose elements form a commutative ring, $A\in\mathcal{R}^{N \times N}, U, V \in\mathcal{R}^{N \times p}$—where $\mathcal{R}$ is a commutative ring whose elements commute under multiplication—

$$(A+UV^\top)^{-1} = A^{-1}-A^{-1}U(I_p + V^\top A^{-1}U)^{-1}V^\top A.$$

Using the resulting expressions, the discrete system can be defined for the DPLR matrices $A_1$ and $A_0$ with $O(N)$ computation:

$$\begin{aligned} x_k =& A_1A_0x_{k-1} + 2A_{1}Bu_k \\ y_k =& C^\top x_k \end{aligned}.$$

At this point the authors treat the former row vector $C$ as a column vector so that its shape matches the other parameters $B, P,$ and $Q$. Following them, from here onward I will write the $C$ in the original system equation as $C^\top$.

We now understand the discretized system matrices. The most important task is to compute the convolution filter $\overline{K}$ of the recurrent system quickly. DPLR improves the efficiency of matrix discretization by enabling the Woodbury identity, but it does not help much with the repeated multiplications needed to compute the convolution filter itself. Exploiting DPLR requires an inverse, so we can instead move to the spectral domain with a generating function of the coefficients. Suppose, for example, that we have an infinite convolution-filter signal:

$$\mathcal{K}(\overline{A}, \overline{B}, \overline{C}) = \{\overline{C}^\top\overline{B},\overline{C}^\top \overline{A}\overline{B},\overline{C}^\top \overline{A}^2\overline{B},\ldots\}$$

What we actually have is a finite filter of length $L$, truncated from that ideal convolution filter:

$$\mathcal{K}_L(\overline{A}, \overline{B}, \overline{C}) = \{\overline{C}^\top\overline{B},\overline{C}^\top \overline{A}\overline{B},\overline{C}^\top \overline{A}^2\overline{B},\ldots,\overline{C}^\top \overline{A}^{L-1}\overline{B}\}$$

A discrete signal of length $L$ can be decomposed into components with frequencies from $2\pi\times\frac{0}{L}$ through $2\pi \times \frac{L-1}{L}$. If a variable $z$ represents a frequency unit, the signal can be replaced by a set of coefficients of a function of $z$; this is called the $z$-transform. In general, $z$ is a complex number (Real + Imaginary). In frequency units, Euler’s angular conversion $e^{-i\Omega}$ lets us express it over $\Omega =\{\frac{2\pi l}{L}\}_{l=0}^{L-1}$.

$$\hat{K}_L(z; \overline{A}, \overline{B}, \overline{C}) := \sum_{i=0}^{L-1} \overline{C}^\top\overline{A}^i\overline{B}z^i = \overline{C}^\top(I-\overline{A}^L)(I-\overline{A}z)^{-1}\overline{B}$$

This transform is the DFT (Discrete Fourier Transform), commonly used to convert a signal on a discrete time axis into one on a discrete frequency axis. Setting the leading factor to $\overline{C}^\top(I-\overline{A}^L) = \tilde{C}$ gives

$$\hat{K}_L(z; \overline{A}, \overline{B}, \overline{C}) = \tilde{C} (1-\overline{A}z)^{-1}\overline{B}$$

We can substitute the closed forms of discretized matrices $\overline{A}, \overline{B}$ to express this in terms of $A, B$:

$$\begin{aligned} \overline{A} =& (I-\Delta/2 \cdot A)^{-1} (I+\Delta/2\cdot A) \\ \overline{B} =& (1-\Delta/2\cdot A)^{-1}\Delta B \end{aligned}$$

Substituting directly into the expression above yields

$$\begin{aligned} \tilde{C}^(I-\overline{A}z)^{-1}\overline{B} =& \tilde{C}\left(I-(I-\Delta/2 \cdot A)^{-1} (I+\Delta/2\cdot A)z \right)^{-1}\overline{B} \\ =& \tilde{C}\overline{B}\left(I-\frac{\Delta}{2}A\right)\left(\left(I-\frac{\Delta}{2}A\right)-\left(I+\frac{\Delta}{2}A\right)z\right)^{-1} \\ =& \tilde{C} \Delta B\left( I(1-z) - \frac{\Delta}{2}A(1+z)\right) \\ =& \frac{\Delta}{1-z}\tilde{C} \left( I - \frac{\Delta A}{2\frac{1-z}{1+z}} \right)^{-1}B \\ =& \frac{2}{1+z}\tilde{C}\left(\frac{2}{\Delta}\frac{1-z}{1+z}I -A \right)^{-1}B \end{aligned}$$

Because we showed that the system matrix is DPLR (Diagonal Plus Low-Rank), we can rewrite this as

$$\frac{2}{1+z}\tilde{C}\left(\frac{2}{\Delta}\frac{1-z}{1+z}I -A \right)^{-1}B = \frac{2}{1+z}\tilde{C}\left(\frac{2}{\Delta}\frac{1-z}{1+z}I -\Lambda+PQ^\top \right)^{-1}B$$

Now we can use the Woodbury identity introduced earlier. Defining the following quantity to simplify the expression,

$$R(z) = \left( \frac{2}{\Delta}\frac{1-z}{1+z} - \Lambda \right)^{-1}$$

leads to the final expansion

$$\tilde{C}(1-\overline{A}z)^{-1}\overline{B} = \frac{2}{1+z}\left(\tilde{C}R(z)B - \tilde{C}R(z)P(1+Q^\top R(z)P)^{-1} Q^\top R(z)B \right).$$

One question concerns the earlier definition $\overline{C}^\top(I-\overline{A}^L) = \tilde{C}$. Ordinarily, $\overline{A}^L$ would have to be computed repeatedly. By reparameterizing the trainable parameter so that it is initialized directly as $\tilde{C}$, however, we can reduce the associated computational cost.

We have reached the final step. To summarize, the operation that computes $\overline{K}$ has been replaced by the generating function $\hat{K}_L(\Omega; \overline{A}, \overline{B}, \overline{C})$, and the result was expanded by exploiting the diagonal component of $\overline{A}$. The resulting expression matches a Cauchy kernel exactly, and efficient algorithms exist for computing Cauchy kernels. A Cauchy matrix or kernel has the structure

$$M \in \mathbb{C}^{M \times N} = M(\Omega, \Lambda) = (M_{ij})_{0<=i < M,~0<=j<N},\quad M_{ij} = \frac{1}{\omega_i - \lambda_j}$$

Looking at the $Q^\top R(\Omega, \Lambda)P$ part of the final expression—where the original $z$ is considered at each element of unit $\Omega$—we see that it can be computed with the cost of a Cauchy matrix-vector multiplication. Ordinarily, computing a length-$L$ convolution kernel for $N$ hidden states costs $O(LN)$. Allowing a small approximation error reduces this to $O((L+N) \log(L+N) \log \frac{1}{\epsilon})$. The proof is mostly unnecessary: computing $Q^\top R(\Omega, \Lambda)P$ means obtaining $\sum_{j} \frac{q_j^\top p_j}{\omega-\lambda_j}$ for every $\omega$ in $\Omega$, which is precisely the form of a Cauchy kernel. The condition $\omega \neq \lambda_j$ always holds.

That completes the algorithm. Let us revisit the brief outline from earlier.

Beginning from the fact that $A$ is DPLR, we transformed into the Fourier domain to compute the generating function $\hat{K}$ rather than $\overline{K}$ directly. We reparameterized $\overline{C}$ in the resulting expression. We then wrote every discretized matrix in the closed form of the state matrix and expanded it with the Woodbury identity, obtaining

$$\begin{bmatrix}\tilde{C}^\top & Q\end{bmatrix}^\top \left(\frac{2}{\Delta} \frac{1-\omega}{1+\omega} -\Lambda \right)^{-1} \begin{bmatrix}B&P\end{bmatrix}$$

One slight difference between my derivation and the paper’s actual algorithm is that I have treated $\tilde{C}^\top = \tilde{C}$ throughout. Because the structure is the same, I do not believe this creates a serious problem. Applying an inverse Fourier transform to the efficiently computed $\hat{K}$ gives $\overline{K}$. At last, we have understood every algorithm behind the paper’s core idea.

Deep-layer architecture built with S4

The preceding derivation shows which parameters S4 needs. $A$ is initialized as some form of HiPPO matrix, and according to that form it consists of the diagonal and vector components $\Lambda, P, Q, B,$ and $C$. Since a diagonal matrix need not store parameters outside its diagonal entries, S4 has a total of $5N$ trainable parameters for state dimension $N$. S4 itself is a linear mapping—it emits a sequence of the same length—but stacking several S4 layers and adding nonlinearities turns it into a deep layer.

In conclusion

I think it is best to end here rather than separately present the experimental results. I expected them to demonstrate efficiency while retaining long-range effectiveness, and that is exactly what the paper shows. What I found most important was seeing how one can continually identify improvements over the preceding paper that established the modeling foundation for SSM—for example, decomposing a structure to simplify matrix multiplication and proving the result mathematically—and produce a sequence of contributions through that process.