Interactive Stochastic Differential Equations

Definition. \(S_n=\sum_{i=1}^{n}\xi_i\), where the \(\xi_i\) are independent and take \(\pm 1\) with probability 1/2 each. \(S_n\) is the position after \(n\) steps, and the essential point is that each step's decision is independent of the past.

Mean 0, variance n. By linearity of expectation \(\mathbb{E}[S_n]=0\). Because the steps are independent the variances add, giving \(\mathrm{Var}(S_n)=n\), so the standard deviation is \(\sqrt{n}\). The distance is proportional not to time but to the square root of time — walk 100 steps and you are on average only 10 steps from the start. This \(\sqrt{t}\) law is the signature of diffusion.

A Gaussian appears. By the central limit theorem \(S_n/\sqrt{n}\) converges to a standard normal as \(n\) grows. Even though the step distribution was the crude discrete \(\pm 1\), adding many of them gives a smooth bell shape. Increasing the number of paths in the demo shows the histogram of endpoints filling into a Gaussian.

It meets the diffusion equation. Running a random walk with lattice spacing \(\Delta x\) and time step \(\Delta t\) makes the probability distribution satisfy a difference equation in \(p(t,x)\). Taking the limit while holding \(\Delta x^2/\Delta t\to 2\alpha\) gives exactly the diffusion equation \(p_t=\alpha p_{xx}\). That the fundamental solution of diffusion in the previous section was a Gaussian of variance \(2\alpha t\) corresponds precisely to \(\mathrm{Var}=n\) here — revealing that microscopic random motion and the macroscopic diffusion PDE are two faces of the same thing.

In machine learning. The parameter trajectory of SGD is modelled as a random walk driven by minibatch noise superimposed on the drift of gradient descent. This is why parameter fluctuations late in training have scale \(\sqrt{\eta}\), and it is the starting point when discussing implicit regularization effects. The forward noise injection of a diffusion model is also essentially this random walk.

\[\begin{aligned} S_n &= \textstyle\sum_{i=1}^n \xi_i, \quad \xi_i = \pm 1 \\ \mathbb{E}[S_n] = 0, \quad &\mathrm{Var}(S_n) = n, \quad \sigma = \sqrt{n} \\ S_n/\sqrt{n} &\;\xrightarrow{d}\; \mathcal{N}(0,1) \end{aligned}\]

Increase the number of paths. Individual paths differ wildly, yet the whole is confined within a fan of width √n (the grey curves are ±σ = ±√n). The histogram on the right filling into a Gaussian is the central limit theorem, and its width growing as the square root of time is diffusion.

Definition. Noise whose value at each instant follows a Gaussian distribution and is independent across instants. In discrete time it appears as \(\xi_t\sim\mathcal{N}(0,\sigma^2)\), and in continuous time as the increment of Brownian motion \(dW\sim\mathcal{N}(0,dt)\).

Stability: adding Gaussians gives a Gaussian. The sum of two independent Gaussians is again Gaussian, and the variances simply add: \(\mathcal{N}(0,\sigma_1^2)+\mathcal{N}(0,\sigma_2^2)=\mathcal{N}(0,\sigma_1^2+\sigma_2^2)\). Thanks to this property, many stages of noise injection can be skipped in one go. That a diffusion model samples an arbitrary \(x_t\) from \(x_0\) in a single step without passing through the intermediate stages — \(x_t=\sqrt{\bar\alpha_t}x_0+\sqrt{1-\bar\alpha_t}\epsilon\) — rests exactly on this property. It is the fundamental reason training is practical.

Maximum entropy. As seen in the information theory article, the distribution maximizing entropy for a given variance is the Gaussian. That is, it is noise that "injects no structure other than a magnitude," making it the least biased tool for destroying data.

The central limit theorem guarantees it. Even if the real noise is not Gaussian, accumulating many small independent fluctuations brings the result close to a Gaussian. The \(\pm 1\) steps producing a Gaussian in the previous section is an instance. So the Gaussian assumption is both a convenience and, in the limit, the correct choice.

The score becomes simple. The gradient of a Gaussian's log-density is linear: \(\nabla_x\log p=-(x-\mu)/\sigma^2\). This simplicity is used decisively later — it is why, when a diffusion model learns the score, the target essentially reduces to "guess the noise that was added."

\[\begin{aligned} \mathcal{N}(0,\sigma_1^2) + \mathcal{N}(0,\sigma_2^2) &= \mathcal{N}(0,\sigma_1^2+\sigma_2^2) \\ x_t &= \sqrt{\bar\alpha_t}\,x_0 + \sqrt{1-\bar\alpha_t}\,\epsilon \\ \nabla_x \log \mathcal{N}(x;\mu,\sigma^2) &= -\frac{x-\mu}{\sigma^2} \end{aligned}\]

Increase the number of stages and add the noise in several instalments. However many instalments you split it into, the final distribution is identical — the grey (added at once) and the orange (added in instalments) overlap. Without this property, training a diffusion model would have required walking hundreds of stages in sequence every time.

