Chaos, predictability, Koopman, DMD and SINDy

Before this lesson: linear algebra 05 (eigenvalues, complex pairs), linear algebra 08 (the SVD and the pseudoinverse) and linear algebra 12 (least squares); lesson 02 (RK4, stability) and lesson 03 (the tangent-linear model). Helpful: optimization 09 (sparse regression). Time: about 100 minutes, plus the notebook.

---

Why this matters #

Weather is chaotic. Lorenz found sensitive dependence on initial conditions in 1963 with the three-variable model of this lesson, a truncated model of convection. Today, skilful forecasts of day-to-day midlatitude weather reach about 10 days. Zhang et al. (2019) estimate that a 10 times smaller initial-condition error would add at most about 5 days. Lyapunov exponents explain why the gain is so small: they set how fast errors grow. They also explain three things you meet in this track:

• gradients through long rollouts explode (lesson 03);
• 4D-Var assimilation windows must be short (lesson 06);
• ML weather models are judged by ensembles and statistics at long lead times, not by one trajectory (lesson 16).

The second half is about learning dynamics from data. The Koopman operator turns nonlinear dynamics into linear dynamics on functions of the state. DMD is the best-fit linear model of snapshot data, and it approximates the Koopman operator. SINDy finds sparse, readable equations. These are standard tools in fluid dynamics and climate data analysis. They are also the ideas behind Koopman autoencoders and the linear latent dynamics of state-space models in ML.

1. Predict first #

2. The picture #

Errors that double: the simplest chaotic map #

Iterate a map $x_{k+1} = f(x_k)$ from two nearby starting points. Their distance $\delta_k$ obeys $\delta_{k+1} \approx |f'(x_k)|\,\delta_k$ (a first-order Taylor step). After $k$ steps,

The number $\lambda$ is the Lyapunov exponent (per step): the average logarithmic stretching along the orbit. It is an average, not the local rate at one point.

Example. The doubling map $x \mapsto 2x \bmod 1$ has $f' = 2$ everywhere, so $\lambda = \ln 2$. Every step doubles the error. In binary, one step shifts all digits of $x$ one place to the left, so each step reveals one more binary digit of $x_0$. In float64, the orbit of $x_0 = 0.1$ becomes exactly 0 after 55 steps: all the information in the starting value has been used up. The logistic map $x \mapsto 4x(1 - x)$ also has $\lambda = \ln 2$ (the notebook checks it).

The predictability time. An error $\delta_0$ reaches a tolerance $\Delta$ when $\delta_0 e^{\lambda k} = \Delta$, that is, after

The starting error enters only through a logarithm. A better initial condition helps, but only a little.

Fixed points and their type #

Now a flow $\dot u = f(u)$, with $u \in \R^n$ (we use $u$ for the state, as in lessons 02 and 03). A fixed point $u^*$ has $f(u^*) = 0$: a state that does not move. Near it, write $u = u^* + \delta$. A Taylor step gives $\dot\delta \approx J(u^*)\,\delta$, with the Jacobian $J = \partial f/\partial u$. This is a linear ODE, and its solutions are combinations of $e^{\lambda_i t}v_i$ for the eigenpairs $(\lambda_i, v_i)$ of $J(u^*)$ (linear algebra 05). So:

• the real part of an eigenvalue decides growth ($> 0$) or decay ($< 0$) along its direction;
• the imaginary part decides rotation: a complex pair $a \pm bi$ gives a spiral with angular frequency $b$.

In the plane, the two eigenvalues are fixed by the trace $T = \operatorname{tr} J$ and the determinant $D = \det J$: $\lambda_{1,2} = T/2 \pm \sqrt{T^2/4 - D}$. So one picture classifies every case.

Interactive figure: loading…

Eigenvalues of $J(u^*)$	Type	In the $(T, D)$ plane
real, both negative	stable node	$D > 0$, $T < 0$, below the parabola
real, opposite signs	saddle	$D < 0$
complex pair, negative real part	stable spiral	$D > 0$, $T < 0$, above the parabola
complex pair, positive real part	unstable spiral	$D > 0$, $T > 0$, above the parabola
purely imaginary pair	centre (for the linear system)	$T = 0$, $D > 0$

In three dimensions there are three eigenvalues, and the same rules apply direction by direction. A real negative eigenvalue together with a complex pair of positive real part is a saddle-focus: trajectories come in along one direction and spiral out in a plane.

Two runs of a chaotic flow #

Lorenz-63 is the classic chaotic flow:

Every fixed point of it is unstable (section 3), yet trajectories stay bounded. They wander around two "wings", the Lorenz attractor, and never repeat. Start two trajectories very close together and watch.

Interactive figure: loading…

Chaos is stretching and folding at the same time. Stretching: some directions grow, so nearby states separate. Folding: the attractor is bounded, so the stretched set must fold back onto itself. Together they mix states: after a while, the state could be anywhere on the attractor.

Three ways to look at dynamics #

• States and nearby trajectories: fixed points, linearisation, Lyapunov exponents. This is the first half of the lesson.
• Observables: instead of moving the state, move functions $g(u)$ of the state, for example the temperature at one place. The rule that moves them, the Koopman operator, is linear even when $f$ is not. The price: it acts on a space of functions, which is infinite-dimensional.
• Data: DMD fits the best linear map between snapshots of chosen observables. SINDy fits a sparse combination of candidate terms to the vector field $f$.

The next figure shows the observable view on a nonlinear system with a known answer. In the state $(x_1, x_2)$ the trajectory bends, so no matrix can follow it. With one extra observable, $x_1^2$, the same motion is exactly linear.

Interactive figure: loading…

3. Definitions and main results #

Fixed points and Lorenz-63 #

A fixed point $u^*$ is hyperbolic if no eigenvalue of $J(u^*)$ has zero real part. Near a hyperbolic fixed point, the nonlinear flow looks like the linear flow $\dot\delta = J(u^*)\delta$ (the Hartman–Grobman theorem): stable if all real parts are negative, a saddle if the signs are mixed, with spiralling where there are complex pairs. A centre is not hyperbolic, and there the linearisation can mislead: small nonlinear terms decide.

The Jacobian of Lorenz-63 is

Fixed points. $\dot x = 0$ gives $y = x$. Then $\dot z = 0$ gives $z = x^2/\beta$, and $\dot y = 0$ gives $x(\rho - 1 - z) = 0$. So either $x = 0$ (the origin), or $z = \rho - 1$ and $x^2 = \beta(\rho - 1)$:

• At the origin, the eigenvalues are $11.828$, $-22.828$ and $-8/3$: a saddle.
• At $C_\pm$, they are $-13.855$ and $0.094 \pm 10.195i$: a saddle-focus (unstable spiral), because $\rho = 28$ is above the Hopf value $\rho_H = \sigma(\sigma + \beta + 3)/(\sigma - \beta - 1) = 24.74$.

Every fixed point is unstable, yet trajectories stay bounded: they wander around $C_+$ and $C_-$ on the attractor. The trace of the Jacobian is $-(\sigma + 1 + \beta) = -13.667$ at every point. So volumes in state space shrink by the factor $e^{-13.667}$ per time unit (this follows from Liouville's formula, step 2 of the proof of Theorem 1 below).

Lyapunov exponents #

A small perturbation $\delta$ of a trajectory $u(t)$ obeys the tangent-linear model $\dot\delta = J(u(t))\,\delta$ (lesson 03). Collect all directions at once in the tangent propagator $M(t) = \partial u(t)/\partial u(0)$:

Definition. The Lyapunov exponents $\lambda_1 \ge \dots \ge \lambda_n$ are

where $s_1 \ge \dots \ge s_n$ are the singular values of $M(t)$ (linear algebra 08). $M(t)$ maps a small ball to an ellipsoid, and $s_i$ are its semi-axes, so $\lambda_i$ is the average exponential growth rate of the $i$-th axis.

Oseledets' theorem says that these limits exist and are the same for almost every starting point on the attractor (with respect to an ergodic invariant measure, here the natural measure of the attractor). For almost every $\delta(0)$, $\tfrac1t\ln\big(\|\delta(t)\|/\|\delta(0)\|\big) \to \lambda_1$: a random perturbation turns towards the fastest-growing direction. A bounded system with $\lambda_1 > 0$ is chaotic. For Lorenz-63 the exponents are $(0.906,\ 0,\ -14.572)$ (Sprott, 2003).

Theorem 1 (sum of the exponents). $\displaystyle\sum_i\lambda_i = \lim_{t\to\infty}\frac1t\int_0^t \operatorname{tr}J(u(s))\,ds$.

Why this is true, in one sentence: $\operatorname{tr}J$ is the rate at which small volumes grow, and the volume of the ellipsoid is the product of its axes, $\prod_i s_i$, which grows at the rate $\sum_i\lambda_i$.

Proof.

1. Jacobi's formula: $\frac{d}{dt}\det M = \det M\cdot\operatorname{tr}(M^{-1}\dot M)$.
2. $M^{-1}\dot M = M^{-1}JM$ has the same trace as $J$ (similar matrices). So $\frac{d}{dt}\det M = \operatorname{tr}(J)\det M$, and $\det M(t) = \exp\int_0^t\operatorname{tr}J(u(s))\,ds$ (Liouville's formula).
3. $|\det M| = \prod_i s_i(M)$. Take the logarithm and divide by $t$: $\sum_i\frac1t\ln s_i(M(t)) = \frac1t\int_0^t\operatorname{tr}J\,ds$.
4. Let $t \to \infty$. $\blacksquare$

For Lorenz-63, $\operatorname{tr}J$ is constant, so $0.906 + 0 - 14.572 = -13.667$ exactly. The strongly negative $\lambda_3$ squeezes volumes, so the attractor has zero volume. Its Kaplan–Yorke dimension $2 + \lambda_1/|\lambda_3| = 2.06$ says that it is a thin set, a little more than two-dimensional.

Theorem 2 (a flow has a zero exponent). Let $u(t)$ be a bounded trajectory of $\dot u = f(u)$ that does not converge to a fixed point, and assume that $\|f(u(t))\|$ does not decay to zero exponentially fast. Then one Lyapunov exponent is $0$. The reason: pushing the state along the trajectory only shifts it in time, and a time shift neither grows nor shrinks on average.

Proof of Theorem 2

The assumption means $\frac1t\ln\|f(u(t))\| \to 0$. It holds when the trajectory keeps coming back to a region where $\|f\| \ge c > 0$, away from the fixed points.

1. Differentiate $f(u(t))$ along the trajectory: $\frac{d}{dt}f(u(t)) = J(u(t))\,\dot u(t) = J(u(t))\,f(u(t))$.
2. So $\delta(t) = f(u(t))$ solves the tangent-linear model, with $\delta(0) = f(u(0))$. It is the perturbation that shifts the trajectory in time: $u(t + \varepsilon) \approx u(t) + \varepsilon f(u(t))$.
3. $\|f(u(t))\|$ is bounded above (continuous $f$ on a bounded set), and by the assumption it does not go to zero exponentially. So $\frac1t\ln\|f(u(t))\| \to 0$: the growth rate in this direction is $0$.
4. By Oseledets' theorem, every direction grows at the rate of one of the exponents. So one $\lambda_i$ equals $0$. $\blacksquare$

Computing $\lambda_1$ (Benettin's method). Integrate the state together with one tangent vector $\delta$ (the tangent-linear model, not a second trajectory). Every $\tau$ time units:

1. add $\ln\|\delta\|$ to a running sum;
2. rescale $\delta$ to norm 1 (keep its direction).

After $K$ intervals, $\lambda_1 \approx \frac{1}{K\tau}\sum_k\ln\|\delta_k\|$. The rescaling is needed: $\delta$ grows like $e^{0.906t}$ and would overflow float64 near $t \approx 780$.

All exponents, and why two trajectories need care

To get all $n$ exponents, carry $n$ tangent vectors as the columns of $Q$. Every step, propagate them and re-orthonormalise with a QR factorisation $Q_{\text{new}}R = M_\tau Q$ (linear algebra 12). Then $\lambda_i \approx \frac{1}{K\tau}\sum_k\ln|R_{ii}^{(k)}|$. Without QR, every vector would turn towards the fastest-growing direction and all estimates would approach $\lambda_1$.

The two-trajectory version (as in the figure) must also rescale the separation before it stops being small; otherwise the estimate saturates at the size of the attractor. A single pair also gives a noisy slope, because the local growth rate changes along the attractor.

Predictability. An error grows like $\delta_0 e^{\lambda_1 t}$. For Lorenz-63:

• the Lyapunov time is $1/\lambda_1 = 1.10$ time units (growth by the factor $e$);
• the doubling time is $\ln 2/\lambda_1 = 0.77$ time units;
• the time to reach a tolerance $\Delta$ is the predictability time $T_p = \frac{1}{\lambda_1}\ln(\Delta/\delta_0)$.

Three consequences for this track:

• Training through rollouts. The gradient of a loss at time $T$ with respect to $u(0)$ contains $M(T)^\top$, so it grows like $e^{\lambda_1 T}$ (lesson 03).
• Data assimilation. The 4D-Var cost function gets many local minima when the window is long, so windows are short (lesson 06).
• Forecasts depend on the day. The finite-time growth rate varies along the attractor, so predictability depends on the state of the flow. This is why ensemble forecasts show more spread on some days than on others.

The Koopman operator #

Take a discrete-time system $u_{k+1} = F(u_k)$. For a flow, $F$ is the flow map over a time step $\Delta t$. An observable is a function $g$ of the state, for example the temperature at one point.

Definition. The Koopman operator $\mathcal K$ acts on observables by composition with $F$:

In words: $\mathcal K g$ is the observable "$g$, one step later".

Proposition (linearity and eigenfunctions).

1. $\mathcal K$ is linear: $\mathcal K(ag + bh) = a\,\mathcal Kg + b\,\mathcal Kh$ for any $F$, linear or not.
2. If $\mathcal K\varphi = \mu\varphi$ (a Koopman eigenfunction with eigenvalue $\mu$), then $\varphi(u_k) = \mu^k\varphi(u_0)$: along every trajectory, $\varphi$ evolves exactly linearly.
3. If $\varphi_1, \varphi_2$ are eigenfunctions with eigenvalues $\mu_1, \mu_2$, then $\varphi_1\varphi_2$ is an eigenfunction with eigenvalue $\mu_1\mu_2$.

Proof of the proposition

1. $(ag + bh)(F(u)) = a\,g(F(u)) + b\,h(F(u))$. Nothing here uses linearity of $F$: we only evaluate functions at the point $F(u)$.
2. $\varphi(u_k) = (\mathcal K\varphi)(u_{k-1}) = \mu\,\varphi(u_{k-1})$. Repeat $k$ times.
3. $\mathcal K(\varphi_1\varphi_2)(u) = \varphi_1(F(u))\,\varphi_2(F(u)) = \mu_1\varphi_1(u)\,\mu_2\varphi_2(u)$. $\blacksquare$

The price: $\mathcal K$ acts on a space of functions, which is infinite-dimensional. In continuous time, an eigenfunction with $f\cdot\nabla\varphi = \omega\varphi$ satisfies $\varphi(u(t)) = e^{\omega t}\varphi(u(0))$, and the flow map over $\Delta t$ has the eigenvalue $\mu = e^{\omega\Delta t}$.

Linear systems. For $F(u) = Au$ and a left eigenvector $w$ of $A$ ($w^\top A = \mu w^\top$), the observable $\varphi(u) = w^\top u$ is an eigenfunction: $\varphi(Au) = w^\top Au = \mu\,w^\top u$.

A finite invariant subspace for a nonlinear system (Brunton and Kutz, §7.4). This is the system of the figure in section 2. With continuous-time rates $\mu$ and $\lambda$ (here $\lambda$ is a rate of the model, not a Lyapunov exponent):

Take the observables $y = (x_1, x_2, x_1^2)$. Step by step:

• $\dot y_1 = \dot x_1 = \mu y_1$;
• $\dot y_2 = \lambda x_2 - \lambda x_1^2 = \lambda y_2 - \lambda y_3$;
• $\dot y_3 = 2x_1\dot x_1 = 2\mu x_1^2 = 2\mu y_3$ (chain rule).

So

a linear system, exactly, and $y(t + \Delta t) = e^{K\Delta t}y(t)$. The nonlinear term $x_1^2$ is now a coordinate. The eigenvalues of $K$ are $\mu$, $\lambda$ and $2\mu$, with the eigenfunctions $x_1$, $x_2 - b\,x_1^2$ with $b = \lambda/(\lambda - 2\mu)$, and $x_1^2$.

Check of the middle one: $\frac{d}{dt}(x_2 - bx_1^2) = \lambda x_2 - \lambda x_1^2 - 2b\mu x_1^2 = \lambda x_2 - (\lambda + 2b\mu)x_1^2$. This equals $\lambda(x_2 - bx_1^2)$ exactly when $\lambda + 2b\mu = \lambda b$, that is, $b(\lambda - 2\mu) = \lambda$. The curve $x_2 = b x_1^2$, where this eigenfunction is 0, is the slow manifold of the figure: the fast rate $\lambda$ pulls every trajectory onto it.

When no finite subspace exists. The quiz shows the first obstacle: for the logistic map the degree of the polynomials doubles at each step. For chaos the obstacle is deeper.

Mixing and the continuous Koopman spectrum

The Lorenz attractor is mixing (Luzzatto, Melbourne and Paccaut, 2005): any two regions of the attractor become, after a long time, statistically independent. On a mixing attractor the only eigenfunctions of $\mathcal K$ in $L^2$ of the invariant measure are the constants, and the rest of the spectrum is continuous. A finite linear model of chaos can therefore only be approximate, and only over short horizons. Practical methods (extended DMD with many observables, delay coordinates, Koopman autoencoders) give good short-term predictions and useful coherent patterns, not an exact finite linear system.

Koopman modes. Suppose each state component is a combination of eigenfunctions, $u = \sum_j v_j\varphi_j(u)$. Then $u_k = \sum_j v_j\,\mu_j^k\,\varphi_j(u_0)$: the state is a sum of fixed spatial patterns $v_j$ (the Koopman modes), each with one growth rate and one frequency. DMD estimates the pairs $(v_j, \mu_j)$ from data, with the state coordinates as observables.

Dynamic mode decomposition (DMD) #

Data. Snapshot pairs $(x_k, x_k')$ with $x_k' = F(x_k)$, usually $x_k' = x_{k+1}$ at spacing $\Delta t$. (From here on we write $x$ for a data vector, as in the DMD literature.) Arrange them as columns:

with $n$ the number of measured values per snapshot (for example grid points) and $m$ the number of snapshots.

The best-fit linear map. $A = \operatorname*{argmin}_A\|X' - AX\|_F = X'X^+$, with $X^+$ the pseudoinverse (linear algebra 08). Each row of $A$ is a least-squares problem, and $X^+$ gives the minimum-norm solution. For a field with $n = 10^6$ grid points, $A$ has $10^{12}$ entries. Exact DMD never forms it.

Exact DMD (Tu et al., 2014, Algorithm 2):

1. Reduced SVD $X = U\Sigma V^*$, truncated to rank $r$: $U_r$ ($n \times r$), $\Sigma_r$, $V_r$ ($m \times r$).
2. Projected operator $\tilde A = U_r^*X'V_r\Sigma_r^{-1}$ ($r \times r$). It equals $U_r^*AU_r$ for $A = X'V_r\Sigma_r^{-1}U_r^*$: the map $A$ seen in the coordinates of the top $r$ singular vectors.
3. Eigendecomposition $\tilde AW = W\Lambda$. The diagonal of $\Lambda$ holds the DMD eigenvalues $\mu_j$.
4. Exact DMD modes $\Phi = X'V_r\Sigma_r^{-1}W$ ($n \times r$). A mode is defined only up to a scale factor (Tu et al. divide column $j$ by $\mu_j$; Brunton and Kutz do not).
5. Continuous-time rates $\omega_j = \ln(\mu_j)/\Delta t$: growth rate $\operatorname{Re}\omega_j$ and angular frequency $\operatorname{Im}\omega_j$. Amplitudes $b = \Phi^+x_1$, and the model $x(t) \approx \sum_j\phi_j\,e^{\omega_j t}\,b_j$, where $\phi_j$ is column $j$ of $\Phi$.

The cost is one thin SVD, $O(nm^2)$, plus small $r \times r$ problems.

Reading a DMD eigenvalue. Write $\mu = |\mu|e^{i\theta}$. Then $\omega = (\ln|\mu| + i\theta)/\Delta t$:

• $|\mu| < 1$: the mode decays; $|\mu| = 1$: it neither grows nor decays; $|\mu| > 1$: it grows;
• the angle $\theta$ is the rotation per step, so the angular frequency is $\theta/\Delta t$.

Interactive figure: loading…

Theorem 3 (exact DMD modes are eigenvectors of $A$; Tu et al., Theorem 1). Let $A = X'V_r\Sigma_r^{-1}U_r^*$. If $\tilde Aw = \mu w$ with $\mu \ne 0$, then $\phi = X'V_r\Sigma_r^{-1}w$ is nonzero and $A\phi = \mu\phi$. Conversely, every nonzero eigenvalue of $A$ is an eigenvalue of $\tilde A$.

Proof of Theorem 3

Write $B = X'V_r\Sigma_r^{-1}$, so that $A = BU_r^*$ and $\tilde A = U_r^*B$. (The two products of the same two matrices in the two orders: this is the whole trick.)

1. $A\phi = BU_r^*Bw = B\tilde Aw = \mu Bw = \mu\phi$.
2. If $\phi = Bw = 0$, then $\tilde Aw = U_r^*Bw = 0$, so $\mu = 0$, a contradiction. So $\phi \ne 0$.
3. Conversely, let $A\phi = \mu\phi$ with $\mu \ne 0$, and $w = U_r^*\phi$. Then $\tilde Aw = U_r^*BU_r^*\phi = U_r^*A\phi = \mu w$.
4. $w \ne 0$: if $U_r^*\phi = 0$, then $A\phi = BU_r^*\phi = 0$, so $\mu = 0$. $\blacksquare$

The older projected modes $U_rW$ are the projections of the exact modes onto the column space of $X$. They are not eigenvectors of $A$ in general.

Theorem 4 (DMD recovers a linear system). If $x_k' = A_{\text{true}}x_k$ for all $k$ and $X$ has rank $n$, then $X'X^+ = A_{\text{true}}$, and exact DMD with $r = n$ returns the eigenvalues of $A_{\text{true}}$.

Proof.

1. $X' = A_{\text{true}}X$, so $X'X^+ = A_{\text{true}}XX^+$.
2. $XX^+$ is the orthogonal projector onto the column space of $X$, which is all of $\R^n$ when $\operatorname{rank}X = n$. So $XX^+ = I$ and $X'X^+ = A_{\text{true}}$.
3. With $r = n$, Theorem 3 gives all nonzero eigenvalues of $A$. ($A_{\text{true}}$ is invertible when it comes from a flow map, so no eigenvalue is 0.) $\blacksquare$

If the data span only a subspace $S$ (for example one trajectory inside an invariant subspace), then $X'X^+ = A_{\text{true}}P_S$, with $P_S$ the projector onto $S$: DMD sees only the modes that the data excite.

DMD in practice.

• Rank $r$. Too large with noisy data gives spurious eigenvalues; too small misses dynamics. Look for a gap in the singular values. Gavish–Donoho hard thresholding is a principled choice.
• Noise. Least squares treats $X$ as exact, so sensor noise biases the eigenvalues towards the origin (too much damping; in the figure, try $\sigma = 0.3$). Forward–backward and total-least-squares DMD reduce this bias (Dawson et al., 2016).
• Few observables. A scalar signal $\cos(\omega k\Delta t)$ needs the two eigenvalues $e^{\pm i\omega\Delta t}$, but a $1 \times 1$ matrix has one. Stack time-delayed copies (a Hankel matrix) to get enough rows: this is time-delay embedding.
• Sampling. $\omega = \ln\mu/\Delta t$ is unique only for $|\operatorname{Im}\omega| < \pi/\Delta t$; faster oscillations alias (they look like slower ones).
• DMD or POD (PCA)? Proper orthogonal decomposition (POD), also called PCA or EOF analysis, ranks spatial patterns by variance and ignores the time order of the snapshots (linear algebra 09). A DMD mode is a pattern with one frequency and one growth rate. So DMD can separate, for example, an annual cycle from a slower oscillation that has a similar variance, where POD would mix them.

SINDy: sparse identification of nonlinear dynamics #

Model. Assume that each component of $f$ is a sparse combination of known candidate functions: $\dot u \approx \Theta(u)\,\Xi$. Here $\Theta(u)$ is a row of $p$ library functions of the state (for the Lorenz state $(x, y, z)$: $1, x, y, z, x^2, xy, \dots$), and $\Xi \in \R^{p \times n}$ is sparse: most entries are 0. With $m$ samples, stack the states as rows of $X$ ($m \times n$), the time derivatives as rows of $\dot X$, and the library as $\Theta(X)$ ($m \times p$):

Column $k$ of $\Xi$ is a regression problem for $\dot u_k$. We want few nonzero coefficients, so that the result is a short equation. This is like the lasso, $\min\|\dot X_k - \Theta\xi\|^2 + \alpha\|\xi\|_1$ (optimization 09), but SINDy usually uses a simpler method.

Sequentially thresholded least squares (STLSQ) (Brunton, Proctor and Kutz, 2016):

1. $\Xi \leftarrow$ the least-squares solution of $\Theta(X)\Xi = \dot X$ (dense: every term is used).
2. Set every entry with $|\xi_{jk}| < \eta$ to zero. The threshold $\eta$ is called $\lambda$ in the book.
3. For each column $k$, refit by least squares using only the library terms that are still active.
4. Repeat steps 2–3 until the set of active terms no longer changes.

Proposition (STLSQ stops). The active set of each column never grows, so STLSQ reaches a fixed point after at most $p$ changes per column.

Proof. An entry set to zero in step 2 is excluded from the refit in step 3, so it stays zero: entries can leave the active set but never return. Each change removes at least one of the $p$ terms. When no entry falls below $\eta$, steps 2–3 change nothing. $\blacksquare$

For Lorenz-63 with all monomials up to degree 2 ($p = 10$), $\Xi$ has 30 entries and 7 of them are nonzero: $-10, 10$ (for $\dot x$: $x$, $y$); $28, -1, -1$ (for $\dot y$: $x$, $y$, $xz$); $-8/3, 1$ (for $\dot z$: $z$, $xy$).

Interactive figure: loading…

Choices that decide success:

• The library must contain the true terms, in good coordinates. SINDy cannot invent a term that is not in the list.
• The threshold must lie between the largest spurious coefficient and the smallest true one (here $|-1| = 1$). Normalise the library columns, or the threshold means different things for $x$ and $x^2$.
• Derivatives. Finite differences amplify noise by about $1/\Delta t$. Use smoothed or total-variation-regularised derivatives, or the weak (integral) form of SINDy.
• Relation to DMD: with a linear library, a discrete-time target $x_{k+1}$ and no threshold, SINDy is DMD.

Answers to section 1. (1) The separation grows like $10^{-8}e^{0.906t}$. It reaches order 1 after $\ln(10^8)/0.906 \approx 20$ time units, and the size of the attractor (about 20) about 3 time units later. A 1000 times smaller error gains only $\ln(1000)/0.906 \approx 7.6$ time units. (2) Not in the state coordinates: the fixed points of $x \mapsto Ax$ form a linear subspace (the null space of $A - I$), so a linear model cannot have the three isolated fixed points of Lorenz-63. With observables it can work exactly for some systems (the slow-manifold example), but never for a chaotic one, whose Koopman spectrum is continuous.

4. Worked example #

(a) DMD by hand. Data from $A_{\text{true}} = \begin{pmatrix} 0.5 & 0.2\\ 0 & 0.8\end{pmatrix}$ with $x_0 = (1, 1)$ and $\Delta t = 0.1$.

Step 1, the snapshots: $x_1 = A_{\text{true}}x_0 = (0.5 + 0.2,\ 0.8) = (0.7, 0.8)$ and $x_2 = A_{\text{true}}x_1 = (0.35 + 0.16,\ 0.64) = (0.51, 0.64)$.

Step 2, the data matrices (the snapshots $x_0, x_1$ and their successors $x_1, x_2$):

Step 3, the pseudoinverse. $X$ is square with $\det X = 0.8 - 0.7 = 0.1 \ne 0$, so $X^+ = X^{-1}$:

Step 4, the product $X'X^{-1}$, row by row:

This is $A_{\text{true}}$ (Theorem 4). Step 5, the eigenvalues: the matrix is triangular, so $\mu_1 = 0.5$ and $\mu_2 = 0.8$, with the eigenvectors (modes) $(1, 0)$ and $(2, 3)$. The continuous-time rates are $\ln 0.5/0.1 = -6.93$ and $\ln 0.8/0.1 = -2.23$. Both modes decay; the second one more slowly.

(b) Lifting makes DMD exact. In the slow-manifold system take $\mu = -0.05$, $\lambda = -1$ and $\Delta t = 0.1$.

• The Koopman eigenvalues of the flow map are $e^{\mu\Delta t} = 0.99501$, $e^{\lambda\Delta t} = 0.90484$ and $e^{2\mu\Delta t} = 0.99005$.
• $b = \lambda/(\lambda - 2\mu) = -1/(-0.9) = 1.111$.
• On a trajectory from $(1.5, -1)$ over $0 \le t \le 10$, DMD on the state $(x_1, x_2)$ returns the rates $-0.05$ and $-1.22$. The rate $\mu$ is right, because the $x_1$ equation is linear. The rate $-1.22$ is wrong: it is neither $\lambda = -1$ nor $2\mu = -0.1$. A $2 \times 2$ matrix has room for only two rates, so the least-squares fit mixes the effects of $\lambda$ and $2\mu$.
• DMD on $(x_1, x_2, x_1^2)$ returns $-0.05$, $-1$ and $-0.1$, and the fitted matrix equals $e^{K\Delta t}$ to about $10^{-14}$ (section 5).

(c) Lorenz numbers without simulation.

• Sum of the exponents: $\sum_i\lambda_i = \operatorname{tr}J = -(10 + 1 + 8/3) = -13.667$.
• With $\lambda_1 = 0.906$ and $\lambda_2 = 0$: $\lambda_3 = -13.667 - 0.906 = -14.573$ (Sprott: $-14.5723$).
• Kaplan–Yorke dimension: $2 + 0.906/14.572 = 2.062$.
• Doubling time: $0.693/0.906 = 0.77$ time units.

(d) One STLSQ step by hand. The first, dense least-squares fit in section 5 (the array dense) gives, for $\dot y$, the coefficients $27.991$ ($x$), $-0.998$ ($y$), $-1.000$ ($xz$), and every other coefficient below $0.005$ in absolute value.

• With $\eta = 0.1$, step 2 keeps $\{x, y, xz\}$ and deletes the other 7 terms.
• Step 3 refits $\dot y$ on these three columns only.
• The next step 2 deletes nothing, so STLSQ stops.

Any $\eta$ between about $0.005$ and $1$ gives the right active set here; $\eta = 2$ deletes the true terms $-y$ and $-xz$ (and $xy$ in $\dot z$).

5. In code #

The first block finds the Lorenz fixed points and estimates $\lambda_1$ with Benettin's method:

import numpy as np

sigma, rho, beta = 10.0, 28.0, 8 / 3

def f(u):                                   # Lorenz-63
    x, y, z = u
    return np.array([sigma * (y - x), x * (rho - z) - y, x * y - beta * z])

def jac(u):                                 # its Jacobian J = df/du
    x, y, z = u
    return np.array([[-sigma, sigma, 0.0], [rho - z, -1.0, -x], [y, x, -beta]])

c = np.sqrt(beta * (rho - 1))               # C± = (±c, ±c, ρ - 1)
for p in [(0, 0, 0), (c, c, rho - 1)]:
    print(np.round(np.sort_complex(np.linalg.eigvals(jac(np.array(p)))), 3))
# [-22.828+0.j  -2.667+0.j  11.828+0.j]             saddle at the origin
# [-13.855 +0.j  0.094-10.195j  0.094+10.195j]      saddle-focus at C+
print(np.trace(jac(np.zeros(3))))           # -13.666666666666666 = -(σ + 1 + β)

def rk4(F, u, h):
    k1 = F(u); k2 = F(u + h / 2 * k1); k3 = F(u + h / 2 * k2); k4 = F(u + h * k3)
    return u + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4)

# Benettin: state + one tangent vector d (d' = J(u) d), renormalised every step
G = lambda w: np.concatenate([f(w[:3]), jac(w[:3]) @ w[3:]])
h, u = 0.01, np.array([1.0, 1.0, 1.0])
for _ in range(2000):                       # 20 time units of transient
    u = rk4(f, u, h)
w, total, n = np.concatenate([u, [1.0, 0.0, 0.0]]), 0.0, 50_000
for _ in range(n):
    w = rk4(G, w, h)
    g = np.linalg.norm(w[3:]); total += np.log(g); w[3:] /= g
lam1 = total / (n * h)
print(f"{lam1:.3f}  doubling time {np.log(2) / lam1:.2f}")   # 0.908  doubling time 0.76

Over 500 time units the estimate is $0.908$; longer runs converge to $0.906$. The second block is exact DMD from the SVD, on the worked examples (a) and (b):

import numpy as np
from scipy.linalg import expm

def exact_dmd(X, Xp, r):
    """Exact DMD (Tu et al. 2014): eigenvalues mu and modes Phi of A = X' X^+, rank r."""
    U, s, Vh = np.linalg.svd(X, full_matrices=False)
    U, s, V = U[:, :r], s[:r], Vh[:r].conj().T
    At = U.conj().T @ Xp @ V / s            # r x r projected operator
    mu, W = np.linalg.eig(At)
    Phi = Xp @ V / s @ W                    # exact DMD modes
    return mu, Phi

# (a) the 2x2 worked example: DMD returns the matrix that made the data
A_true = np.array([[0.5, 0.2], [0.0, 0.8]])
snaps = [np.array([1.0, 1.0])]
for _ in range(2):
    snaps.append(A_true @ snaps[-1])
S = np.column_stack(snaps)
mu, Phi = exact_dmd(S[:, :-1], S[:, 1:], r=2)
print(np.sort(mu.real), np.sort(np.log(mu.real) / 0.1).round(3))   # [0.5 0.8] [-6.931 -2.231]
print(np.abs(A_true @ Phi - Phi * mu).max() < 1e-12)              # True: A φ = μ φ

# (b) the slow-manifold system, sampled exactly: state vs lifted observables
mu_, lam, dt = -0.05, -1.0, 0.1
t = np.arange(0, 10 + 1e-9, dt)
x1 = 1.5 * np.exp(mu_ * t)
b = lam / (lam - 2 * mu_)
x2 = (-1.0 - b * 1.5**2) * np.exp(lam * t) + b * x1**2   # exact solution from (1.5, -1)
rates = lambda Z, r: np.sort(np.log(exact_dmd(Z[:, :-1], Z[:, 1:], r)[0]).real / dt)
print(rates(np.vstack([x1, x2]), 2).round(3))            # [-1.222 -0.05 ]: one rate is wrong
lifted = np.vstack([x1, x2, x1**2])
print(rates(lifted, 3).round(3))                         # [-1.   -0.1  -0.05]: λ, 2μ, μ
K = np.array([[mu_, 0, 0], [0, lam, -lam], [0, 0, 2 * mu_]])
A_fit = lifted[:, 1:] @ np.linalg.pinv(lifted[:, :-1])
print(f"{np.abs(A_fit - expm(K * dt)).max():.0e}")       # 3e-15

The third block is SINDy on Lorenz-63, with derivatives estimated from the data:

import numpy as np
from scipy.integrate import solve_ivp

rhs = lambda t, u: [10 * (u[1] - u[0]), u[0] * (28 - u[2]) - u[1], u[0] * u[1] - 8 / 3 * u[2]]
dt = 0.002
t = np.arange(0, 20, dt)
X = solve_ivp(rhs, (0, t[-1]), [-8.0, 8.0, 27.0], t_eval=t, method="DOP853",
              rtol=1e-12, atol=1e-12).y.T               # 10,000 states, shape (10000, 3)
dX = np.gradient(X, dt, axis=0, edge_order=2)           # derivatives from the data

x, y, z = X.T
names = ["1", "x", "y", "z", "x^2", "xy", "xz", "y^2", "yz", "z^2"]
Theta = np.column_stack([np.ones_like(x), x, y, z, x*x, x*y, x*z, y*y, y*z, z*z])

def stlsq(Theta, dX, eta, iters=10):
    Xi = np.linalg.lstsq(Theta, dX, rcond=None)[0]      # step 1: dense least squares
    for _ in range(iters):
        small = np.abs(Xi) < eta                         # step 2: threshold
        Xi[small] = 0
        for k in range(dX.shape[1]):                     # step 3: refit the active terms
            big = ~small[:, k]
            Xi[big, k] = np.linalg.lstsq(Theta[:, big], dX[:, k], rcond=None)[0]
    return Xi

dense = np.linalg.lstsq(Theta, dX, rcond=None)[0]
true = np.zeros((10, 3))
true[[1, 2], 0] = [-10, 10]; true[[1, 2, 6], 1] = [28, -1, -1]; true[[3, 5], 2] = [-8 / 3, 1]
print(f"largest spurious |xi| in the dense fit: {np.abs(dense[true == 0]).max():.4f}")  # 0.0040
Xi = stlsq(Theta, dX, eta=0.1)
for name, row in zip(names, Xi):
    if np.any(row):
        print(f"{name:>4}", "  ".join(f"{v:9.4f}" for v in row))
#    x    -9.9992    27.9924     0.0000
#    y     9.9992    -0.9986     0.0000
#    z     0.0000     0.0000    -2.6664
#   xy     0.0000     0.0000     0.9999
#   xz     0.0000    -0.9998     0.0000
rel = np.abs(Xi[true != 0] - true[true != 0]) / np.abs(true[true != 0])
print(f"same terms: {np.array_equal(Xi != 0, true != 0)}, largest relative error {rel.max():.2%}")
# same terms: True, largest relative error 0.14%

The dense fit already has the seven true coefficients within 0.25% (the error comes from the second-order differences), and every spurious one is below 0.005. STLSQ removes the spurious terms, and the refit brings the largest relative error to 0.14%. PyDMD and PySINDy are the standard libraries for these methods; the notebook builds them from scratch.

6. Where it shows up #

• Weather and climate. Lorenz-63 and Lorenz-96 are the test beds of data assimilation (lesson 06). Lyapunov exponents explain the predictability limit, the need for ensembles, and why a perfect emulator still diverges from the truth. So long rollouts are evaluated by statistics, spectra and probabilistic scores (CRPS, lesson 07), not by pointwise error at long lead times.
• Climate patterns. DMD-type methods extract oscillating modes from geophysical fields: multi-resolution DMD on sea-surface temperature data from 1990–2010 isolates the 1997 El Niño as its own mode (Kutz, Fu and Brunton, 2016).
• Seismology. Simple spring–slider models of a fault with friction are chaotic (Huang and Turcotte, 1990). This is one argument for why the exact timing of earthquakes is hard to predict, and why seismic hazard is stated as a probability.
• Training through dynamics. Backpropagation through a rollout of length $T$ of a chaotic model multiplies gradients by about $e^{\lambda_1 T}$ (lesson 03). Short rollouts keep these gradients and the memory cost under control: GraphCast trains with rollouts of at most 12 six-hour steps (Lam et al., 2023).
• Machine learning. Koopman autoencoders learn coordinates in which the dynamics are linear (Lusch, Kutz and Brunton, 2018); SINDy autoencoders learn coordinates in which the dynamics are sparse (Champion et al., 2019). State-space layers (S4, Mamba) also evolve a latent state by a linear recurrence and read it out nonlinearly (deep learning 16).
• Equation discovery and closures. SINDy and its PDE version identify reduced models of fluid and plasma flows and subgrid closure terms.

7. Common mistakes #

• Estimating $\lambda_1$ from two trajectories without rescaling, or over a too short time. The separation saturates at the attractor size, and finite-time rates fluctuate along the attractor.
• Mixing units. An exponent per sampling step equals $\lambda\Delta t$, where $\lambda$ is the exponent per unit time. In the same way, a DMD eigenvalue $\mu$ is per step and $\omega = \ln\mu/\Delta t$ is per unit time.
• Reading $\operatorname{Re}\mu$ as growth. Growth or decay is decided by $|\mu|$ compared with 1 (or by $\operatorname{Re}\omega$), not by the sign of $\operatorname{Re}\mu$.
• Too coarse sampling for DMD. Frequencies above $\pi/\Delta t$ alias.
• Calling the projected modes $U_rW$ exact DMD modes, or choosing the rank $r$ by eye on noisy data.
• Believing that a finite linear model can be exact for chaos. Truncated Koopman models of chaotic systems are short-horizon approximations.
• SINDy: a library without the true terms; a threshold above the smallest true coefficient; unnormalised columns; raw finite differences of noisy data.
• Reading a DMD eigenvalue near 1 as dynamics when it is just the mean of the data (a constant observable has $\mu = 1$). Subtract the mean first, or keep this in mind.

8. Check yourself #

1. For Lorenz-63, $\lambda_1 = 0.906$ and $\lambda_2 = 0$. Find $\lambda_3$ without any simulation.
   Answer

   Theorem 1 with the constant trace: $\lambda_1 + \lambda_2 + \lambda_3 = -(\sigma + 1 + \beta) = -13.667$, so $\lambda_3 = -14.573$.
2. Suppose a single exponent $\lambda$ described weather error growth. What $\lambda$ makes a 10 times smaller initial error add 5 days of skill (the estimate of Zhang et al.)? What is the doubling time?
   Answer

   $T_p = \frac1\lambda\ln(\Delta/\delta_0)$, so dividing $\delta_0$ by 10 adds $\ln(10)/\lambda = 5$ days: $\lambda = 0.46$ per day, and the doubling time is $\ln 2/0.46 = 1.5$ days. (Real error growth is faster at small scales and saturates there, so one exponent is a simplification.)
3. Data come from $x_{k+1} = Ax_k$ with $A = \operatorname{diag}(0.9, 0.5)$ and $x_0 = (1, 0)$. What does DMD return, and why?
   Answer

   Every snapshot is a multiple of $e_1$, so $\operatorname{rank}X = 1$ and $X'X^+ = 0.9\,XX^+ = \operatorname{diag}(0.9, 0)$. DMD finds only $\mu = 0.9$. The mode with $0.5$ is never excited, so the data cannot reveal it (Theorem 4 needs $\operatorname{rank}X = n$).
4. Show that $x_1^2$ is a Koopman eigenfunction of the slow-manifold system, and give its eigenvalue in continuous time and for the flow map over $\Delta t$.
   Answer

   $\frac{d}{dt}x_1^2 = 2x_1\dot x_1 = 2\mu x_1^2$, so $x_1^2(t) = e^{2\mu t}x_1^2(0)$: the continuous-time eigenvalue is $2\mu$, and the flow-map eigenvalue is $e^{2\mu\Delta t}$. It is the product of the eigenfunction $x_1$ with itself, so its eigenvalue is the square of $e^{\mu\Delta t}$.
5. Why does DMD on the scalar signal $x_k = \cos(\omega k\Delta t)$ fail, and how do time delays fix it?
   Answer

   A $1 \times 1$ fit has one real eigenvalue, but $\cos$ needs the pair $e^{\pm i\omega\Delta t}$. The signal satisfies $x_{k+2} = 2\cos(\omega\Delta t)\,x_{k+1} - x_k$, so the delay vector $(x_k, x_{k+1})$ evolves by an exact $2 \times 2$ linear map with eigenvalues $e^{\pm i\omega\Delta t}$. DMD on the Hankel matrix recovers them.
6. You run STLSQ on noisy Lorenz data and get extra terms such as $0.03\,z^2$ in $\dot x$. Name two fixes.
   Answer

   Better derivatives (smoothing, total-variation regularisation, or the weak/integral form); a larger threshold (still below the smallest true coefficient, 1) after normalising the library columns; more data; or a shorter library. Then check the fitted model by simulating it and comparing statistics with the data.

Recall cardsAnswer in your head first, then reveal. Add the cards to your review deck to see them again just before you would forget.

1. Define the largest Lyapunov exponent $\lambda_1$ of a trajectory.

   Show answer

2. What is $\lambda_1$ for Lorenz-63 ($\sigma = 10$, $\rho = 28$, $\beta = 8/3$), and what is the doubling time of an error?

   Show answer

3. What is the sum of all Lyapunov exponents, and why?

   Show answer

4. Give the predictability time for an error $\delta_0$ and a tolerance $\Delta$.

   Show answer

5. Why must Benettin's method renormalise the tangent vector?

   Show answer

6. Define the Koopman operator, and say why it is linear.

   Show answer

7. Give the exact DMD formula for the projected operator and the modes.

   Show answer

8. A DMD eigenvalue is $\mu$ with time step $\Delta t$. What are its growth rate and angular frequency?

   Show answer

9. What problem does SINDy solve?

   Show answer

10. DMD or POD: what is the difference in what a mode means?

    Show answer

Add these cards to my reviewOpen the review →

9. Interview angle #

Typical questions and short model answers:

• "What is a Lyapunov exponent, and what does it mean for forecasting?" — The average exponential rate at which nearby trajectories separate. The predictability time grows only like $\ln(1/\delta_0)/\lambda_1$, so better observations help little; hence ensembles and probabilistic scores.
• "What is the Koopman operator? What is the catch?" — Composition with the dynamics; it is linear on observables. It is infinite-dimensional, and for chaos its spectrum is continuous, so finite models are approximate.
• "DMD or PCA/POD?" — POD ranks patterns by variance; DMD gives patterns with one frequency and one growth rate each, from the best-fit linear map. Exact DMD never forms $A$.
• "How would you discover equations from data, and what breaks it?" — SINDy: a library plus sparse regression. It breaks with noisy derivatives, a library without the right terms, or bad coordinates.
• "Your emulator matches the truth for 5 days and then drifts, although its one-step error is tiny. Is it broken?" — Not necessarily. With $\lambda_1 > 0$ any small error grows. Evaluate statistics, spectra and ensemble scores instead of pointwise error at long lead times.

10. Practice #

The notebook for this lesson runs in Colab or in the browser. It has five exercises, each with a check cell and a folded solution:

1. Find the fixed points of Lorenz-63 and classify them with the eigenvalues of the Jacobian.
2. Estimate the largest Lyapunov exponent with the tangent-linear model and renormalisation (≈ 0.906).
3. Write exact DMD from the SVD and check that it recovers the eigenvalues of a known linear system.
4. Compare DMD on the state and on the lifted observables of the slow-manifold system.
5. Write STLSQ and recover the Lorenz-63 coefficients from data within 1%.

11. Watch and read #

Chaotic Dynamical SystemsSteve Brunton · 44 min

Lorenz-63, sensitive dependence and why chaos limits prediction.

Lyapunov Exponents - Dynamical Systems | Lecture 31Jason Bramburger · 17 min

How the exponents are defined and computed, as in section 3.

Dynamic Mode Decomposition (Overview)Steve Brunton · 18 min

The SVD-based DMD algorithm and how to read the eigenvalues and modes.

Koopman Spectral Analysis (Overview)Steve Brunton · 28 min

Observables, Koopman eigenfunctions and the slow-manifold example.

Sparse Identification of Nonlinear Dynamics (SINDy)Steve Brunton · 27 min

The library, the sparse regression and the Lorenz-63 example.

• Brunton and Kutz, Data-Driven Science and Engineering, 2nd edition, free draft from the authors, chapter 7: §7.1 (Lorenz and the challenges), §7.2 (DMD: the exact DMD algorithm and the spectral expansion), §7.3 (SINDy, with code for STLSQ and the Lorenz example), §7.4 (Koopman operator theory, the slow-manifold example and the logistic map), §7.5 (data-driven Koopman analysis: EDMD and delay coordinates).
• Tu, Rowley, Luchtenburg, Brunton and Kutz, On dynamic mode decomposition: theory and applications (2014): §2.2 (exact DMD, Algorithm 2, Theorems 1–2), §2.3 (projected versus exact modes), §4.1 (Koopman spectral analysis).
• Brunton, Proctor and Kutz, Discovering governing equations from data by sparse identification of nonlinear dynamical systems, PNAS (2016): the original SINDy paper (free arXiv version).
• Strogatz, Nonlinear Dynamics and Chaos (book): chapters 5–6 (linear systems and phase plane, classification of fixed points) and chapter 9 (the Lorenz equations).
• Sprott, Chaos and Time-Series Analysis (Oxford, 2003): the Lyapunov exponents of Lorenz-63 used here.
• Zhang, Sun, Magnusson, Buizza, Lin, Chen and Emanuel, What is the predictability limit of midlatitude weather?, Journal of the Atmospheric Sciences (2019).
• Papers cited in the lesson: Luzzatto, Melbourne and Paccaut, The Lorenz attractor is mixing (2005); Kutz, Fu and Brunton, Multi-resolution dynamic mode decomposition (2016); Lusch, Kutz and Brunton, Deep learning for universal linear embeddings of nonlinear dynamics (2018); Champion, Lusch, Kutz and Brunton, Data-driven discovery of coordinates and governing equations (2019); Lam et al., GraphCast: Learning skillful medium-range global weather forecasting (2023); Huang and Turcotte, Are earthquakes an example of deterministic chaos?, Geophysical Research Letters 17 (1990).
• Libraries: PyDMD and PySINDy.