Linear ODEs as state-space systems

1.1 What is a state-space system?

A linear state-space system is the pair of equations

where $\statevec(t) \in \R^N$ is the state at time $t \ge 0$, $u(t) \in \R^D$ is the input signal, $y(t) \in \R^M$ is the output, and the four matrices $(\statemat, \inputmat, \outputmat, \feedmat)$ are real-valued with shapes $N \times N$, $N \times D$, $M \times N$, $M \times D$. The first equation is the state equation — a first-order linear ODE describing how the state evolves. The second is the output equation — a static linear readout. The integer $N$ is the state dimension; it controls how much history the model can carry forward.

The pair is “linear” because both equations are linear in $\statevec$ and $u$ separately. It is “state-space” because $\statevec(t)$ summarizes everything the system needs from the past in order to predict the future given future inputs — a property called the Markov property of state. Given the present state and the input from now onward, the system’s entire future is determined.

A few conventions before going further. We treat $\statevec(t)$ as a column vector and apply $\statemat$ on the left: $\statemat \statevec$, not $\statevec \statemat$. We absorb the feedthrough term $\feedmat u(t)$ into the output equation only — it appears nowhere in the state dynamics — and we will often set $\feedmat = 0$ in examples to keep notation light. When the input is absent ($u \equiv 0$), the state equation reduces to

the homogeneous linear ODE, whose solution is the focus of §1.2.

Why this object

Three reasons the state-space form is the right starting point for a book on sequence-model architectures:

1. It separates dynamics from input. Whatever the input signal looks like, the dynamics matrix $\statemat$ alone determines the system’s long-run behavior — its stability, oscillation frequencies, and decay rates. Most of the analysis in Chapter 2 looks only at $\statemat$.
2. It generalizes immediately. Replacing the scalar derivative $\ddt$ with a finite difference yields a discrete recurrence (Chapter 4); replacing $\R$ with $\C$ admits oscillatory modes (Chapter 6, Chapter 10); replacing the constant $\statemat$ with an input-dependent $\statemat(u_t)$ yields the selective-SSM family (Chapter 9). Each generalization preserves the state-space skeleton.
3. It’s the form in which continuous-time models are discretized. Every SSM derivation in the literature starts from a continuous state equation and applies a numerical integration scheme to get a discrete recurrence. If you don’t know the continuous object, the discrete one is opaque.

1.2 The matrix exponential and fundamental solutions

The homogeneous ODE $\ddt \statevec = \statemat \statevec$ with initial condition $\statevec(0) = \statevec_0$ has a unique solution for every $\statevec_0 \in \R^N$ — and that solution is

where $e^{\statemat t}$ is the matrix exponential of $\statemat t$, defined by the same series as the scalar exponential:

The fundamental solution property — that $t \mapsto e^{\statemat t} \statevec_0$ solves the homogeneous ODE — is verified by differentiating the series term-by-term and using $\frac{d}{dt} \frac{(\statemat t)^k}{k!} = \frac{\statemat^k t^{k-1}}{(k-1)!}$. The resulting derivative is $\statemat \cdot e^{\statemat t} \statevec_0$, which matches $\statemat \statevec(t)$. Existence and uniqueness for the inhomogeneous case ($u \not\equiv 0$) follow the variation of parameters formula

an integral that recurs (in discretized form) as the convolution kernel of every LTI SSM in Chapter 8.

Eigenvalue structure determines everything

The matrix exponential’s behavior is governed entirely by the eigenvalues of $\statemat$. If $\statemat$ is diagonalizable, $\statemat = V \Lambda V^{-1}$ with $\Lambda = \diag(\lambda_1, \ldots, \lambda_N)$, then

So $e^{\statemat t}$ acts on the eigenbasis as independent scalar exponentials. The real parts of $\lambda_i$ control decay or growth: $\Re(\lambda_i) < 0$ means $e^{\lambda_i t} \to 0$; $\Re(\lambda_i) > 0$ means blow-up. The imaginary parts control oscillation frequency. Three regimes recur throughout the book:

1. All eigenvalues with negative real part — the system is asymptotically stable: every trajectory $\statevec(t) \to 0$ as $t \to \infty$. This is the regime in which a recurrence “forgets” its initial state, and it is the design target for the $\statemat$ matrices of S4, Mamba, and friends.
2. At least one eigenvalue with positive real part — the system is unstable: some initial conditions blow up. A trained SSM whose $\statemat$ drifts into this regime is the dominant numerical failure mode (Chapter 5).
3. Purely imaginary eigenvalues — the system is marginally stable and oscillates without decay. These are the energy-preserving modes of Hamiltonian systems and are the natural home of symplectic discretization (Chapter 6).

When $\statemat$ is not diagonalizable — when it has repeated eigenvalues with deficient eigenspaces — the matrix exponential picks up polynomial-in-$t$ factors on top of the $e^{\lambda_i t}$ terms. The Jordan normal form (Chapter 3, §3.1) gives the precise statement. In practice almost every $\statemat$ arising from random initialization or from physical systems is diagonalizable, and the polynomial-factor case shows up only at codimension-one boundaries in parameter space.

