Lorenz, Lyapunov and why a weather model shows up in deep learning

00The only calculus you need

You can read this whole page knowing one idea. A differential equation is a rule that says how fast something is changing, given where it is right now. Write the thing as $y$ and its rate of change as $\dot{y}$ (the dot means "derivative with respect to time", the same $\frac{dy}{dt}$ from first-year calculus). The simplest possible rule is

which says "the speed at which $y$ grows is proportional to $y$ itself". Money in a savings account obeys this with $\lambda$ the interest rate. A cup of coffee cooling toward room temperature obeys it with $\lambda < 0$. The solution is the function $y(t)$ that satisfies the rule at every instant, and for this rule it is $y(t) = y_0 e^{\lambda t}$. You can check it: differentiate $y_0 e^{\lambda t}$ and you get $\lambda y_0 e^{\lambda t}$, which is $\lambda y$. Positive $\lambda$ means exponential growth, negative $\lambda$ means exponential decay. This is the entire reason $e^{\lambda t}$ appears throughout the slide's caption. Whenever the rate of change of an error is proportional to the error, the error grows or shrinks exponentially, and $\lambda$ is the rate.

How a computer solves one

Almost no differential equation has a closed-form solution like the one above. Instead you simulate. If you know $y$ right now and you know the rule for $\dot{y}$, then a tiny time $\Delta t$ later $y$ will be approximately

Take the step, recompute $\dot{y}$ at the new position, step again. This is the Euler method, and if the update rule looks familiar it should. Gradient descent is exactly Euler's method applied to the rule $\dot{\theta} = -\nabla L(\theta)$ with step size $\eta$ playing the role of $\Delta t$. The simulations on this page use RK4 (fourth-order Runge–Kutta), which is the Euler method with a smarter estimate of the slope over the step. It is more accurate for the same $\Delta t$ but the idea is identical. Every "trajectory" you see below was produced by stepping a rule forward thousands of times.

More than one variable

Real systems have several quantities changing at once, each with its own rule, and the rules refer to each other. Stack the quantities into a vector $\mathbf{s} = (s_1, \dots, s_n)$ and the rules into a vector-valued function $f$. Then $\dot{\mathbf{s}} = f(\mathbf{s})$ is a system of $n$ coupled equations and the solution is a path through $n$-dimensional space. That is an ordinary differential equation (ODE). "Ordinary" just means the unknowns depend on a single variable, time. A partial differential equation (PDE) is the case where the unknown is a whole field that varies over space as well as time, like the temperature at every point in a room. You will only meet PDEs once on this page, as the thing Lorenz started from and then simplified away.

One last piece of linear algebra you already have. If the rule is linear, $\dot{\mathbf{s}} = J\mathbf{s}$ for a matrix $J$, then along each eigenvector $\mathbf{v}$ of $J$ with eigenvalue $\mu$ the system behaves like the one-variable case, $\dot{(\cdot)} = \mu\,(\cdot)$, so the component along $\mathbf{v}$ grows like $e^{\mu t}$. A real negative $\mu$ decays, a real positive $\mu$ grows, and a complex $\mu = a + bi$ gives $e^{at}$ growth or decay multiplied by a rotation at frequency $b$, which is a spiral. That is all the eigenvalue talk in the next section means.

01Dynamical systems and stability

A dynamical system is any rule of the form above, $\dot{\mathbf{s}} = f(\mathbf{s})$ with $\mathbf{s} \in \mathbb{R}^n$. Some vocabulary. $\mathbf{s}$ is the state, everything you need to know about the system right now. $f$ is the vector field, an arrow attached to every point in state space saying which way the state moves next. The space of all possible states is phase space. A solution $\mathbf{s}(t)$ is a trajectory, a curve in phase space that follows the arrows. Sim 1 below draws the arrows so you can see all of this at once.

A fixed point is a state $\mathbf{s}^*$ where $f(\mathbf{s}^*) = 0$, so nothing moves. The important question is what happens near it. Zoom in close enough and any smooth function looks linear, so replace $f$ by its tangent plane at $\mathbf{s}^*$ (this is called linearizing). Writing $\boldsymbol{\delta} = \mathbf{s} - \mathbf{s}^*$ for the small displacement,