Definition. The Wiener process \(W_t\) is defined by four conditions. ① \(W_0=0\). ② Increments over non-overlapping intervals are independent. ③ \(W_t-W_s\sim\mathcal{N}(0,\,t-s)\). ④ The paths are continuous. Note in particular that the variance of an increment is the time interval itself.

Continuous yet differentiable nowhere. Since the size of an increment is of order \(\sqrt{\Delta t}\), the difference quotient is \(\Delta W/\Delta t\sim 1/\sqrt{\Delta t}\), which diverges as \(\Delta t\to 0\). So \(\dot{W}\) does not exist. This is why a stochastic differential equation cannot be written in the form \(\frac{dX}{dt}=\ldots\) and is written in the differential form \(dX=\ldots\,dW\). Zooming in on the demo shows that however closely you look it never becomes smooth, repeating the same roughness — self-similarity.

The scaling property. \(W_{ct}\) and \(\sqrt{c}\,W_t\) have the same distribution. It means that stretching time by \(c\) spreads space by \(\sqrt{c}\) — the \(\sqrt{n}\) law of the random walk carried into the continuous version.

It is the same thing as the diffusion equation. The probability density of \(W_t\) is exactly \(p(t,x)=\frac{1}{\sqrt{2\pi t}}e^{-x^2/2t}\), which is the fundamental solution of the diffusion equation with \(\alpha=1/2\). That is, Brownian motion is the "particle version" of the diffusion equation — the most important instance of "From Particle to Density Dynamics" three sections back, generalized later into Fokker–Planck.

In machine learning. The forward process of a diffusion model is Brownian motion (with a drift added). What the noise schedule \(\sqrt{1-\bar\alpha_t}\) does is set the time scale of this process, and "how much noise to add" is the same question as "how long to diffuse for."

\[\begin{aligned} W_0 &= 0, \quad W_t - W_s \sim \mathcal{N}(0, t-s) \\ \mathbb{E}[W_t] = 0, \quad &\mathrm{Var}(W_t) = t, \quad \mathbb{E}[W_sW_t] = \min(s,t) \\ W_{ct} &\overset{d}{=} \sqrt{c}\,W_t \end{aligned}\]

Raise the zoom factor and peer into part of a path. However far you zoom it never becomes smooth — you are looking at non-differentiability with your own eyes. The grey curves are the ±√t range, and the paths spread out within it.

Definition. A stochastic process is a collection of random variables \(\{X_t\}_{t\in T}\) indexed by a set \(T\). If \(T\) is discrete it is a time series; if continuous, a continuous-time process. Fixing a single \(\omega\) gives a function of time \(t\mapsto X_t(\omega)\), called a sample path.

It can be cut two ways. Fixing a time \(t\) makes \(X_t\) an ordinary random variable whose distribution we may ask about (the marginal). Fixing \(\omega\) instead gives a single path. The vertical cross-section and the horizontal curve of the demo are these two respectively. Understanding a stochastic process is moving freely between these two views — and indeed training a diffusion model deals with the marginals while sampling deals with paths.