1.3 The damped harmonic oscillator

The canonical small example. A unit mass on a spring with stiffness $k > 0$ and damping coefficient $c > 0$, displacement $q(t)$, evolves according to Newton’s second law,

This is a second-order scalar ODE, not first-order, so it does not yet fit the state-space template. We lift it by treating velocity as an extra state coordinate: let $\statevec(t) = (q(t), \dot q(t))^\top \in \R^2$. Then $\ddt \statevec = \statemat \statevec$ with

This lifting trick — introduce auxiliary states to reduce ODE order to one — is the same procedure by which every higher-order ODE becomes a first-order state-space system, and it foreshadows how implicit higher-order schemes (Chapter 6) work internally.

The eigenvalues of $\statemat$ are the roots of the characteristic polynomial $\lambda^2 + c \lambda + k = 0$:

Three sub-cases follow from the sign of the discriminant $c^2 - 4k$:

• Overdamped ($c^2 > 4k$): two real negative eigenvalues. Both modes decay monotonically; the displacement approaches zero without oscillation.
• Critically damped ($c^2 = 4k$): a repeated negative eigenvalue $\lambda = -c/2$. The matrix $\statemat$ is in this case not diagonalizable, and the solution contains a $t \cdot e^{\lambda t}$ term.
• Underdamped ($c^2 < 4k$): a complex-conjugate eigenpair $\lambda_\pm = -c/2 \pm i\omega$, with $\omega = \tfrac{1}{2}\sqrt{4k - c^2}$. Trajectories spiral toward the origin: oscillation at angular frequency $\omega$ with envelope decaying as $e^{-ct/2}$.

The total energy $E(t) = \tfrac{1}{2}(\dot q^2 + k q^2)$ is a useful diagnostic. Its time derivative is $\dot E = -c \dot q^2 \le 0$, i.e. energy decreases monotonically whenever $\dot q \ne 0$. Both stationary and oscillating trajectories satisfy this. The companion damped_oscillator.py simulates the underdamped regime and plots $E(t)$ alongside the trajectory; the energy curve gives an at-a-glance check that the numerical integration preserves the dissipation structure of the continuous system.

1.4 Coupled systems and the Jacobian

The damped oscillator is a 2-dimensional toy. Real systems — and trained SSMs — live in much higher dimensions, and they typically arise from coupling many simpler subsystems. The bookkeeping is straightforward once you accept that the state vector $\statevec$ stacks all subsystem states and $\statemat$ encodes both intra-subsystem dynamics and inter-subsystem coupling.

A useful larger example is a ring of coupled oscillators: $n$ identical damped oscillators arranged on a circle, each coupled to its two nearest neighbors by springs of stiffness $\kappa$. The displacements $q_1, \ldots, q_n$ and velocities $\dot q_1, \ldots, \dot q_n$ stack into a state vector $\statevec \in \R^{2n}$. The $\statemat$ matrix has a $2 \times 2$ block-diagonal structure for the per-oscillator dynamics, plus an off-diagonal coupling pattern that links each oscillator’s velocity equation to its neighbors’ positions:

with indices taken mod $n$. The bracketed term $q_{i-1} - 2 q_i + q_{i+1}$ is the discrete Laplacian of $q$ around the ring; it’s the same discretization of $\partial^2/\partial x^2$ you would use in a finite-difference method for the wave equation.

The $2n \times 2n$ Jacobian matrix $\jacobian$ — which for a linear system is just $\statemat$ — has eigenvalues that are easy to compute analytically because of the ring’s circulant structure. They come in complex-conjugate pairs

(for the underdamped regime). Each $m$ indexes a standing wave mode on the ring with wavenumber $2\pi m / n$; the $m=0$ mode is uniform translation (no spatial structure), $m = n/2$ (when $n$ is even) is the maximally-oscillating mode where adjacent oscillators move in antiphase. The companion coupled_oscillators.py constructs $\statemat$ for $n = 8$ and visualizes the eigenvalues in the complex plane.

The key qualitative lesson: once you have the right state-space lift, even a system with rich spatial structure reduces to “compute eigenvalues of $\statemat$, classify by real and imaginary parts.” The eigenvalue distribution of a trained SSM’s $\statemat$ — pre-training, mid-training, and converged — is one of the most informative diagnostics you can compute (Chapter 2, §2.2; Chapter 5).

The Jacobian beyond linear systems

For a nonlinear system $\ddt \statevec = f(\statevec)$, the Jacobian $\jacobian(\statevec) = \partial f / \partial \statevec$ is the matrix of partial derivatives at $\statevec$. The linearization $\ddt \delta\statevec = \jacobian(\statevec^*) \delta\statevec$ around a fixed point $\statevec^*$ (where $f(\statevec^*) = 0$) describes how small perturbations $\delta\statevec$ evolve. The eigenvalues of $\jacobian(\statevec^*)$ then determine local stability of $\statevec^*$ by the same classification as in §1.2.