where $J$ is the Jacobian, the matrix of partial derivatives of $f$, the same object you backpropagate through. Now we are in the linear case from section 0. Along each eigenvector of $J$ with eigenvalue $\mu$, the displacement grows like $e^{\mathrm{Re}(\mu)\,t}$. So a fixed point is stable if every eigenvalue has negative real part (all displacements die out) and unstable if any eigenvalue has positive real part (some displacement grows). Complex eigenvalues add rotation, which is what makes spirals.

With $c > 0$ this is a damped spring. Both points spiral into the origin and their separation shrinks exponentially, so a small error in where you started does not matter for long. This is the world the first two mental models on the slide ("GD plus noise", "Brownian motion") live in. Small perturbations keep you close.

With $c < 0$ the origin is unstable and the separation grows exponentially. But notice this is boring too. Both points fly off to infinity along nearly the same path. Exponential growth of errors alone is not chaos. Chaos needs one more ingredient, which the next sections supply.

02The Lorenz equations and where they came from

In 1963 the meteorologist Edward Lorenz published "Deterministic Nonperiodic Flow". He was studying Rayleigh–Bénard convection, which is what happens when you heat a layer of fluid from below. Warm fluid rises, cools at the top, sinks, and the layer organizes into rolling cells. The full description is a PDE, because temperature and velocity are fields defined at every point in the fluid. Lorenz did something that should feel familiar from ML. He wrote the fields as a sum of spatial patterns (Fourier modes, which are just sine waves of different wavelengths), kept the three most important patterns, and threw the rest away. The coefficients of the three surviving patterns became three ordinary numbers that change in time, and their rules of motion are the Lorenz equations:

Loosely, $x$ is the speed of the convection roll, $y$ is the horizontal temperature difference between the rising and sinking fluid, and $z$ is how much the vertical temperature profile departs from a straight line. The parameters are $\sigma$ (the Prandtl number, a ratio of how fast momentum diffuses to how fast heat diffuses), $\rho$ (the Rayleigh number, how hard you are heating relative to the threshold where convection starts), and $\beta$ (a geometric aspect ratio of the cells). The standard values $\sigma = 10$, $\rho = 28$, $\beta = 8/3$ are the ones Lorenz used and the ones the slide uses.

You do not need the fluid story to use the equations. For the rest of this page, treat them as a three-dimensional vector field that is simple to write down and does something astonishing.

The fixed points

Setting all three derivatives to zero gives the origin, and, when $\rho > 1$, a symmetric pair

At the standard parameters, $\sqrt{\beta(\rho-1)} = \sqrt{72} \approx 8.49$ and $z = 27$. Those are the two horizontal lines on the slide. Physically $C_\pm$ are steady convection rolls turning clockwise or anticlockwise. Linearize at $C_\pm$ as in section 1. The Jacobian has one negative real eigenvalue (a direction that is pulled in) and a complex pair (a spiral). For $\rho$ below about $24.74$ the complex pair has negative real part, so the spiral winds inward and the rolls are stable. Above it the real part turns positive, the steady rolls are unstable, and yet nothing else is stable either. The trajectory has nowhere to settle.

03The attractor, one knob at a time

An attractor is the set that trajectories end up on after transients die away. A stable fixed point is an attractor. A stable periodic orbit (a closed loop) is an attractor. The Lorenz system at $\rho = 28$ has a strange attractor: trajectories are pulled onto a bounded, butterfly-shaped set, never leave it, never repeat, and never settle. Turn $\rho$ up from zero and watch the three regimes appear.

Three regimes to notice. For $\rho < 1$ the heating is too weak for convection and everything decays to the origin. For $1 < \rho < 24.74$ the trajectory spirals in to one of the two rolls $C_\pm$, with oscillations that shrink. The slide's first 17 time units are this same spiral, but with $\rho = 28$ the spiral is unstable so the oscillations grow instead of shrinking until the trajectory is thrown across to the other lobe. From then on it circulates around one lobe for a while, hops, circulates, hops, with no pattern to the counts.