Finite-dimensional distributions determine the process. Knowing the joint distribution at arbitrary times \(t_1<\cdots<t_k\) determines the process (by Kolmogorov's extension theorem). That Brownian motion could be defined by only four conditions is because those conditions pin down all its finite-dimensional distributions.

Filtration. \(\mathcal{F}_t\) means "the information observed up to time \(t\)," and it grows as time passes. This notion is needed to state rigorously that a decision does not peek at the future. It explains why the Itô integral of the coming sections uses the left endpoint of each interval — because the decision must be made now, without knowing the future.

In machine learning. The \(\{x_t\}\) of a diffusion model is a stochastic process, and training is the problem of matching the score of the marginal \(p_t(x)\) at each time. Reinforcement learning state trajectories, time series forecasting, and Gaussian process regression all use this language.

\[\begin{aligned} \{X_t\}_{t\in T}, \quad &X_t : \Omega \to \mathbb{R} \\ t \text{ fixed} &\;\Rightarrow\; \text{random variable } X_t \sim p_t \\ \omega \text{ fixed} &\;\Rightarrow\; \text{sample path } t\mapsto X_t(\omega) \\ \mathcal{F}_s &\subseteq \mathcal{F}_t \quad (s \le t) \end{aligned}\]

Move the time slider and the vertical cross-section at that instant is drawn as a distribution on the right. One curve running horizontally is a sample path, and a vertical cut is the random variable at that instant — two ways of seeing the same object.

Definition. A process \(\{X_t\}\) is Markov when \(p(X_t\mid X_s, X_u\;(u<s))=p(X_t\mid X_s)\) for all \(s<t\). Knowing the present, the past adds no information about the future.

Why this is decisive. Without the Markov property, drawing the next state would require carrying the entire history. With it, only the current state is needed, so memory is constant and simulation proceeds one step at a time. This is why both the Euler method and Euler–Maruyama can take the form "make the next value from the current one."

Transition probabilities and Chapman–Kolmogorov. Writing \(p(t,y\mid s,x)\) for the transition density, we have \(p(t,y\mid s,x)=\int p(t,y\mid u,z)\,p(u,z\mid s,x)\,dz\) for \(s<u<t\). It says "marginalize over every path passing through an intermediate time," and differentiating this relation as \(u\to t\) gives the Kolmogorov equations later on.

Brownian motion is Markov. Because its increments are independent, knowing only the current position determines the future distribution. How it got here does not matter. Combining this section with the previous one explains why a process described by an SDE is Markov — because its increments are built only from the current state and independent noise.

In machine learning. The forward process of DDPM is explicitly defined as a Markov chain: \(q(x_t\mid x_{t-1})\). The reverse process is also assumed Markov and \(p_\theta(x_{t-1}\mid x_t)\) is learned. The MDPs of reinforcement learning, state space models for language (Mamba and others), and MCMC samplers all stand on this property. Interestingly, an autoregressive language model is not Markov — it conditions on the entire history, which is why context length becomes a cost.

\[\begin{aligned} p(X_t \mid X_s, \text{past}) &= p(X_t \mid X_s) \\ p(t,y\mid s,x) &= \int p(t,y\mid u,z)\,p(u,z\mid s,x)\,dz \\ q(x_{1:T}\mid x_0) &= \textstyle\prod_t q(x_t \mid x_{t-1}) \end{aligned}\]

Switch to "with memory" and the next step is pulled toward an average of several past steps — the path becomes smoother, but simulating it requires carrying the whole history. In Markov mode the next point is determined by the current point alone, which is the premise of every step-by-step sampler.

Definition. Define the change over a short time \(\Delta t\) by conditional moments. The drift \(f(x,t)=\lim_{\Delta t\to 0}\frac{1}{\Delta t}\mathbb{E}[\Delta X\mid X=x]\) is the average rate of change, and the diffusion coefficient \(g^2(x,t)=\lim_{\Delta t\to 0}\frac{1}{\Delta t}\mathbb{E}[(\Delta X)^2\mid X=x]\) is the rate of growth of the variance.

Why is one \(\Delta t\) and the other \(\sqrt{\Delta t}\)? The displacement produced by drift is proportional to \(\Delta t\), while the random displacement is proportional to \(\sqrt{\Delta t}\) (the \(\sqrt{n}\) law of the random walk). So when \(\Delta t\) is very small the random term is overwhelmingly larger. This is why the zigzag dominates on short views while direction emerges on long ones, and this difference of scale is the source of every peculiarity of Itô calculus in the next section.

Their competition creates a stationary distribution. If drift pulls toward the origin while diffusion spreads outward, the distribution settles where the two balance. The Ornstein–Uhlenbeck process \(dX=-\theta X\,dt+\sigma\,dW\) is the example, with stationary distribution \(\mathcal{N}(0,\sigma^2/2\theta)\). The forward process of a diffusion model has exactly this form, which is why running it long enough reaches a standard Gaussian whatever the data was.

Seen in the language of densities. Drift carries the density bodily (a transport term) while diffusion blurs it (a diffusion term). The two terms of the Fokker–Planck equation \(\partial_t p=-\nabla\cdot(fp)+\frac{1}{2}\nabla^2(g^2p)\) later on are exactly these two — think of it as the continuity equation and the diffusion equation combined into a single equation.

In machine learning. The VP-SDE (corresponding to DDPM) takes \(f=-\frac{1}{2}\beta(t)x\) and \(g=\sqrt{\beta(t)}\), so the drift contracts toward the origin while diffusion fills in noise. The VE-SDE (SMLD) takes \(f=0\), growing the variance with no drift at all. Designing a noise schedule is designing these two coefficients.

\[\begin{aligned} f(x,t) &= \lim_{\Delta t\to 0}\tfrac{1}{\Delta t}\mathbb{E}[\Delta X\mid x] &&\text{(drift)} \\ g^2(x,t) &= \lim_{\Delta t\to 0}\tfrac{1}{\Delta t}\mathbb{E}[(\Delta X)^2\mid x] &&\text{(diffusion)} \\ \Delta X &\approx f\,\Delta t + g\sqrt{\Delta t}\,\xi, \quad \xi\sim\mathcal{N}(0,1) \end{aligned}\]

Move the drift and diffusion strengths separately. With drift alone you get a deterministic curve, with diffusion alone it spreads without direction, and with both it spreads while being pulled. On the OU preset the two forces balance and the distribution stops at a fixed width — the principle by which a diffusion model reaches a Gaussian.

Definition. A stochastic differential equation is written in the form \(dX_t=f(X_t,t)\,dt+g(X_t,t)\,dW_t\), where \(f\) is the drift, \(g\) the diffusion coefficient and \(W_t\) Brownian motion. It is not written as the quotient \(\frac{dX}{dt}\) because \(W_t\) cannot be differentiated.

It is really an integral equation. The differential form is only shorthand; the real meaning is the integral equation \(X_t=X_0+\int_0^t f\,ds+\int_0^t g\,dW_s\). The first integral is an ordinary Riemann integral, but the second is an object that must be newly defined — which is the Itô integral of the next section.

The solution is a distribution, not a path. An ODE matched one trajectory to one initial value. An SDE gives a different path every time from the same initial value. So "solution" usually means one of two things — a strong solution with respect to a given Brownian path, or a weak solution agreeing only in distribution. What matters in practice is usually the marginal distribution at each time.

An ODE is a special case. With \(g=0\) it reduces to an ordinary ODE. So an SDE is an extension of an ODE, and the numerical methods and stability concepts learned earlier carry over (with modifications).

In machine learning. The entire forward process of a diffusion model is a single SDE. The score-based viewpoint of Song et al. unified DDPM and SMLD as two SDEs, the VP-SDE and the VE-SDE, and at that moment the discrete choice of "number of noise steps" became the continuous object of a "time interval," opening the way to any numerical method. Improvements in sampling speed poured out from here.

\[\begin{aligned} dX_t &= f(X_t,t)\,dt + g(X_t,t)\,dW_t \\ X_t &= X_0 + \int_0^t f\,ds + \int_0^t g\,dW_s \\ g \equiv 0 &\;\Rightarrow\; \text{ODE} \end{aligned}\]

Draw several times from the same initial value. The ODE (grey) is always the same curve, while the SDE (coloured) gives a different path every time. Turn the diffusion strength down to 0 and every path collapses onto the grey curve — meaning an SDE contains an ODE.

Definition. The Itô integral \(\int_0^T H_s\,dW_s\) is defined as a limit of Riemann sums, but always using the value at the left endpoint of each interval: \(\sum_i H_{t_i}(W_{t_{i+1}}-W_{t_i})\). This choice is what defines the Itô integral.

Why the left endpoint? There are two reasons. First, so as not to peek at the future — on the interval \([t_i,t_{i+1}]\) we must decide \(H\) without knowing what the increment will be. In finance, "decide how many shares to buy now without knowing whether the price will rise" is exactly this situation. Second, as a consequence \(H_{t_i}\) becomes independent of the increment, so \(\mathbb{E}\left[\int H\,dW\right]=0\) holds (the martingale property). The computation simplifies dramatically.

The choice changes the answer. In ordinary integration the limit is the same whether one uses the left point, the midpoint or the right point. But a Brownian path is so rough that each choice converges to a different value. Using the midpoint gives the Stratonovich integral, which satisfies the ordinary chain rule but loses the martingale property. In the demo you can see the three rules applied to the same path and the values separating.

A concrete example. Computing \(\int_0^T W\,dW\), ordinary calculus intuition suggests \(\frac{1}{2}W_T^2\), but the Itô integral gives \(\frac{1}{2}W_T^2-\frac{1}{2}T\). An extra \(-T/2\) appears, and the identity of that term is the quadratic variation of the next section.

In machine learning. Diffusion models and SDE-based generative models conventionally use the Itô convention. The Euler–Maruyama discretization taking the form "compute the coefficients at the current state and step to the next" corresponds exactly to Itô's left-endpoint rule, which is why implementation is natural.

\[\begin{aligned} \int_0^T H\,dW &= \lim \sum_i H_{t_i}\,(W_{t_{i+1}}-W_{t_i}) \\ \mathbb{E}\Big[\int_0^T H\,dW\Big] &= 0 \\ \int_0^T W\,dW &= \tfrac{1}{2}W_T^2 - \tfrac{1}{2}T \end{aligned}\]

Increase the number of subdivisions. For an ordinary function the three rules gather to the same value, but on a Brownian path they remain separated to the end. Check in the readout that the difference between Itô (left) and Stratonovich (midpoint) converges to exactly T/2.

Definition. For a partition \(0=t_0<\cdots<t_n=T\), \([X]_T=\lim_{\|\Delta\|\to 0}\sum_i (X_{t_{i+1}}-X_{t_i})^2\). For Brownian motion, remarkably, \([W]_T=T\) — and this value is deterministic, not random.

Compare with a smooth function. If \(f\) is differentiable then \(\Delta f\approx f'\Delta t\), so \(\sum(\Delta f)^2\approx\sum (f')^2\Delta t^2\to 0\). It vanishes as the partition is refined. But Brownian motion has \(\Delta W\sim\sqrt{\Delta t}\), so \((\Delta W)^2\sim\Delta t\), and adding \(n\) of them leaves \(n\cdot\frac{T}{n}=T\). In the demo you can see the value refusing to go to 0 as the partition is refined, clinging to \(T\) instead.