Trained recurrent networks — including SSMs in inference mode — are formally nonlinear (the state update is state = f(state, input) where f includes the input-dependent matrices). The linearization-around-a-state gives a local state-space system at each time step, and tracking how that local $\statemat$ varies along a trajectory is the foundation of the Lyapunov-exponent analysis in Chapter 2.

1.5 Phase portraits and the eigenvalue classification

A phase portrait is the geometric picture of all trajectories of an autonomous system in state space, drawn as oriented curves. For 2-dimensional linear systems $\ddt \statevec = \statemat \statevec$ on $\R^2$, the topology of the portrait is determined entirely by the trace $\tr(\statemat)$ and the determinant $\det(\statemat)$ — equivalently, by the eigenvalues — and the standard classification gives six qualitative behaviors:

The damped oscillator’s three sub-cases (over/critical/under-damped) realize three of these: stable node, stable degenerate node, stable spiral. The undamped limit $c \to 0$ slides the eigenvalues onto the imaginary axis, turning the spiral into a center; this is the Hamiltonian limit where energy is conserved exactly.

For higher-dimensional systems the topology is richer (think: 3-D systems can have spirals around line-shaped invariant manifolds), but the building blocks are the same: each pair of complex-conjugate eigenvalues contributes a 2-D spiral subspace, each real eigenvalue contributes a 1-D direction of decay or growth. This decomposition by invariant subspaces — equivalently, by the Jordan blocks of $\statemat$ — is the geometric content of the matrix exponential’s eigenvalue structure (§1.2).

1.6 A preview of frequency response

The variation-of-parameters formula

expresses the state as a convolution of the input with the kernel $K(t) := e^{\statemat t} \inputmat$. The output $y(t) = \outputmat \statevec(t)$ (taking $\feedmat = 0$) is therefore

The function $h(t) := \outputmat e^{\statemat t} \inputmat$ for $t \ge 0$ (and zero otherwise) is the impulse response of the system: the output you would observe if the input were a Dirac delta $u(t) = \delta(t)$.

Taking the Laplace transform turns convolution into multiplication. With $\widehat y(s) := \int_0^\infty y(t) e^{-st} dt$ and similarly for $\widehat u$, you get

where $H(s) := \outputmat (s I - \statemat)^{-1} \inputmat$ is the transfer function. The eigenvalues of $\statemat$ — already shown to govern $e^{\statemat t}$ — appear as the poles of $H(s)$, since $(sI - \statemat)^{-1}$ blows up exactly when $s$ equals an eigenvalue of $\statemat$. The poles’ location in the complex $s$-plane mirrors the eigenvalue classification of §1.5: left-half-plane poles ⇔ stable system; right-half-plane poles ⇔ unstable; imaginary-axis poles ⇔ marginally stable / oscillatory.

This connection is the bridge to Chapter 8, where you’ll see that the S4 layer’s convolutional view is literally a numerical evaluation of the impulse response $h(t)$ at uniformly spaced time points, and to Chapter 11, where the spectral structure of long convolution kernels (Hyena, RetNet) is analyzed via the same transfer-function framing.

1.7 What’s next

The next two chapters develop the analytical machinery that this chapter has only sketched. Chapter 2 makes the eigenvalue-stability story rigorous (Lyapunov, A-stability, BIBO) and adds the QR-based computational method for tracking eigenvalue trajectories during training. Chapter 3 introduces the structured-matrix vocabulary (Toeplitz, semiseparable, Vandermonde, Cauchy) needed to discuss the kernel constructions of S4 and friends. Chapters 4–6 then take the continuous state-space system here and discretize it — first crudely (ZOH, bilinear), then with full numerical-analysis machinery (Butcher tableau, A-stability), then via implicit and structure-preserving methods (Gauss–Legendre, symplectic integrators), which is where the active research pilot’s anchor (the C1 pilot) lives.

You can also profitably jump ahead: Chapter 8’s LTI SSM presentation reuses everything in this chapter and is self-contained for readers who want to see the SSM application before working through the full math foundation.

1.8 Exercises

Seven problems mixing computation and theory. Short/numerical exercises (1–4) have inline collapsible solutions; long/proof exercises (5–7) have full worked solutions in §1.9.

Exercise 1.1 (computation)

Compute the matrix exponential of $\statemat = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix}$ by hand, term-by-term, and verify it equals the rotation matrix $R(t) = \begin{pmatrix} \cos t & \sin t \\ -\sin t & \cos t \end{pmatrix}$.

Exercise 1.2 (computation)

For the damped oscillator $\ddot q + c\dot q + k q = 0$ with $k = 1$ and $c = 0.5$, classify the damping regime and find the eigenvalues of $\statemat$.

Exercise 1.3 (computation)