Two things make this an attractor rather than an explosion. The divergence of the vector field, $\nabla\cdot f = \partial\dot x/\partial x + \partial\dot y/\partial y + \partial\dot z/\partial z$, measures whether the flow spreads a small blob of states out (positive) or compresses it (negative). For Lorenz it is the constant $-(\sigma + 1 + \beta) \approx -13.67$, so any blob of initial conditions shrinks in volume exponentially fast, which is why trajectories are pulled onto a thin set. And nonlinear terms $xz$ and $xy$ fold the stretched-out directions back on themselves so nothing escapes to infinity. Stretch in one direction, squeeze in the others, fold. That combination is the recipe for chaos. The result has zero volume but is not a surface either. Its fractal dimension is about $2.06$.

04The butterfly experiment

Now reproduce the slide. Run the system to $t = 20$, make a copy, nudge the copy's $x$ by $\delta_0$, and let both run. The two trajectories are both solutions of the exact same deterministic equations, with no noise anywhere. The only difference is $\delta_0$.

Read the lower plot. The separation climbs along a roughly straight line on the log axis, which means exponential growth, $\delta(t) \approx \delta_0 e^{\lambda t}$. It wobbles because the local stretching rate depends on where on the attractor you are (near the center of the butterfly the stretching is strongest). Then it saturates at about the size of the attractor, around 20 to 40 units, because two points on a bounded set cannot get further apart than that. After saturation the orange and blue curves are as unrelated as two trajectories started from random points.

The table below is the real punchline. On average, each factor of ten in precision buys only about $\ln 10 / \lambda \approx 2.5$ time units of agreement. The spacing is uneven from row to row because the stretching rate varies along the attractor, but the trend is unmistakable. That is the sense in which "infinitesimal at the input becomes macroscopic at the output". Linear improvements in forecast horizon cost exponential improvements in measurement.

Weather is the original application. Atmospheric models have a largest Lyapunov exponent that corresponds to errors roughly doubling every couple of days, and the initial state is only known to finite precision from sparse measurements. That is why the useful forecast horizon is around two weeks regardless of compute, and why forecasters run ensembles of slightly perturbed initial conditions instead of a single run. The ensemble spread is the orange-blue divergence, measured.

05Lyapunov exponents, properly

The largest Lyapunov exponent $\lambda_1$ is the average exponential growth rate of an infinitesimal perturbation, averaged along a trajectory for a long time:

The order of limits matters. You shrink the perturbation first so it never saturates, then take the time average. Numerically you do the equivalent trick of renormalizing the perturbation back to tiny size every few steps and summing the logs of the growth factors (the Benettin algorithm). For Lorenz at the standard parameters, $\lambda_1 \approx 0.906$ per unit time. The Lyapunov time $1/\lambda_1 \approx 1.1$ is the characteristic time for errors to grow by a factor of $e$. Errors double every $\ln 2 / \lambda_1 \approx 0.77$ time units.

Relation to the eigenvalues from section 1

At a fixed point, the Lyapunov exponents are exactly the real parts of the Jacobian eigenvalues. Everywhere else they generalize that idea. Instead of one Jacobian you have a product of many Jacobians along the trajectory, $M(t) = \prod_k e^{J(\mathbf{s}(t_k))\Delta t}$, and the Lyapunov exponents are the long-run growth rates of the singular values of $M(t)$, divided by $t$. The Oseledets theorem guarantees these limits exist for almost every starting point and do not depend on it. So $\lambda_1$ is a property of the attractor, not of where you began.

The full spectrum

An $n$-dimensional system has $n$ exponents, one per independent direction. For Lorenz they are approximately