The symbolic rule \(dW^2=dt\). In practice this fact is written as the multiplication rules \((dW)^2=dt\), \(dW\,dt=0\), \((dt)^2=0\). That is, \(dW\) is an object of size \(\sqrt{dt}\), so its square survives as a first-order term. Second-order terms that could ordinarily be discarded in a Taylor expansion can no longer be discarded here, and the consequence is Itô's lemma in the next section.

Volatility can be measured. \([X]_T\) is a quantity computable from a single observed path. So in finance it is used to estimate realized volatility from high-frequency data, and it becomes the route by which the diffusion coefficient of a stochastic model is identified from data.

In machine learning. This property governs the order of accuracy of SDE numerical methods. The reason Euler–Maruyama converges more slowly than deterministic Euler (strong order 1/2) is the term of size \(\sqrt{\Delta t}\), and it connects to the phenomenon that stochastic samplers degrade faster than deterministic ones when the number of sampling steps in a diffusion model is reduced.

\[\begin{aligned} [X]_T &= \lim \sum_i (X_{t_{i+1}} - X_{t_i})^2 \\ [W]_T &= T \quad \text{(deterministic!)} \\ (dW)^2 = dt, \quad &dW\,dt = 0, \quad (dt)^2 = 0 \end{aligned}\]

Increase the number of subdivisions. The quadratic variation of the smooth curve (grey) vanishes to 0, while that of the Brownian path (blue) converges to T. The first variation (the sum of absolute values) does the opposite and diverges to infinity for the Brownian path — these two facts together make the strange property "non-differentiable yet of finite quadratic variation."