Lift the third-order ODE $q''' + 2 q'' + q = u(t)$ into a 3-dimensional state-space system $\ddt \statevec = \statemat \statevec + \inputmat u$. Write $\statemat$ and $\inputmat$ explicitly.

Exercise 1.4 (computation)

For the ring of $n = 4$ identical damped oscillators with $k = 1$, $c = 0$, $\kappa = 0.5$ (no damping, symmetric coupling), use the formula in §1.4 to find the eigenvalues of $\statemat$ analytically. Confirm they lie on the imaginary axis (since $c = 0$) and verify by running companions/ch01/jax/coupled_oscillators.py with $n = 4$.

Exercise 1.5 (theory) — solution in §1.9

Prove that the matrix-exponential series $\sum_k M^k/k!$ converges absolutely (entrywise) for every square matrix $M$, with bound $\norm{e^M}_F \le e^{\norm{M}_F}$ where $\norm{\cdot}_F$ is the Frobenius norm.

Exercise 1.6 (theory) — solution in §1.9

Show that if $\statemat$ is skew-symmetric ($\statemat^\top = -\statemat$), then $e^{\statemat t}$ is orthogonal for every $t \in \R$ (that is, $(e^{\statemat t})^\top e^{\statemat t} = I$). Use this to give a geometric reason why the rotation matrix of Exercise 1.1 is orthogonal.

Exercise 1.7 (theory) — solution in §1.9

Prove the variation of parameters formula: the unique solution of $\ddt \statevec = \statemat \statevec + \inputmat u(t)$ with $\statevec(0) = \statevec_0$ is

1.9 Full solutions to theory exercises

Solution to Exercise 1.5

Let $\norm{\cdot}_F$ denote the Frobenius norm, which satisfies the sub-multiplicative property $\norm{AB}_F \le \norm{A}_F \norm{B}_F$. Then $\norm{M^k}_F \le \norm{M}_F^k$ by induction. Therefore

The partial sums are entrywise bounded by terms of a convergent scalar series, so each matrix entry is Cauchy in $K$ and converges. The limit is what we call $e^M$. The bound $\norm{e^M}_F \le e^{\norm{M}_F}$ follows by passing to the limit. ∎

Solution to Exercise 1.6

We have $\ddt e^{\statemat t} = \statemat e^{\statemat t}$ and $\ddt (e^{\statemat t})^\top = (e^{\statemat t})^\top \statemat^\top$. Let $Q(t) := (e^{\statemat t})^\top e^{\statemat t}$. Then

since $\statemat^\top + \statemat = 0$ by skew-symmetry. So $Q(t) = Q(0) = I$ for all $t$, i.e. $e^{\statemat t}$ is orthogonal.

Geometric interpretation: the rotation matrix $R(t)$ of Exercise 1.1 was generated by $\statemat = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix}$, which is skew-symmetric. So $R(t)$ is orthogonal — which we already knew, but this derivation pins down why: orthogonality follows structurally from skew-symmetry of the generator, not from any property of the cos/sin functions. ∎

Solution to Exercise 1.7

Define $\statevec(t)$ by the formula on the right-hand side. We verify (a) the initial condition and (b) the ODE.

Initial condition. At $t = 0$ the integral vanishes (zero-length interval) and $e^{\statemat \cdot 0} = I$, so $\statevec(0) = \statevec_0$. ✓

ODE. Differentiate using the product rule on the integral (Leibniz’s rule for differentiating under the integral sign):

Uniqueness follows from the classical Picard–Lindelöf theorem applied to the right-hand side $f(\statevec, t) = \statemat \statevec + \inputmat u(t)$, which is Lipschitz in $\statevec$ uniformly in $t$ with constant $\norm{\statemat}$ (see Hairer–Nørsett–Wanner Hairer et al. (1993) for the standard proof). ∎

1.10 Companion code

Three runnable companions in companions/ch01/jax/:

• damped_oscillator.py — simulates the underdamped oscillator, plots energy decay (Figure 1.1)
• coupled_oscillators.py — constructs the ring-Laplacian state matrix, plots eigenvalue spectrum (Figure 1.2)
• matrix_exponential.py — compares scipy.linalg.expm against truncated series sums; the truncated series catastrophically diverges for matrices with large spectral radius (auxiliary, used by Exercise 1.5)

All companions import their plotting style from companions/_shared/plot_utils.py. To run from the repo root:

PYTHONPATH=. python companions/ch01/jax/damped_oscillator.py
PYTHONPATH=. python companions/ch01/jax/coupled_oscillators.py
PYTHONPATH=. python companions/ch01/jax/matrix_exponential.py

Figures are written to public/figures/ch01/ and referenced from the <Figure> components above.

Stability theory: Lyapunov, A-stability, BIBO

2.1 Three notions of stability

For the rest of this chapter, the system under analysis is the LTI state-space model from Chapter 1:

with $\statemat \in \R^{N \times N}$. The three notions:

1. Lyapunov stability asks: with $u \equiv 0$, do trajectories $\statevec(t) = e^{\statemat t} \statevec_0$ stay bounded as $t \to \infty$? Asymptotically stable if $\statevec(t) \to 0$.
2. A-stability is a property of a numerical integrator, not of the system. It asks: does the integrator preserve Lyapunov stability for every stable LTI test problem, at every step size $\stepsize > 0$? If yes, the method is A-stable.
3. BIBO stability asks: for the input–output mapping $u \mapsto y$, does every bounded input ($\sup_t \abs{u(t)} < \infty$) produce a bounded output? The system is BIBO-stable if yes.

When $\statemat$ has all eigenvalues in the open left half-plane (LHP), all three coincide for the LTI case. But:

• Lyapunov-stable systems can be Lyapunov-stable but not asymptotically stable (centers: pure imaginary eigenvalues, oscillations don’t decay).
• A-stability is a integrator property; the explicit Euler method is not A-stable, while the implicit Euler method is. The same continuous system therefore may need a stable integrator to preserve its stability.
• BIBO depends on the transfer function $H(s) = \outputmat (sI - \statemat)^{-1} \inputmat$ and on the choice of $(\statemat, \inputmat, \outputmat)$ realization. Two systems with the same $\statemat$ can have different BIBO properties if their $\outputmat$ “hides” unstable modes.

The rest of the chapter develops each notion in turn.

2.2 Lyapunov stability: theory

We restrict to the homogeneous case $\ddt \statevec = \statemat \statevec$ (zero input). The eigenvalue criterion is essentially the content of Chapter 1, §1.2, formalized:

The proof reduces to the matrix-exponential formula. If $\statemat = V \Lambda V^{-1}$ is diagonalizable, then $\statevec(t) = V e^{\Lambda t} V^{-1} \statevec_0$ and the $i$-th eigencomponent decays like $e^{\Re(\lambda_i) t}$. The Jordan-block caveat in case (1) handles the case where a defective imaginary eigenvalue produces polynomial growth in $t$ even though the real part is zero.

For nonlinear systems, the local-linearization picture from Chapter 1, §1.4 makes the same theorem applicable at each fixed point: the fixed point is asymptotically stable if the Jacobian’s eigenvalues all sit in the LHP. This is the principle that makes “check the eigenvalues” the first move in nearly every SSM stability analysis.

The Lyapunov equation

There’s an alternative characterization of asymptotic stability that doesn’t require computing eigenvalues. The continuous-time Lyapunov equation

(with $Q$ a given symmetric positive-definite matrix) has a unique symmetric positive-definite solution $P$ if and only if $\statemat$ is asymptotically stable. The function $V(\statevec) = \statevec^\top P \statevec$ is then a Lyapunov function for the system: $\dot V = -\statevec^\top Q \statevec < 0$ for all nonzero $\statevec$, so $V$ decreases along trajectories and forces them to the origin.

For SSM analysis, the Lyapunov-equation view is occasionally useful — for example, the controllability and observability Gramians of an LTI realization are Lyapunov solutions Antoulas (2005) — but the eigenvalue criterion is the more direct tool.

2.3 Lyapunov exponents: computation

The eigenvalues of $\statemat$ work as a stability diagnostic for LTI systems, where $\statemat$ is constant. For systems where the linearization $\jacobian(t)$ varies with time (e.g. a trained recurrent network in inference mode, where the per-step Jacobian depends on the input), eigenvalues at any single $t$ don’t tell the long-run story. The right tool is the Lyapunov exponent, which generalizes the eigenvalue’s real part to time-varying systems.

Three things make Lyapunov exponents the right object:

1. They generalize eigenvalues to time-varying systems. For a constant Jacobian $\jacobian_t \equiv \statemat$, $\lyapexp_i = \log\abs{\mu_i(\statemat)}$ where $\mu_i$ are the eigenvalues’ magnitudes (in the discrete case) or $\Re(\lambda_i)$ (in the continuous case). The connection back to Chapter 1’s eigenvalue classification is exact in the autonomous case.
2. The maximal Lyapunov exponent $\lyapexp_1$ measures sensitive dependence. A positive $\lyapexp_1$ means infinitesimally close initial conditions diverge exponentially (chaos); $\lyapexp_1 = 0$ is the edge-of-chaos regime; $\lyapexp_1 < 0$ means contracting dynamics.
3. They are robust to coordinate changes. Eigenvalues of the per-step Jacobian depend on the basis you write the state in; Lyapunov exponents do not.

The QR algorithm

Computing $\sigma_i(\Phi_T)$ for large $T$ directly is numerically impossible: even for a well-behaved system with $\lyapexp_1 = 0.1$, after $T = 100$ steps the largest singular value is $e^{10} \approx 22000$ and the smallest may be $e^{-10} \approx 5 \times 10^{-5}$. The product matrix $\Phi_T$ becomes either numerical zero or numerical infinity in finite-precision arithmetic, and finite-precision singular values lose all meaning. The fix, due to Benettin et al. Benettin et al. (1980) , is to factor the growth out at every step using QR decomposition.

The QR-based Lyapunov algorithm (the version that has become standard in dynamical-systems software):

Initialize Q_0 = I (identity), running sums L_i = 0 for i = 1..N.
For t = 1, 2, ..., T:
    M_t = J_t @ Q_{t-1}             # one-step forward propagation
    Q_t, R_t = qr(M_t)              # extract orthogonal + upper triangular
    Update L_i += log |R_t[i,i]| for each i
At the end:
    Lyapunov_i = L_i / T

The key insight: by re-orthonormalizing the propagated frame at every step, we factor the exponential growth into the diagonal of the upper-triangular $R_t$ rather than letting it accumulate in $\Phi_T$. The diagonal entries of $R_t$ are well-conditioned numbers near 1, and their logarithms sum nicely to give the Lyapunov exponents.

The companion lyapunov_qr.py implements this algorithm and applies it to the ring of coupled oscillators from Chapter 1, §1.4. The resulting spectrum’s structure — all exponents negative for the damped case, all near zero for the undamped case — matches the eigenvalue analysis exactly, validating the implementation.

2.4 A-stability: integrator regions of absolute stability

Switch perspectives. We now have a continuous system $\ddt \statevec = \statemat \statevec$ and we discretize it with some numerical integrator at step size $\stepsize$. The integrator approximates the true continuous trajectory by a discrete recurrence

where $R(z)$ is the integrator’s stability function, evaluated at the matrix argument $\stepsize \statemat$ via the matrix-functional-calculus convention. Different integrators give different $R$:

• Forward (explicit) Euler: $\statevec_{n+1} = \statevec_n + \stepsize \statemat \statevec_n$, so $R_{\text{FE}}(z) = 1 + z$.
• Backward (implicit) Euler: $\statevec_{n+1} = \statevec_n + \stepsize \statemat \statevec_{n+1}$, giving $R_{\text{BE}}(z) = 1/(1-z)$.
• Bilinear (Tustin / trapezoidal): $R_{\text{bil}}(z) = (1 + z/2)/(1 - z/2)$.
• Exact (zero-order hold): $R_{\text{ZOH}}(z) = e^z$. This is the integrator S4 uses by default.

The discrete recurrence is Lyapunov-stable if and only if every eigenvalue $\mu_i$ of $R(\stepsize \statemat)$ satisfies $\abs{\mu_i} \le 1$ (the discrete-time analog of the LHP criterion is the unit disk). Since the eigenvalues of $R(\stepsize \statemat)$ are $R(\stepsize \lambda_i)$ where $\lambda_i$ are the eigenvalues of $\statemat$, the question reduces to: for which $z = \stepsize \lambda_i$ in the complex plane is $\abs{R(z)} \le 1$?

A-stability matters because it means: if the continuous system is asymptotically stable, the integrator preserves that stability at every step size $\stepsize > 0$. Without A-stability, you have to pick $\stepsize$ small enough that $\stepsize \lambda_i$ lands inside the (bounded) stability region.

The four examples above:

• Forward Euler: $\abs{1 + z} \le 1$ ⇔ $z$ in the closed disk of radius 1 around $-1$. Tiny region; not A-stable. For a continuous system with eigenvalue $\lambda = -100$ (fast decay), forward Euler requires $\stepsize < 2/100$ to avoid blowup.
• Backward Euler: $\abs{1/(1-z)} \le 1$ ⇔ $\abs{1-z} \ge 1$ ⇔ $z$ outside the open disk of radius 1 around $1$. This contains the entire closed LHP. A-stable.
• Bilinear (trapezoidal): $\abs{R_{\text{bil}}(z)} \le 1$ ⇔ $\Re(z) \le 0$ exactly. The stability region is exactly the closed LHP. A-stable, and the boundary coincides with the imaginary axis — a property called L-stability when combined with $R(z) \to 0$ as $\Re(z) \to -\infty$, which bilinear does not satisfy ($R_{\text{bil}}(-\infty) = -1$).
• Zero-order hold (matrix exponential): $\abs{e^z} = e^{\Re(z)} \le 1$ ⇔ $\Re(z) \le 0$. The stability region is the closed LHP. A-stable and L-stable. This is the gold standard for LTI discretization — and it’s what S4 uses internally.

Why this matters for SSMs

The S4 family uses ZOH discretization, and is therefore stable at every step size for which the continuous system is stable. Mamba-1 uses ZOH as well. Mamba-3’s exponential-trapezoidal rule (Chapter 10) is a second-order generalization of bilinear that retains A-stability. The C1 research pilot (see Chapter 6) asks whether symplectic integrators — A-stable methods designed to preserve geometric invariants — can do better for complex-state SSMs whose $\statemat$ has near-imaginary eigenvalues.

When a stability-region analysis fails, it usually fails dramatically: a forward-Euler discretization of a stiff system blows up in a handful of steps. This is exactly the failure mode that careful integrator choice prevents.

For full coverage of A-stability theory and a wider zoo of integrator stability functions, Hairer–Wanner’s stiff-systems volume Hairer & Wanner (1996) is the canonical reference.

2.5 BIBO stability

The third notion looks at the input–output mapping rather than the autonomous dynamics. For the system $\ddt \statevec = \statemat \statevec + \inputmat u$, $y = \outputmat \statevec$, the output under a unit-impulse input is the impulse response

The system is BIBO-stable if its impulse response satisfies $\int_0^\infty \norm{h(t)} \, dt < \infty$, equivalently if $h \in L^1$. Under this condition, any bounded input $u$ produces a bounded output $y$ — Young’s inequality on convolutions gives the precise bound $\sup_t \norm{y(t)} \le \norm{h}_{L^1} \cdot \sup_t \abs{u(t)}$.

Poles of $H(s)$ are the values of $s$ where $\det(sI - \statemat) = 0$ — i.e. the eigenvalues of $\statemat$ — unless a pole–zero cancellation occurs because $\outputmat$ doesn’t “see” some eigenmode. This cancellation matters: a system can be internally unstable (some $\lambda_i$ in the right half-plane) but BIBO-stable if those unstable modes are unobservable from $\outputmat$ and unreachable from $\inputmat$. The Kalman canonical-decomposition theorem makes this precise; for our purposes, the eigenvalues of $\statemat$ in the LHP are sufficient for both internal and BIBO stability.

For an SSM viewed as a sequence-to-sequence map, BIBO is the natural input–output stability concept: bounded inputs (typical of token embeddings) should produce bounded outputs (typical of pre-softmax logits). When training an SSM produces an $\statemat$ whose eigenvalues drift into the right half-plane, both internal (Lyapunov) and input–output (BIBO) stability fail simultaneously, and the model’s outputs blow up. Tracking the eigenvalues of $\statemat$ during training is therefore a single diagnostic for both failure modes.

2.6 Connection to Lyapunov functions

The Lyapunov-function method gives a non-spectral way to prove stability. For nonlinear systems where eigenvalues don’t directly apply, finding a function $V: \R^N \to \R$ with $V(0) = 0$, $V(\statevec) > 0$ for $\statevec \ne 0$, and $\dot V(\statevec) := \nabla V(\statevec) \cdot f(\statevec) < 0$ for all $\statevec \ne 0$ proves asymptotic stability of $\statevec^* = 0$.

For our LTI case $f(\statevec) = \statemat \statevec$, the Lyapunov equation $\statemat^\top P + P \statemat = -Q$ (with $Q$ p.d.) gives $V(\statevec) = \statevec^\top P \statevec$ as a quadratic Lyapunov function — and its existence with $P$ p.d. is equivalent to asymptotic stability (as noted in §2.2). The interest in the function-method picture for our context is mainly conceptual: it explains why “energy” arguments work to prove stability of mechanical systems (the energy is the Lyapunov function), and it generalizes immediately to nonlinear systems where the eigenvalue criterion only gives local information at fixed points.

This connects forward to Chapter 6, where the question of which discretization preserves a Lyapunov function (or, more strongly, an energy invariant) leads naturally into symplectic methods.

2.7 What’s next

Chapter 3 develops the structured-linear-algebra vocabulary (SVD, Jordan form, condition number, Toeplitz/Vandermonde/Cauchy/semiseparable structure) needed for the SSM kernel constructions in Chapters 7–9. Chapter 4 takes the discretization story started in §2.4 and develops it systematically: order conditions, accuracy classes, the Butcher tableau, the bilinear and ZOH derivations in detail. Chapter 6 picks up the A-stability theme and pushes into implicit and structure-preserving integrators — symplectic schemes, geometric integration, and the C1 pilot’s home territory.

2.8 Exercises

Six problems. Inline-collapsible solutions for the shorter ones; full solutions for the longer theory problems in §2.9.

Exercise 2.1 (computation)

Use the eigenvalue criterion (Theorem 2.1) to classify the stability of $\ddt \statevec = \statemat \statevec$ for $\statemat = \begin{pmatrix} -1 & 4 \\ -1 & -1 \end{pmatrix}$.

Exercise 2.2 (computation)

For the matrix $\statemat = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix}$ of Exercise 1.1, find the eigenvalues and classify the Lyapunov stability. Is the system asymptotically stable?

Exercise 2.3 (computation)

The forward Euler method has stability function $R(z) = 1 + z$. For the continuous test problem $\dot x = \lambda x$ with $\lambda = -10$ (a stiff system with fast time scale), find the maximum step size $\stepsize$ for which forward Euler is stable.

Exercise 2.4 (computation)

Show that the bilinear method $R_{\text{bil}}(z) = (1 + z/2)/(1 - z/2)$ maps the closed left half-plane $\{\Re(z) \le 0\}$ onto the closed unit disk $\{\abs{w} \le 1\}$. (This is the Möbius map picture of the bilinear transform.)

Exercise 2.5 (theory) — solution in §2.9

State and prove the Lyapunov equation characterization of asymptotic stability: $\statemat$ is asymptotically stable iff for every symmetric positive-definite $Q$, the equation $\statemat^\top P + P \statemat = -Q$ has a unique symmetric positive-definite solution $P$.

Exercise 2.6 (theory) — solution in §2.9

Prove the BIBO criterion of Theorem 2.3: a causal LTI system is BIBO-stable iff all poles of its transfer function lie in the open left half-plane. (Hint: Use the Laplace-transform / convolution duality, and exploit the representation of the impulse response as a sum of exponentials.)

2.9 Full solutions to theory exercises

Solution to Exercise 2.5

Forward direction (Lyapunov equation ⇒ asymptotic stability): Suppose $P$ is symmetric positive-definite and satisfies $\statemat^\top P + P \statemat = -Q$ for some symmetric positive-definite $Q$. Define $V(\statevec) = \statevec^\top P \statevec$. Then $V > 0$ for $\statevec \ne 0$ (since $P \succ 0$), and along trajectories $\ddt \statevec = \statemat \statevec$:

So $V$ is a strict Lyapunov function and trajectories satisfy $V(\statevec(t)) \le V(\statevec(0)) e^{-\alpha t}$ for $\alpha > 0$ determined by the smallest eigenvalue of $Q$ over the largest of $P$. Hence $\statevec(t) \to 0$ and the origin is asymptotically stable.

Reverse direction (asymptotic stability ⇒ Lyapunov equation has p.d. solution): Suppose $\statemat$ is asymptotically stable. Define

The integral converges absolutely because $e^{\statemat t}$ decays exponentially. $P$ is symmetric (transpose preserves the integrand structure) and positive-definite (since $Q \succ 0$ and $e^{\statemat t}$ is invertible). Direct computation:

using the decay of $e^{\statemat t}$ at $\infty$. Uniqueness of $P$ follows from the linearity of the Lyapunov operator $\mathcal{L}(P) := \statemat^\top P + P \statemat$ and the fact that its kernel is trivial when $\statemat$ has no purely imaginary eigenvalues. ∎

Solution to Exercise 2.6

(⇒) BIBO-stable implies poles in open LHP. Suppose the system is BIBO-stable, i.e. $h \in L^1$. The Laplace transform $H(s) = \int_0^\infty h(t) e^{-st} dt$ then exists and is analytic for $\Re(s) \ge 0$ (since the integrand is bounded by $\abs{h(t)}$ which is integrable, the integral converges uniformly). $H$ is rational (the matrix-inverse formula gives $H(s) = \outputmat (sI - \statemat)^{-1} \inputmat$, a ratio of polynomials with denominator $\det(sI - \statemat)$). A rational function analytic on the closed right half-plane has no poles there — so all poles of $H(s)$ lie in the open LHP.

(⇐) Poles in open LHP implies BIBO-stable. Suppose all poles of $H(s)$ lie in the open LHP. Partial-fraction decomposition gives $H(s) = \sum_k \frac{c_k}{(s - p_k)^{m_k}}$ where $p_k$ are the poles with multiplicity $m_k$, and $\Re(p_k) < 0$ for all $k$. The inverse Laplace transform of each term is $\frac{c_k t^{m_k - 1}}{(m_k - 1)!} e^{p_k t}$, which is $L^1$ on $[0, \infty)$ because the exponential decay beats any polynomial growth in $t$. Sum of $L^1$ functions is $L^1$, so $h \in L^1$ and BIBO holds.

(For the multi-input multi-output case, apply the same argument entrywise to $H(s) \in \C^{M \times D}$.) ∎

2.10 Companion code

Two runnable companions in companions/ch02/jax/:

• lyapunov_qr.py — implements the QR-based Lyapunov-exponent algorithm and applies it to the ring-of-oscillators system from Chapter 1
• stability_regions.py — plots the regions of absolute stability for forward Euler, backward Euler, bilinear (trapezoidal), and ZOH in the complex plane

To run from the repo root:

PYTHONPATH=. python companions/ch02/jax/lyapunov_qr.py
PYTHONPATH=. python companions/ch02/jax/stability_regions.py

Figures land in public/figures/ch02/ (referenced from §2.3 and §2.4 above).

Linear algebra for sequence models: structured matrices and conditioning

3.1 Eigenvalue decomposition and Jordan normal form

Chapter 1’s matrix exponential built directly on the spectral structure of $\statemat$. The fundamental theorem is:

A Jordan block of size $k$ for eigenvalue $\lambda$ looks like

i.e. eigenvalue on the diagonal, ones on the superdiagonal, zeros elsewhere. The matrix exponential of a Jordan block is