Each sign tells a story. The positive one is the stretching that causes sensitivity to initial conditions. The zero one is always present for a continuous-time system on an attractor that is not a fixed point, because a perturbation along the direction of motion just shifts you forward or back in time and neither grows nor shrinks. The large negative one is the squeezing that pulls trajectories onto the thin attractor. Their sum, $-13.67$, equals the divergence $-(\sigma + 1 + \beta)$ from section 3, which is a nice consistency check: total volume contraction rate equals the sum of the per-direction rates.

Predictability horizon

If your initial error is $\delta_0$ and you can tolerate error $\Delta$, the time until you exceed it is

This is the formula behind the table in section 4. The logarithm is why heroic measurement effort gives disappointing returns.

06The double pendulum

The right half of the slide is a mechanical system with the same property. A pendulum with a second pendulum hanging from its tip has a four-dimensional state (two angles, two angular velocities). The equations come from Newton's laws written in terms of the two angles (the standard machinery for that is called Lagrangian mechanics, and you do not need it here). They are not pretty, but they are fully deterministic and energy-conserving. Released from a small angle it swings predictably. Released from a large angle it is chaotic, and two copies that differ by less than you could ever measure disagree completely within seconds.

At a 20° release the two copies stay together indefinitely. The motion is nearly a sum of two sine waves and the Lyapunov exponents are essentially zero (energy conservation forbids them from being negative, and small-angle motion is not chaotic). At 120° or more the upper bob has enough energy to flip over, the dynamics become chaotic, and the trails separate within a few seconds. Same equations, same physics, different region of phase space. Chaos is a property of a regime, not of a system as a whole.

07Chaos in discrete time

Everything above is continuous time. Training a network is an iterated map, $\mathbf{s}_{k+1} = F(\mathbf{s}_k)$, so it helps to see the same ideas with a step counter instead of a clock. For maps, the linearization is $\boldsymbol{\delta}_{k+1} = J_k\,\boldsymbol{\delta}_k$ with $J_k = \partial F / \partial \mathbf{s}$ at $\mathbf{s}_k$, and the largest Lyapunov exponent is