Definition. If \(dX=f\,dt+g\,dW\) and \(Y=\varphi(X,t)\), then \[dY=\Big(\frac{\partial\varphi}{\partial t}+f\frac{\partial\varphi}{\partial x}+\frac{1}{2}g^2\frac{\partial^2\varphi}{\partial x^2}\Big)dt+g\frac{\partial\varphi}{\partial x}\,dW.\] An extra term \(\frac{1}{2}g^2\varphi_{xx}\) is attached to the ordinary chain rule.

Where does it come from? Expanding to second order in a Taylor series gives \(d\varphi=\varphi_t dt+\varphi_x dX+\frac{1}{2}\varphi_{xx}(dX)^2+\cdots\). Ordinarily \((dX)^2\) would be of order \((dt)^2\) and discarded. But since \((dW)^2=dt\) from the previous section, \((dX)^2=g^2(dW)^2=g^2dt\) becomes a first-order term and survives. Itô's lemma is, in effect, "the chain rule with a second-order Taylor term that can no longer be thrown away."

The Jensen effect. The sign of the extra term is set by \(\varphi_{xx}\). If \(\varphi\) is convex (\(\varphi_{xx}>0\)) then randomness lifts the mean of \(Y\). It is the same principle as \(\mathbb{E}[X^2]>0\) when \(X\) is mean-zero noise, and in the demo you can watch the mean of \(Y=X^2\) pushed upward in proportion to the diffusion strength.

The canonical example: geometric Brownian motion. Applying \(\varphi=\log S\) to \(dS=\mu S\,dt+\sigma S\,dW\) gives \(d\log S=(\mu-\frac{\sigma^2}{2})dt+\sigma\,dW\). Ordinary computation would have produced only \(\mu\,dt\), yet an extra \(-\sigma^2/2\) appears. This term, saying that greater volatility lowers the mean log return, is the heart of the Black–Scholes formula.

In machine learning. That the mean and variance of \(x_t\) come out in closed form for the forward SDE of a diffusion model, and the derivation of the reverse SDE later on, both rest on Itô's lemma. It is also used when computing the time evolution of the log-density \(\log p_t(x)\), and the result of that is the Fokker–Planck equation.

\[\begin{aligned} dY &= \Big(\varphi_t + f\varphi_x + \tfrac{1}{2}g^2\varphi_{xx}\Big)dt + g\,\varphi_x\,dW \\ d(W^2) &= dt + 2W\,dW \\ d\log S &= \big(\mu - \tfrac{\sigma^2}{2}\big)dt + \sigma\,dW \end{aligned}\]

Change the function φ and raise the diffusion strength. For a convex function (x²) the mean of Y is pushed upward, and for a concave one (log) downward. The grey dashed line is the prediction from applying only the ordinary chain rule, and the gap to the actual mean (orange) is exactly the ½g²φ″ term.

Definition. \(X_{n+1}=X_n+f(X_n,t_n)\Delta t+g(X_n,t_n)\sqrt{\Delta t}\,\xi_n\), where \(\xi_n\sim\mathcal{N}(0,1)\) is drawn afresh each step. The deterministic term is the Euler method unchanged; only the size of the noise term, \(\sqrt{\Delta t}\), is different.

Why \(\sqrt{\Delta t}\)? Because the standard deviation of a Brownian increment is \(\sqrt{\Delta t}\) (\(W_{t+\Delta t}-W_t\sim\mathcal{N}(0,\Delta t)\)). Using \(\Delta t\) here would make the noise far too small and the scheme would converge, as \(\Delta t\to 0\), to an ODE rather than an SDE. Omitting this square root is the most common bug in implementations.

Two kinds of convergence. Strong convergence measures how close the paths themselves are for the same Brownian path, and Euler–Maruyama has order \(1/2\). Weak convergence measures how close the distributions (expectations) are, and there the order is \(1\). What matters for generative models is weak convergence — there is no need to reproduce a particular path, only to draw samples from the correct distribution. In the graph you can see the two errors descending with different slopes.

Why higher-order methods are hard. For deterministic ODEs we obtained fourth order easily with RK4, but for SDEs one must handle iterated integrals of \(dW\) (Lévy areas and so on), making higher-order methods far more complicated and expensive. The Milstein method raises the strong order to 1 but requires the derivative of \(g\).

In machine learning. The stochastic samplers of diffusion models (ancestral sampling, SDE samplers) are essentially Euler–Maruyama. Raising \(\Delta t\) to reduce the number of steps grows the discretization error of the noise term and degrades sample quality, and this difference in order is precisely why deterministic ODE samplers are favourable at low step counts — the practical motivation for the Probability Flow ODE of the final section.

\[\begin{aligned} X_{n+1} &= X_n + f\,\Delta t + g\sqrt{\Delta t}\,\xi_n, \quad \xi_n\sim\mathcal{N}(0,1) \\ \text{strong convergence} &: \; \mathbb{E}|X_T - \hat{X}_T| = O(\Delta t^{1/2}) \\ \text{weak convergence} &: \; |\mathbb{E}\varphi(X_T) - \mathbb{E}\varphi(\hat{X}_T)| = O(\Delta t) \end{aligned}\]

Raising Δt makes the paths rougher, and using Δt instead of √Δt (the toggle) makes the noise vanish so the trajectory collapses onto a deterministic curve — the discretization has turned the SDE into an ODE. Check the difference between strong and weak order of convergence in the graph below.

The reverse SDE (Anderson, 1982). Reversing the forward \(dX=f\,dt+g\,dW\) in time gives \[dX=\big[f(X,t)-g^2(t)\,\nabla_x\log p_t(X)\big]dt+g(t)\,d\bar{W},\] where \(\bar{W}\) is reverse-time Brownian motion. Only one term, \(-g^2\nabla_x\log p_t\), is added to the forward coefficients.

What is the score? \(s(x,t)=\nabla_x\log p_t(x)\) is called the score, and it is "the vector field pointing where the density increases." Being a gradient, it is orthogonal to the contours and points toward the peaks. If diffusion scattered the data, the reverse process walks back along these arrows to where the data used to be.

Why extra information is needed. We said diffusion is irreversible. Looking at an individual path alone, that is correct — it cannot be undone. But knowing the information of the whole distribution allows the inference "a particle arriving at this position most likely came from over there." The score is exactly that information, and the point of the theorem is that reversing the forward process requires knowing the marginal distribution at every time.

The normalizing constant disappears. The decisive advantage of the score is that \(\nabla_x\log p=\nabla_x\log\tilde{p}-\nabla_x\log Z=\nabla_x\log\tilde{p}\). The intractable normalizing constant \(Z\) vanishes under differentiation. Learning a density directly means wrestling with \(Z\), while the score needs no such thing — the fundamental reason score-based modelling works.

In machine learning: this is a diffusion model. Approximating the score by a neural network \(s_\theta(x,t)\) lets us solve the reverse SDE and generate data from noise. And because the forward process is Gaussian, the score has a linear form, so the problem of matching the score reduces to "a regression problem of guessing the noise \(\epsilon\) that was added" — which is DDPM's loss \(\|\epsilon-\epsilon_\theta(x_t,t)\|^2\), meshing exactly with the fact that MSE is a Gaussian NLL.

\[\begin{aligned} \text{forward: } dX &= f\,dt + g\,dW \\ \text{reverse: } dX &= \big[f - g^2\,\nabla_x\log p_t(X)\big]dt + g\,d\bar{W} \\ s_\theta(x,t) &\approx \nabla_x\log p_t(x) = -\epsilon/\sigma_t \end{aligned}\]

Play the forward process and the two-peaked distribution is blurred into a Gaussian. Then play the reverse and the two peaks are revived by following the score (the arrows). Press "Score off" to run the reversal without that term — the peaks are not restored and it stays blurred. That term is everything.

The marginal. For a solution \(X_t\) of an SDE, the probability density \(p_t(x)\) at each time is called the marginal distribution of the process. Individual paths are random, but \(p_t\) evolves completely deterministically — starting from the same initial distribution always gives the same \(p_t\).

Where randomness disappears. This is the single most important shift of perspective in this article. In the language of paths a different answer comes out every time, but translated into the language of densities the answer is one. It is because taking an average over infinitely many paths washes the chance away, and it is "From Particle to Density Dynamics" repeated in a random version.

Two forces. There are only two factors changing the density. Drift carries it bodily (the continuity equation's term), and diffusion blurs it (the diffusion equation's term). The Fokker–Planck equation of the next section is nothing more than these two terms written in one equation.

Why this viewpoint is needed. The goal of a generative model is "to draw samples from the data distribution," not "to reproduce a particular path." So what we actually want to control is \(p_t\). And the fact that several processes can produce the same \(p_t\) — random ones and deterministic ones alike — becomes the heart of the final section.

In machine learning. The training objective of a diffusion model is to match the score of \(p_t\) at each time, and training is possible because \(p_t\) is given in closed form (a Gaussian convolution). Flow Matching designs the desired \(p_t\) path first and then finds the velocity field realizing it, making the most direct use of this section's viewpoint.

\[\begin{aligned} p_t(x) &= \text{the density of } t \text{ at time } X_t \text{} \\ \text{paths: random} &\quad\longleftrightarrow\quad p_t: \text{deterministic} \\ \partial_t p_t &= \underbrace{-\nabla\cdot(f p_t)}_{\text{transport}} + \underbrace{\tfrac{1}{2}\nabla^2(g^2 p_t)}_{\text{diffusion}} \end{aligned}\]