where $\mathbf{u}_k$ is the unit perturbation direction, renormalized each step. In one dimension $J_k$ is just the derivative $F'(s_k)$, so $\lambda_1$ is the average of $\ln|F'|$ along the orbit. If the map stretches on average ($|F'| > 1$ typically) you get chaos, if it squeezes you get convergence.

A continuous ODE needs at least three dimensions to be chaotic (in two dimensions trajectories cannot cross, which prevents the folding). A map can be chaotic in one dimension because a map can jump. The classic example is the logistic map $x_{k+1} = r\,x_k(1 - x_k)$, a toy population model.

Below $r = 3$ the orbit settles to a fixed point and $\lambda_1 < 0$. Between 3 and about 3.57 it settles into cycles of period 2, 4, 8 and so on (the period-doubling cascade) and $\lambda_1$ stays negative except exactly at the doubling points where it touches zero. Past about 3.57 the exponent becomes positive and the two orbits that started $10^{-8}$ apart stop agreeing after about 20 steps, for the same logarithmic reason as the table in section 4.

Gradient descent is an iterated map too

Take gradient descent on a one-dimensional quadratic $L(w) = \tfrac{h}{2}w^2$ with learning rate $\eta$. The update is $w_{k+1} = (1 - \eta h)\,w_k$. This is a linear map with derivative $1 - \eta h$, so its Lyapunov exponent is simply

For $0 < \eta h < 2$ this is negative and GD converges (monotonically for $\eta h < 1$, with oscillating sign for $1 < \eta h < 2$). At $\eta h = 2$ it is exactly zero. Beyond 2 it is positive and GD diverges. In many dimensions the same thing holds eigenvalue by eigenvalue with $h$ replaced by each eigenvalue of the Hessian, which is where the familiar stability condition $\eta < 2/h_{\max}$ comes from. The "edge of stability" phenomenon (Cohen et al., ICLR 2021, which the course covers in Lecture 6), where sharpness $h_{\max}$ climbs until it sits right at $2/\eta$, is training arranging itself so that the top direction has Lyapunov exponent about zero. Notice that a quadratic can be unstable but never chaotic, because the diverging iterates are unbounded. It takes the nonlinearity of a real loss to fold them back.

08Why this matters for deep learning

The slide puts "SGD is like a chaotic system" next to "SGD is like GD plus noise" and "SGD is like Brownian motion" and then asks, "But also, training deep nets?" The question matters because the three models make different predictions about what a small perturbation does. Under the first two, if you change the seed a little the final weights change a little. Under the chaotic model, any perturbation at all, no matter how small, gets amplified to the full run-to-run spread. Here is the evidence that the chaotic model is the right one for the weights.

One bit is enough

Summers and Dinneen (ICML 2021) isolated every source of nondeterminism in training a ResNet on CIFAR-10 one at a time: initialization seed, data order, data augmentation seed, and the nondeterministic floating-point reductions inside cuDNN. Each source alone produced about as much variability between trained models as all sources together. To explain that they ran the butterfly experiment on training: change a single weight at initialization by one unit in the last place of a float (about $6\times10^{-11}$) and keep everything else bit-identical. The resulting models were nearly as different from each other as models trained from entirely different random initializations. This is exactly the slide's orange-blue curve, with the network's weights as the state. And it explains a CIFAR slide earlier in the same lecture, which says the sources of randomness are not additive and "any single source leads to similar variance". Once the dynamics have a positive Lyapunov exponent, every perturbation saturates to the attractor's size, so the sources cannot add. The smallest one already gets you the full spread.

But the loss is not chaotic in the same way

This is the subtlety that makes the chaotic model useful rather than depressing. The weights diverge, yet the final test accuracy of the diverged runs is nearly the same (within a fraction of a percent on CIFAR). The attractor is enormous in weight space but the loss is nearly constant across it. The lecture's conclusion that "randomness is often not your biggest problem" is a statement about the loss, not the weights. Related evidence from Frankle, Dziugaite, Roy and Carbin (ICML 2020) found that after a short early phase of training, two copies of a network trained with different data orders end up in the same linearly connected basin, with no loss barrier on the straight line between them. Sensitivity to initial conditions in the weights is compatible with insensitivity in the quantities you care about. Keep both facts in your head at once.

Backprop through time is a product of Jacobians

A recurrent network unrolled for $T$ steps is an iterated map $\mathbf{h}_{t+1} = F(\mathbf{h}_t, x_t)$. The gradient of the loss at step $T$ with respect to $\mathbf{h}_0$ is

the same product of Jacobians whose growth rate defines the Lyapunov exponent. So exploding gradients (Pascanu, Mikolov and Bengio 2013) and a positive Lyapunov exponent of the forward dynamics are the same fact seen from two sides. Vanishing gradients are the negative-exponent case. Metz, Freeman, Schoenholz and Kachman ("Gradients are Not All You Need", 2021) push this further: when you differentiate through a chaotic system (a physics simulator, an unrolled optimizer in meta-learning, a long RNN), the gradient is formally correct but its variance explodes exponentially in the horizon, and a black-box estimator can beat the exact gradient. Their diagnostic is the spectrum of the Jacobians, which is this page's $\lambda_1$.

Signal propagation at initialization

A deep feed-forward net at initialization is also an iterated map, layer index playing the role of time. Poole, Lahiri, Raghu, Sohl-Dickstein and Ganguli (NeurIPS 2016) computed how two inputs' correlation evolves through random layers and found an ordered phase (inputs collapse together, negative exponent) and a chaotic phase (nearby inputs are pulled apart, positive exponent) depending on the weight variance. Schoenholz et al. (ICLR 2017) then showed that trainability of very deep nets requires initializing near the boundary, the edge of chaos, where the exponent is near zero and signals neither die nor explode across depth. Residual connections and normalization are ways of keeping a deep net near that boundary without tuning.

09Blind spots and glossary

Ideas that sit right next to the slide and are easy to miss when you first meet it.

Chaos, stochasticity and instability are three different things

Stochastic means randomness is injected. Unstable means a linearized perturbation grows (an eigenvalue beyond the unit circle). Chaotic means deterministic, bounded, and exponentially sensitive. Full-batch GD showing spikes (the slide right after "Chaotic systems" in the same lecture, titled "stochastic behavior in deterministic training") is instability at the edge of stability, not chaos. The slide is careful to say chaotic systems look random.

Lyapunov exponent is not an eigenvalue

Eigenvalues of a single Jacobian describe one point. The Lyapunov exponent is a long-time average over a whole trajectory. The two coincide only at a fixed point. A system can have every local Jacobian contracting in some direction yet still be chaotic because the directions rotate.

Finite-time Lyapunov exponents

Over a finite window the growth rate fluctuates, which is why the log-separation curve in Sim 3 wobbles. People report finite-time exponents to describe where on the attractor sensitivity is worst. The asymptotic exponent is a single number for the whole attractor.

Numerical trajectories are wrong but statistically right

Any computer trajectory of a chaotic system departs from the true one after a few Lyapunov times because of rounding. The shadowing lemma says that for a wide class of systems there is nonetheless a true trajectory that stays close to the computed one. So you cannot trust the specific path, but you can trust averages over the attractor. The same reasoning is why we trust loss curves from a single nondeterministic GPU run as statistically representative.

Ergodicity

Time averages along one trajectory equal space averages over the attractor. This is what lets you estimate $\lambda_1$ from a single long run and what justifies the Oseledets result that $\lambda_1$ does not depend on where you start.

Dissipative vs conservative

Lorenz is dissipative (volume shrinks, attractor exists). The double pendulum is conservative (Hamiltonian, volume preserved by Liouville's theorem, no attractor). Conservative systems can still be chaotic, with Lyapunov exponents summing to zero and coming in $\pm$ pairs. Training with SGD and weight decay is dissipative. Pure Hamiltonian dynamics (as in HMC samplers) are conservative.

Strange attractor and fractal dimension

The Lorenz attractor has Kaplan–Yorke dimension $\approx 2.06$, computed from the exponents as $2 + \lambda_1/|\lambda_3|$. It is more than a surface and less than a volume. The analogous question for training, what is the geometry of the set of weights that SGD visits, is open.

Chaos was proved, not just simulated

Lorenz observed chaos numerically in 1963. A rigorous computer-assisted proof that the attractor exists and is strange came only in 2002, from Warwick Tucker. It is a reminder that "the plot looks chaotic" and "the system is chaotic" are different claims.

Sources

• E. N. Lorenz, "Deterministic Nonperiodic Flow", J. Atmos. Sci. 20, 1963.
• S. Strogatz, Nonlinear Dynamics and Chaos, ch. 9 (the Lorenz chapter) for every formula in sections 2 to 5.
• C. Summers and M. J. Dinneen, "Nondeterminism and Instability in Neural Network Optimization", ICML 2021. arXiv:2103.04514
• L. Metz, C. D. Freeman, S. S. Schoenholz, T. Kachman, "Gradients are Not All You Need", 2021. arXiv:2111.05803
• J. Frankle, G. K. Dziugaite, D. M. Roy, M. Carbin, "Linear Mode Connectivity and the Lottery Ticket Hypothesis", ICML 2020. arXiv:1912.05671
• R. Pascanu, T. Mikolov, Y. Bengio, "On the difficulty of training recurrent neural networks", ICML 2013. arXiv:1211.5063
• B. Poole et al., "Exponential expressivity in deep neural networks through transient chaos", NeurIPS 2016 (arXiv:1606.05340). S. Schoenholz et al., "Deep Information Propagation", ICLR 2017 (arXiv:1611.01232).
• J. Cohen, S. Kaur, Y. Li, J. Z. Kolter, A. Talwalkar, "Gradient Descent on Neural Networks Typically Occurs at the Edge of Stability", ICLR 2021. arXiv:2103.00065
• T. Hashimoto, Stanford CS 312 (Deep Learning Alchemy), Lecture 2 slides, "Mental models of deep learning randomness" and "Chaotic systems". lecture_02.pdf