Increase the number of paths. Individual paths differ every time, yet the histogram (blue) always converges to the same curve (orange). Pressing "Resample" repeatedly does not shake the density — you are watching randomness disappear at the level of the density.

Definition. The marginal density \(p_t(x)\) of the SDE \(dX=f(X,t)\,dt+g(t)\,dW\) satisfies the PDE \[\frac{\partial p_t}{\partial t}=-\nabla\cdot\big(f(x,t)\,p_t\big)+\frac{1}{2}g^2(t)\,\nabla^2 p_t.\] This is called the Fokker–Planck equation (or the Kolmogorov forward equation).

The identity of the two terms. The first term is the continuity equation itself, with drift carrying the density along. The second is the diffusion equation itself, with noise blurring the density. Fokker–Planck is nothing more than the sum of the two equations from Part IV of this article — not something new but the sum of two things already known.

How it is derived. Computing \(\frac{d}{dt}\mathbb{E}[\varphi(X_t)]\) for an arbitrary test function \(\varphi\) via Itô's lemma gives \(\mathbb{E}[f\varphi_x+\frac{1}{2}g^2\varphi_{xx}]\). Setting this equal to \(\int\varphi\,\partial_t p\,dx\) and integrating by parts to move the derivatives onto \(p\) gives the equation above. The \(\frac{1}{2}g^2\varphi_{xx}\) term of Itô's lemma becomes the diffusion term directly — the quadratic variation carrying all the way through to here.

The stationary distribution. A distribution with \(\partial_t p=0\) may exist. If the drift is the gradient of a potential, \(f=-\nabla U\), and \(g\) is constant, the stationary distribution is \(p_\infty\propto e^{-2U/g^2}\). This is the basis on which Langevin dynamics can draw samples from an arbitrary distribution, and setting \(f=\frac{g^2}{2}\nabla\log p_{\text{target}}\) makes the stationary distribution exactly the target.

In machine learning. Because the Fokker–Planck solution corresponding to the forward SDE of a diffusion model comes out in closed form (a Gaussian), \(x_t\) can be sampled at any time \(t\) in one step and the ground-truth score is known. And as the next sections show, the fact that more than one process satisfies this PDE leads to the Probability Flow ODE.

\[\begin{aligned} \frac{\partial p}{\partial t} &= -\nabla\cdot(f p) + \frac{1}{2}g^2\nabla^2 p \\ dX = f\,dt + g\,dW &\;\longleftrightarrow\; \text{the PDE above} \\ f = -\nabla U &\;\Rightarrow\; p_\infty \propto e^{-2U/g^2} \end{aligned}\]

Turn drift and diffusion on and off separately. With drift alone the density is carried along keeping its shape; with diffusion alone it blurs in place. With both it is carried while blurring and finally stops at the stationary distribution (grey dashed). That the histogram of the particle simulation above overlays the PDE solution below is the point of this section.

The forward equation (= Fokker–Planck). An equation in the arrival variable, describing how a density spreads given an initial distribution: \(\partial_t p=-\nabla\cdot(fp)+\frac{1}{2}g^2\nabla^2 p\).

The backward equation. Setting \(u(x,t)=\mathbb{E}[\varphi(X_T)\mid X_t=x]\) gives \[-\frac{\partial u}{\partial t}=f\cdot\nabla u+\frac{1}{2}g^2\nabla^2 u,\qquad u(x,T)=\varphi(x).\] It is an equation in the departure variable, solved backwards in time from a terminal condition.

The division of labour. "Starting from this distribution now, what distribution results later?" is answered by the forward equation. "Starting from here now, what is the expectation of some later quantity?" is answered by the backward equation. The former carries densities forward, the latter expectations backward. The two operators are adjoint to each other, so a single integration by parts moves between them.

Feynman–Kac. The backward equation leads to the remarkable statement that "the solution of a PDE can be written as the expectation of a stochastic process." That is, the answer can be obtained by Monte Carlo without solving the PDE — draw many paths and average, and that is the solution. It is the only route remaining when grid methods collapse in high dimensions, and derivative pricing is computed on this principle.

The Markov property is a premise. Both equations are differentiated forms of the Chapman–Kolmogorov relation. So if the process is not Markov, this local PDE description does not hold at all.

In machine learning. That the reverse process of a diffusion model starts from a terminal condition (a Gaussian) and comes down backwards in time is the same direction as the backward equation. The Bellman equation of reinforcement learning is also a discrete-time cousin of the backward equation in propagating future expectations backward, and indeed in the continuous-time limit the same form appears as the HJB equation.

\[\begin{aligned} \text{forward: } &\partial_t p = -\nabla\cdot(fp) + \tfrac{1}{2}g^2\nabla^2 p, \quad p(\cdot,0)=p_0 \\ \text{backward: } &-\partial_t u = f\cdot\nabla u + \tfrac{1}{2}g^2\nabla^2 u, \quad u(\cdot,T)=\varphi \\ &u(x,t) = \mathbb{E}\big[\varphi(X_T)\mid X_t = x\big] \end{aligned}\]

Toggle between forward and backward. Forward spreads the initial density onward, while backward spreads the terminal condition φ back. In backward mode, click anywhere and the average of φ over paths starting there is computed, and that value agrees with the PDE solution — confirming Feynman–Kac by Monte Carlo.

The correspondence. With \(g=0\) the SDE becomes the ODE \(\dot{x}=v(x,t)\), and Fokker–Planck loses its diffusion term to become the continuity equation \(\partial_t p+\nabla\cdot(pv)=0\). To every ODE there corresponds exactly one continuity equation.

Following the log-density. Riding along with a particle and watching the density at its location gives the simple equation \(\frac{d}{dt}\log p_t(x(t))=-\nabla\cdot v\). Since divergence was the rate of volume change, this is the obvious relation that density thins in proportion to how much the volume expanded. And it is decisive in practice — only a trace need be computed, rather than a Jacobian determinant.

Exact log-likelihood. So a generative model built from a deterministic flow can compute its log-likelihood exactly. Take the log-density at the initial distribution and add the integral of \(-\nabla\cdot v\) along the path, and that is all. Continuous Normalizing Flows work on this principle, and even the trace can be approximated by a Hutchinson estimator to make it usable in high dimensions.

Deterministic, yet the density changes. Worth stressing again. Even with no randomness whatsoever in the individual particles, if \(\nabla\cdot v\ne 0\) the density spreads or gathers. "The density spreads" and "the paths are random" are separate facts. This separation is what makes the reversal of the next section possible.

In machine learning. Neural ODEs plus continuous normalizing flows, and Flow Matching, all have this structure. Flow Matching fixes a target density path first and learns by regression the \(v\) realizing it, which is possible because the continuity equation answers "which \(v\) produces that density path?"

\[\begin{aligned} \dot{x} &= v(x,t) \;\longleftrightarrow\; \partial_t p + \nabla\cdot(pv) = 0 \\ \frac{d}{dt}\log p_t(x(t)) &= -\nabla\cdot v(x(t), t) \\ \log p_T(x_T) &= \log p_0(x_0) - \int_0^T \nabla\cdot v\,dt \end{aligned}\]

The particles move deterministically and yet the density below spreads or gathers. The colour attached to each particle is the change in log-density at its location (−∇·v), and the integral in the readout is the log-likelihood correction term. That a distribution deforms without any randomness — that is the key to the next section.

Theorem. The marginal density \(p_t\) of the SDE \(dX=f\,dt+g\,dW\) coincides exactly, at every time, with the density of the following deterministic ODE: \[\frac{dx}{dt}=f(x,t)-\frac{1}{2}g^2(t)\,\nabla_x\log p_t(x).\]

Why it holds — a one-line calculation. Rewrite the diffusion term of Fokker–Planck in divergence form. Since \(\frac{1}{2}g^2\nabla^2p=\nabla\cdot\big(\frac{1}{2}g^2\nabla p\big)=\nabla\cdot\big(p\cdot\frac{1}{2}g^2\nabla\log p\big)\), we get \[\partial_t p=-\nabla\cdot\Big(p\big[f-\tfrac{1}{2}g^2\nabla\log p\big]\Big).\] This is the continuity equation for the velocity field \(v=f-\frac{1}{2}g^2\nabla\log p\). The diffusion term has been absorbed into the transport term, at the cost of making the velocity field depend on the score.

What is the same and what is different. The marginal distribution at each time is exactly the same. But the paths are entirely different — SDE paths zigzag while ODE paths are smooth. Running the two ensembles side by side in the demo shows the path shapes utterly unlike each other while the histograms overlap. If all we want is the distribution, randomness was never required.

Three practical gains. ① Fast sampling — being an ODE, higher-order solvers such as RK4 or DPM-Solver can be used, and unlike Euler–Maruyama bound to weak order 1, it stays accurate at large steps. This is the background to 1000 steps becoming 10–20. ② Exact log-likelihood — computed by the \(-\nabla\cdot v\) integral of the previous section. ③ Invertibility and a latent space — being deterministic, an image can be run back to noise, edited, and regenerated. DDIM inversion and image editing come from here.

The price. With no randomness, one initial noise corresponds to one result. Sample diversity comes only from the initial distribution, and one loses the "correcting mistakes along the way" property that SDE samplers often have (the self-correction of the Langevin term). So in practice the two are sometimes mixed, trading off quality against speed.

To summarize. This article set out from curves following a vector field, rose to the density of those curves, and then to random paths and their density. What we confirmed at the end is that, in the face of a density as the goal, determinism and randomness are interchangeable. The mathematics of diffusion models stands on that one sentence.

\[\begin{aligned} \text{SDE: } dX &= f\,dt + g\,dW \\ \text{PF-ODE: } \frac{dx}{dt} &= f - \tfrac{1}{2}g^2\,\nabla_x\log p_t(x) \\ &\Rightarrow \text{same } p_t, \text{, different paths} \end{aligned}\]

Play the two ensembles at once. The SDE paths above zigzag while the ODE paths below are smooth. Yet the two histograms on the right overlap — the readout shows the distance between the two distributions, and it approaches 0 as the number of paths grows. This is the conclusion of the article.