# Diffusion Models Generate Images Like Painters:
an Analytical Theory of Outline First, Details Later

Binxu Wang  John J. Vastola

###### Abstract

How do diffusion generative models convert pure noise into meaningful images?
We argue that generation involves first committing to an outline, and then to
finer and finer details. The corresponding reverse diffusion process can be
modeled by dynamics on a (time-dependent) high-dimensional landscape full of
Gaussian-like modes, which makes the following predictions: (i) individual
trajectories tend to be very low-dimensional; (ii) scene elements that vary
more within training data tend to emerge earlier; and (iii) early
perturbations substantially change image content more often than late
perturbations. We show that the behavior of a variety of trained unconditional
and conditional diffusion models like Stable Diffusion is consistent with
these predictions. Finally, we use our theory to search for the latent image
manifold of diffusion models, and propose a new way to generate interpretable
image variations. Our viewpoint suggests generation by GANs and diffusion
models have unexpected similarities.

diffusion model, geometry

Figure 1: Characteristics of image generation by diffusion models. A. Tracking
latent states $G(\boldsymbol{\mathbf{x}}_{t})$ (top row), differences between
nearby time steps
$G(k(\boldsymbol{\mathbf{x}}_{t-1}-\boldsymbol{\mathbf{x}}_{t}))$ (middle
row), and final image estimates
$G(\hat{\boldsymbol{\mathbf{x}}}_{0}(\boldsymbol{\mathbf{x}}_{t}))$ (bottom
row) suggests different measures of progress. B. Individual trajectories are
effectively two-dimensional, with the transition from
$\boldsymbol{\mathbf{x}}_{T}$ to $\boldsymbol{\mathbf{x}}_{0}$ being rotation-
like.

##  1 Introduction

Imagine an artist painting a picture of a natural landscape. We generally
expect higher-level scene elements to appear first, and lower-level details
later: the borders of the land and sky might be drawn, then the largest
objects (like mountains and trees) might be placed, then minor objects (like
rocks and small animals) might be placed, and finally fine details (like
textures and shading) might be filled in. Do diffusion generative models
(Sohl-Dickstein et al., 2015; Song & Ermon, 2019; Song et al., 2021), which
can generate natural landscapes like those of an artist, also construct images
like this? If not, how do they work?

By visualizing images throughout a reverse diffusion trajectory—from pure
noise to the final image—the naive answer appears to be no. One gets the
impression of an image emerging fully-formed from the noise; one is
‘uncovering’ the image, or ‘opening one’s eyes’ to reveal an image that was
always there. But this naive answer is highly misleading, and we will later
argue that it is entirely wrong.

One significant step towards understanding this question was provided by Hertz
et al., who found that models’ cross-attention maps substantially control
scene composition (2022). They noted that different parts of an image may be
‘attended to’ at different times in a reverse diffusion trajectory. But it is
unclear if attention mechanisms are necessary for image generation to be like
painting, or if this is a generic property of diffusion models.

Understanding image generation is important for a few reasons. As Hertz et al.
show, improving understanding can make the process more controllable, and
hence more precise, efficient, and useful. These improvements may be crucial
for generalizing image generation to stable video generation, as some early
work has attempted (Ho et al., 2022). Beyond AI, these details could have some
bearing on the brain; in the same way convolutional neural networks (CNNs) are
used to test hypotheses about how early visual areas in the brain might
behave, diffusion models may also be useful as a toy model of human visual
imagination. Some research in psychology suggests that humans imagine visual
scenes hierarchically, with important scene elements imagined first and
details imagined later (Kosslyn & Shwartz, 1977).

Motivated by a variety of salient observations, in this paper we propose a
simple but insightful theory of reverse diffusion based on approximating image
distributions as Gaussian mixtures. It suggests an outline-first, details-
later view of image generation that begins with a prototype and gradually
refines it. This theoretical view yields several interesting predictions for
how real generation proceeds, including that: (i) individual trajectories tend
to be very low-dimensional; (ii) scene elements that vary more within training
data tend to emerge earlier; and (iii) early perturbations substantially
change image content more often than late perturbations.

This view also sheds light on an open problem: what is the nature of the low-
dimensional latent manifold that, at least in principle, smoothly
parameterizes diffusion-generated images? For generative adversarial networks
(GANs) (Goodfellow et al., 2014), the geometry of this manifold is fairly
well-understood (Wang & Ponce, 2021); for diffusion models, it is unclear
‘where’ it even is. Recently, Wenliang and Moran searched for it by dissecting
the Jacobian of the score vector field (2022). Here, we claim on-manifold
directions that yield meaningful image perturbations can be identified via
principal component analysis (PCA) on trajectories.

The paper is organized as follows. In Section 2, we provide necessary
background on diffusion models. In Section 3, we collect salient observations
about image generation, e.g. about trajectory geometry and feature emergence
order. In Section 4, we explain these observations using an analytical model
of reverse diffusion, and explore some of the model’s consequences. In Section
5, we test these predictions in a variety of trained unconditional diffusion
models. In Section 6, we apply our theoretical view to efficiently generate
interpretable variations of Stable Diffusion images.

##  2 Diffusion generative modeling basics

There are several complementary theoretical viewpoints on diffusion generative
modeling (Sohl-Dickstein et al., 2015; Song & Ermon, 2019; Yang et al., 2022);
here, we review the details of Song et al.’s unifying view (2021).

Diffusion generative models involve mapping one probability distribution
$p(\boldsymbol{\mathbf{x}})$ to another via a stochastic process—typically
pure diffusion or an Ornstein-Uhlenbeck (OU) process. This so-called ‘forward’
process can be inverted by considering its ‘reverse’ process, which is
mathematically guaranteed to exist for reasonable choices of initial
distribution and forward process (Anderson, 1982). If
$p(\boldsymbol{\mathbf{x}})$ is mapped to a distribution that is easy to
sample from, like a normal distribution, by sampling from it and running the
reverse process one can obtain samples from $p(\boldsymbol{\mathbf{x}})$.
(Details follow; Appendix D relates our notation to others’ notation.)

##### Forward/reverse diffusion.

A fairly general forward process is given by the stochastic differential
equation (SDE)

$\dot{\boldsymbol{\mathbf{x}}}=\boldsymbol{\mathbf{f}}(\boldsymbol{\mathbf{x}})+g(t)\
\boldsymbol{\mathbf{\eta}}(t)$ (1)

where $g(t)$ is a time-dependent noise amplitude,
$\boldsymbol{\mathbf{\eta}}(t)$ is a vector of independent Gaussian white
noise terms, and time runs from $t=0$ to $t=T=1$. The reverse process is

$\dot{\boldsymbol{\mathbf{x}}}=\boldsymbol{\mathbf{f}}(\boldsymbol{\mathbf{x}})-g(t)^{2}\boldsymbol{\mathbf{s}}(\boldsymbol{\mathbf{x}},t)+g(t)\
\boldsymbol{\mathbf{\eta}}(t)$ (2)

where
$\boldsymbol{\mathbf{s}}(\boldsymbol{\mathbf{x}},t):=\nabla_{\boldsymbol{\mathbf{x}}}\log
p(\boldsymbol{\mathbf{x}},t)$ is the score function, and where we use the
standard convention that times runs backward, i.e. from $t=1$ to $t=0$. In
this paper, we focus on one common forward process: the variance-preserving
SDE

$\dot{\boldsymbol{\mathbf{x}}}=-\beta(t)\boldsymbol{\mathbf{x}}+g(t)\boldsymbol{\mathbf{\eta}}(t)$
(3)

where we enforce the constraint that $\beta(t)=\frac{1}{2}g^{2}(t)$. The
marginal probabilities of this process are

$p(\boldsymbol{\mathbf{x}}_{t}|\boldsymbol{\mathbf{x}}_{0})=\mathcal{N}(\boldsymbol{\mathbf{x}}_{t}|\alpha_{t}\boldsymbol{\mathbf{x}}_{0},\sigma_{t}^{2}\boldsymbol{\mathbf{I}})$
(4)

where $\alpha_{t}$ and $\sigma_{t}^{2}$ satisfy
$\alpha_{t}^{2}+\sigma_{t}^{2}=1$ and are defined to be

$\begin{split}\alpha_{t}:=e^{-\int_{0}^{t}\beta(t^{\prime})dt^{\prime}}\hskip
14.45377pt\sigma_{t}^{2}:=1-e^{-2\int_{0}^{t}\beta(t^{\prime})dt^{\prime}}\
.\end{split}$ (5)

In practice, to reverse an SDE like this, we do not simulate the corresponding
reverse process; instead, we simulate a probability flow ODE with the same
marginal probabilities (Song et al., 2021). To reverse Eq. 3, we would use

$\dot{\boldsymbol{\mathbf{x}}}=-\beta(t)\boldsymbol{\mathbf{x}}-\frac{1}{2}g^{2}(t)\
\boldsymbol{\mathbf{s}}(\boldsymbol{\mathbf{x}},t)$ (6)

where time again runs backward from $t=1$ to $t=0$.

##### Learning the score function.

The score function, which is required to reverse the forward process, can be
learned via gradient descent on a score-matching objective

$\mathbb{E}_{\boldsymbol{\mathbf{x}}_{0}\sim
p(\boldsymbol{\mathbf{x}}_{0}),\boldsymbol{\mathbf{\epsilon}}\sim\mathcal{N}(\boldsymbol{\mathbf{0}},\boldsymbol{\mathbf{I}})}\int_{0}^{1}\gamma_{t}\|\boldsymbol{\mathbf{\epsilon}}_{\theta}(\alpha_{t}\boldsymbol{\mathbf{x}}_{0}+\sigma_{t}\boldsymbol{\mathbf{\epsilon}},t)-\boldsymbol{\mathbf{\epsilon}}\|_{2}^{2}\
dt$ (7)

where
$\boldsymbol{\mathbf{\epsilon}}_{\theta}(\boldsymbol{\mathbf{x}}_{t},t)\approx-\sigma_{t}\nabla_{\boldsymbol{\mathbf{x}}}\log
p(\boldsymbol{\mathbf{x}}_{t})$ is a neural network approximation of the
negative scaled score function, and $\gamma_{t}$ is a positive weighting
function (Song et al., 2021).

##### Latent diffusion models.

Pixel space representations of images are extremely high-dimensional, and
movement through this space often corresponds to changing imperceptible high-
frequency details, making diffusion inefficient. One alternative uses an auto-
encoder ($\boldsymbol{\mathbf{x}}\leftrightarrow G(\boldsymbol{\mathbf{x}})$)
to map images into an expressive latent space, and trains a diffusion model in
latent rather than pixel space. Models like Stable Diffusion leverage this
efficiency increase to scale to massive data sets at a modest cost (Rombach et
al., 2022).

##### DDIM/PNDM samplers.

Sampling the reverse process is somewhat independent of score function
learning, enabling researchers to separately study its efficiency. Ho et al.
use the DDPM sampler (2020), which is effectively the same as discretizing a
reverse SDE (Eq. 2). Song et al. propose the DDIM sampler (2020), a
deterministic sampler equivalent to solving the probability flow ODE (c.f. Eq.
6), which dramatically reduced the required number of steps. The idea of
sampling via solving an ODE enables researchers to exploit various numerical
methods for integrating ODEs; one powerful example is Liu et al.’s PNDM
sampler (2022), the current default sampler for Stable Diffusion. In this
work, we focus on the DDIM/PNDM samplers, and comment on how other samplers
affect our results in Section B.6.

##  3 Salient observations about image generation

###  3.1 How should we measure generation progress?

A common way to monitor image generation progress is to observe how
$\boldsymbol{\mathbf{x}}_{t}$ (or, in the case of latent diffusion, the
decoded image $G(\boldsymbol{\mathbf{x}}_{t})$ that corresponds to it) changes
over time. As previously mentioned, this approach tends to show a fully-formed
image unveiled from noise (Fig. 1A, top row). But is this what is ‘actually’
happening?

A simple but useful alternative is to observe (appropriately scaled)
differences $k(\boldsymbol{\mathbf{x}}_{t-k}-\boldsymbol{\mathbf{x}}_{t})$
between close time points—the next ‘layer of paint’ that has been added to the
canvas. These often appear to be like a sketch of the final outcome early in
generation, and are increasingly contaminated by noise towards the end (Fig.
1A, middle row).

A more principled alternative, first proposed by Ho et al. (2020), is to
consider the sequence of endpoint estimates
$\hat{\boldsymbol{\mathbf{x}}}_{0}$ of the reverse diffusion trajectory. In
particular, the weighted combination of the latent state and score given by

$\hat{\boldsymbol{\mathbf{x}}}_{0}(\boldsymbol{\mathbf{x}}(t)):=\frac{\boldsymbol{\mathbf{x}}(t)-\sigma_{t}\boldsymbol{\mathbf{\epsilon}}_{\theta}(\boldsymbol{\mathbf{x}}(t),t)}{\alpha_{t}}$
(8)

provides an endpoint estimate that improves over time. We found tracking this
statistic throughout reverse diffusion provides substantial insight into the
generation process: as early as the first time step, a rough outline is often
visible. As time goes on, one tends to see progressively finer details filled
in (Fig. 1A, bottom). Intriguingly, the endpoint estimate
$\hat{\boldsymbol{\mathbf{x}}}_{0}(\boldsymbol{\mathbf{x}}_{t})$ is initially
similar to the scaled state difference
$k(\boldsymbol{\mathbf{x}}_{t-k}-\boldsymbol{\mathbf{x}}_{t})$. While Fig. 1
depicts an example of conditional generation in a latent diffusion model, we
observed similar results for many unconditional diffusion models. For MNIST,
the estimate $\hat{\boldsymbol{\mathbf{x}}}_{0}$ is initially ambiguous, then
starts to resemble the digit eventually generated (Fig. 9). For the relatively
unimodal face data set CelebA-HQ, the estimate appears to be interpretable
even at the first time step (Fig.7)—specifically, it is an average face on a
generic background (Langlois & Roggman, 1990).

###  3.2 When do different image features tend to emerge?

What is the order of feature emergence according to the endpoint estimate
$\hat{\boldsymbol{\mathbf{x}}}_{0}(\boldsymbol{\mathbf{x}}_{t})$? As Ho et al.
(2020) note, high-level features tend to emerge before low-level ones. For
example, when generating “a portrait of an aristocrat” (Fig. 1A), a generic
face-like shape appears, and then is refined to include blobs that correspond
to hat/hair and body. Coarse facial features and various image colors emerge,
facial hair appears, and the blob above the face is ‘reinterpreted’ into a
hat. Finally, high-frequency details of the face and clothes are added.
Similar findings are observed in unconditional models (e.g. of faces,
churches, etc.; Fig. 7,8).

###  3.3 What is the shape of individual trajectories?

We also studied the geometry of reverse diffusion trajectories. Although the
dimensionality of image/latent space is quite large, individual trajectories
are effectively two-dimensional: the average variance explained by the top two
principal components (PCs) is 99.98% for our CelebA model, and 99.54% for
Stable Diffusion (see Table B.2 and Fig. 11). To good approximation,
$\boldsymbol{\mathbf{x}}_{t}$ always remains in the plane defined by the
initial noise $\boldsymbol{\mathbf{x}}_{T}$ and the final state
$\boldsymbol{\mathbf{x}}_{0}$: the variance explained by a projection onto
this 2D plane is higher than 99.2% for all models. Furthermore, the reverse
diffusion trajectory is well-approximated by a rotation within this plane
(Fig. 1B), i.e.

$\boldsymbol{\mathbf{x}}_{t}\approx
K_{t}\boldsymbol{\mathbf{x}}_{0}+\sqrt{1-K_{t}^{2}}\
\boldsymbol{\mathbf{x}}_{T}$ (9)

where $0\leq K_{t}\leq 1$, $K_{T}=0$, and $K_{0}=1$. When $K_{t}=\alpha_{t}$,
this equation explains almost all trajectory variance: 97.51% for Stable
Diffusion and 98.93% for CelebA diffusion.

##  4 Theoretical analysis of reverse trajectories

How should we understand the previous observations? The explicit form of the
probability flow ODE (Eq. 6) suggests an obvious but surprisingly powerful
geometric view of image generation: a particle climbs down the changing
landscape defined by the minus log-likelihood function and forward process.
Here, we describe a simplified but useful version of this picture that both
informs our intuition, and makes quantitative predictions about reverse
diffusion.

Figure 2: Conceptual view of unimodal landscape dynamics.

###  4.1 Summary of landscape dynamics view of generation

Assume that the training set of images is well-approximated by a vast
collection of Gaussian-like modes. The forward process smears these modes out,
repeatedly grouping nearby modes until only one large and nondescript mode
remains. The reverse process does exactly the opposite: one mode splits into
several, those modes split into more, and so on, until one recovers a Gaussian
mixture like the training set.

Image generation can be characterized by two related processes: the dynamics
of a particle driven by the constantly changing shape of this landscape, and
the dynamics of the continuously updated estimate of that particle’s endpoint.
Particle dynamics correspond to the conversion of noise into an image, while
estimate dynamics represent the layer-by-layer (outline-first, details-later)
painting of that image.

To understand these processes, we must know: (i) their dynamics on short time
scales, when the number of modes is fixed; (ii) the frequency and effects of
mode-splitting. Both aspects cause higher-level features to appear earlier.

###  4.2 One Gaussian mode: dynamics of particle

We will begin by studying a simpler problem: reverse diffusion given a
training set well-described by a single Gaussian mode. In this setting, we can
be quantitatively precise. In Sections 4.4 and 4.5, we will show how to apply
this theory to the general multi-mode setting.

Suppose that the mean and covariance of this mode are
$\boldsymbol{\mathbf{\mu}}$ and $\boldsymbol{\mathbf{\Sigma}}$. The score
function at time $t>0$ is the score of a Gaussian
$\mathcal{N}(\alpha_{t}\boldsymbol{\mathbf{\mu}},\sigma_{t}^{2}\boldsymbol{\mathbf{I}}+\alpha_{t}^{2}\boldsymbol{\mathbf{\Sigma}})$,
which is

$\boldsymbol{\mathbf{s}}(\boldsymbol{\mathbf{x}},t)=\nabla_{\boldsymbol{\mathbf{x}}}p(\boldsymbol{\mathbf{x}},t)=(\sigma_{t}^{2}\boldsymbol{\mathbf{I}}+\alpha_{t}^{2}\boldsymbol{\mathbf{\Sigma}})^{-1}(\alpha_{t}\boldsymbol{\mathbf{\mu}}-\boldsymbol{\mathbf{x}})\
.$ (10)

The corresponding probability flow ODE

$\dot{\boldsymbol{\mathbf{x}}}=-\beta(t)\boldsymbol{\mathbf{x}}-\beta(t)(\sigma_{t}^{2}\boldsymbol{\mathbf{I}}+\alpha_{t}^{2}\boldsymbol{\mathbf{\Sigma}})^{-1}(\alpha_{t}\boldsymbol{\mathbf{\mu}}-\boldsymbol{\mathbf{x}})$
(11)

can be solved exactly (Appendix G).
$\boldsymbol{\mathbf{y}}:=\boldsymbol{\mathbf{x}}-\alpha_{t}\boldsymbol{\mathbf{\mu}}$
satisfies

$-\dot{\boldsymbol{\mathbf{y}}}/\beta(t)=-\left[(\sigma_{t}^{2}\boldsymbol{\mathbf{I}}+\alpha_{t}^{2}\boldsymbol{\mathbf{\Sigma}})^{-1}-\boldsymbol{\mathbf{I}}\right]\boldsymbol{\mathbf{y}}\
$ (12)

The covariance matrix has a compact singular value decomposition (SVD)
$\boldsymbol{\mathbf{\Sigma}}=\boldsymbol{\mathbf{U}}\boldsymbol{\mathbf{\Lambda}}\boldsymbol{\mathbf{U}}^{T}$,
where
$\boldsymbol{\mathbf{U}}=[\boldsymbol{\mathbf{u}}_{1},...,\boldsymbol{\mathbf{u}}_{r}]$
is a $D\times r$ semi-orthogonal matrix. The columns of
$\boldsymbol{\mathbf{U}}$ are the principal axes along which the mode varies,
and together comprise an ‘image manifold’. If we define the matrix

$[\tilde{\boldsymbol{\mathbf{\Lambda}}}_{t}]_{ij}:=\delta_{ij}\frac{\alpha_{t}^{2}\lambda_{i}}{\sigma_{t}^{2}+\alpha_{t}^{2}\lambda_{i}}\
,$ (13)

where $\lambda_{i}:=\boldsymbol{\mathbf{\Lambda}}_{ii}$, we can rewrite Eq. 12
as

$-\dot{\boldsymbol{\mathbf{y}}}/\beta(t)=\frac{1}{\sigma_{t}^{2}}\left[-{\alpha_{t}^{2}}\boldsymbol{\mathbf{I}}+\boldsymbol{\mathbf{U}}\tilde{\boldsymbol{\mathbf{\Lambda}}}_{t}\boldsymbol{\mathbf{U}}^{T}\right]\boldsymbol{\mathbf{y}}\
.$ (14)

Eq. 12 is linear, and the dynamics of $\boldsymbol{\mathbf{y}}_{t}$ along each
principal axis are independent of dynamics along the others. We can write
$\boldsymbol{\mathbf{y}}_{t}$ as a sum of on- and off-manifold components:

$\begin{split}\boldsymbol{\mathbf{y}}_{t}&=\psi(t,0)\boldsymbol{\mathbf{y}}^{\perp}_{T}+\sum_{k=1}^{r}\psi(t,\lambda_{k})c_{k}(T)\boldsymbol{\mathbf{u}}_{k}\\\
c_{k}(T)&=\boldsymbol{\mathbf{u}}_{k}^{T}(\boldsymbol{\mathbf{x}}_{T}-\alpha_{T}\boldsymbol{\mathbf{\mu}})\\\
\boldsymbol{\mathbf{y}}^{\perp}_{T}&=(\boldsymbol{\mathbf{I}}-\boldsymbol{\mathbf{U}}\boldsymbol{\mathbf{U}}^{T})(\boldsymbol{\mathbf{x}}_{T}-\alpha_{T}\boldsymbol{\mathbf{\mu}})\\\
\psi(t,\lambda)&:=\sqrt{\frac{1+(\lambda-1)\alpha_{t}^{2}}{1+(\lambda-1)\alpha_{T}^{2}}}\end{split}$
(15)

where $T>0$ marks the beginning of dynamics. Since
$\alpha_{t}\xrightarrow{t\to 0}1$, the components of $\boldsymbol{\mathbf{y}}$
orthogonal to the image manifold decay toward zero. The coefficients of on-
manifold components grow if $\lambda_{k}<1$, remain fixed if $\lambda_{k}=1$,
and shrink if $\lambda_{k}>1$; in all cases, $\psi(t;\lambda_{k})$ approaches
$\sqrt{\lambda_{k}}$, the standard deviation along
$\boldsymbol{\mathbf{u}}_{k}$, as $t\to 0$.

Eq. 15 implies high variance image manifold features (i.e. directions with
high $\lambda_{k}$) appear first, and that low variance features appear later.
It also implies an increasing commitment to image features as time goes on.
Suppose at time $t<t^{\prime}<T$ the off-manifold directions are perturbed by
$\delta\boldsymbol{\mathbf{y}}^{\perp}$, and the on-manifold direction
coefficients are perturbed by amounts $\delta c_{k}$; the total differences
$\Delta\boldsymbol{\mathbf{y}}^{\perp}(t)$ and $\Delta c_{k}(t)$ in these
components at time $t$ are

$\begin{split}\Delta\boldsymbol{\mathbf{y}}^{\perp}_{t}&=\delta\boldsymbol{\mathbf{y}}^{\perp}\frac{\psi(t,0)}{\psi(t^{\prime},0)}=\delta\boldsymbol{\mathbf{y}}^{\perp}\sqrt{\frac{1-\alpha_{t}^{2}}{1-\alpha_{t^{\prime}}^{2}}}\\\
\Delta c_{k}(t)&=\delta
c_{k}\frac{\psi(t,\lambda_{k})}{\psi(t^{\prime},\lambda_{k})}=\delta
c_{k}\sqrt{\frac{1+(\lambda_{k}-1)\alpha_{t}^{2}}{1+(\lambda_{k}-1)\alpha_{t^{\prime}}^{2}}}\
.\end{split}$ (16)

Both on- and off-manifold perturbations have smaller effects as $t^{\prime}$
gets closer to $t$, since $\alpha_{t^{\prime}}$ is closer to $\alpha_{t}$.
Off-manifold perturbations decay, while on-manifold perturbations have a
lasting effect (and perturbations in high variance directions have larger
effects). Eq. 15 also implies (see Appendix H)

$\begin{split}\boldsymbol{\mathbf{x}}_{t}\approx\
&\alpha_{t}\boldsymbol{\mathbf{x}}_{0}+\sqrt{1-\alpha_{t}^{2}}\
\boldsymbol{\mathbf{x}}_{T}\\\
&+\sum_{k=1}^{r}\left\\{\sqrt{\sigma_{t}^{2}+\lambda_{k}\alpha_{t}^{2}}-\alpha_{t}\sqrt{\lambda_{k}}-\sigma_{t}\right\\}c_{k}(T)\boldsymbol{\mathbf{u}}_{k}\
,\end{split}$

which can be used to explain our observation that
$\boldsymbol{\mathbf{x}}_{t}$ dynamics tend to look like a rotation within a
2D plane (Fig. 1B). It is a rotation when the sum term vanishes, which happens
when $\alpha_{t}\approx 0$ and $\alpha_{t}\approx 1$. The sum term
approximately vanishes when $\lambda_{k}$ is sufficiently large or small—but
in general, it does not, meaning $\boldsymbol{\mathbf{x}}_{t}$ is only
approximately a rotation. Interestingly, since $\alpha_{t}$ does not depend on
$\boldsymbol{\mathbf{\mu}}$ or $\boldsymbol{\mathbf{\Sigma}}$, the rotation
does not depend on any properties of the mode. We empirically find, however,
that it does depend on the sampler (Section B.6.1).

###  4.3 One Gaussian mode: dynamics of endpoint estimate

The endpoint estimate $\hat{\boldsymbol{\mathbf{x}}}_{0}$ (see Eq. 8 and
Appendix F) is

$\hat{\boldsymbol{\mathbf{x}}}_{0}(\boldsymbol{\mathbf{x}}_{t}):=\frac{\boldsymbol{\mathbf{x}}_{t}+\sigma_{t}^{2}(\sigma_{t}^{2}\boldsymbol{\mathbf{I}}+\alpha_{t}^{2}\boldsymbol{\mathbf{\Sigma}})^{-1}(\alpha_{t}\boldsymbol{\mathbf{\mu}}-\boldsymbol{\mathbf{x}}_{t})}{\alpha_{t}}\
.$ (17)

Initially, when $\sigma_{t}^{2}\sim 1$, we have
$\hat{\boldsymbol{\mathbf{x}}}_{0}\approx\boldsymbol{\mathbf{\mu}}$, which we
can interpret as something like the ‘average face’ (Langlois & Roggman, 1990)
if the training set consists of face images. Later, this might look like the
average of a more specialized set of faces, e.g. the average older male face.
Close to the end of reverse diffusion, when $\sigma_{t}^{2}\sim 0$, we have
$\hat{\boldsymbol{\mathbf{x}}}_{0}\approx\boldsymbol{\mathbf{x}}_{t}$, i.e.
the endpoint estimate and true trajectory endpoint match.

Using Eq. 15, note that $\hat{\boldsymbol{\mathbf{x}}}_{0}$ can be written

$\begin{split}\hat{\boldsymbol{\mathbf{x}}}_{0}(\boldsymbol{\mathbf{x}}_{t})&=\boldsymbol{\mathbf{\mu}}+\left[\boldsymbol{\mathbf{I}}-\sigma_{t}^{2}(\sigma_{t}^{2}\boldsymbol{\mathbf{I}}+\alpha_{t}^{2}\boldsymbol{\mathbf{\Sigma}})^{-1}\right]\boldsymbol{\mathbf{y}}_{t}/\alpha_{t}\\\
&=\boldsymbol{\mathbf{\mu}}+\boldsymbol{\mathbf{U}}\tilde{\boldsymbol{\mathbf{\Lambda}}}_{t}\boldsymbol{\mathbf{U}}^{T}\boldsymbol{\mathbf{y}}_{t}/\alpha_{t}\\\
&=\boldsymbol{\mathbf{\mu}}+\sum_{k=1}^{r}\xi(t,\lambda_{k})c_{k}(T)\boldsymbol{\mathbf{u}}_{k}\\\
\xi(t,\lambda)&:=\frac{\alpha_{t}\lambda}{\sqrt{(\sigma_{t}^{2}+\lambda\alpha_{t}^{2})(\sigma_{T}^{2}+\lambda\alpha_{T}^{2})}}\end{split}$
(18)

which quantifies the fraction of $\hat{\boldsymbol{\mathbf{x}}}_{0}$ captured
by the mean $\mu$ and each of the principal component. Interestingly, this
equation implies that $\hat{\boldsymbol{\mathbf{x}}}_{0}$ always remains on
the image manifold, which explains why visualized estimates always look like
images, rather than images contaminated with noise.

Eq. 18 implies that image manifold features appear in
$\hat{\boldsymbol{\mathbf{x}}}_{0}$ in order of descending variance: first the
highest variance features appear, then the next highest, and so on. There is
also increasing commitment to features that appear earlier, in the sense that
perturbations are less likely to change them later in reverse diffusion. To
quantify this, we can derive an equation analogous to Eq. 16:

$\begin{split}\Delta\hat{\boldsymbol{\mathbf{x}}}_{0}&=\sum_{k=1}^{r}\frac{\alpha_{t}\lambda_{k}}{\sqrt{\left[\sigma_{t}^{2}+\alpha_{t}^{2}\lambda_{k}\right]\left[\sigma_{t^{\prime}}^{2}+\alpha_{t^{\prime}}^{2}\lambda_{k}\right]}}\delta
c_{k}\boldsymbol{\mathbf{u}}_{k}\end{split}$ (19)

where $\Delta\hat{\boldsymbol{\mathbf{x}}}_{0}$ is the change in
$\hat{\boldsymbol{\mathbf{x}}}_{0}$ at time $t$ due to perturbations $\delta
c_{k}$ at time $t^{\prime}$. Our conceptual view of particle and estimate
dynamics in a unimodal landscape is depicted in Fig. 2.

###  4.4 Short time scale dynamics in general

Moving beyond the unimodal case, on short time scales the landscape looks like
a Gaussian mixture whose modes are shrinking in size. A priori, we might
expect that multiple modes influence the particle’s movement at any given
time. However, a quirk of high-dimensional Gaussians is that almost all of
their probability is concentrated in a thin shell around their mode (see
Appendix E), and we can use this fact to approximate the score. Note that if
$p(\boldsymbol{\mathbf{x}},t)$ represents a Gaussian mixture whose components
have weights $\pi_{k}$, means $\boldsymbol{\mathbf{\mu}}_{k}(t)$, and
covariance matrices $\boldsymbol{\mathbf{\Sigma}}_{k}(t)$,

$\boldsymbol{\mathbf{s}}(\boldsymbol{\mathbf{x}},t)=\frac{\sum_{k}\pi_{k}\boldsymbol{\mathbf{\Sigma}}^{-1}(t)(\boldsymbol{\mathbf{\mu}}(t)-\boldsymbol{\mathbf{x}})\mathcal{N}(\boldsymbol{\mathbf{x}};\boldsymbol{\mathbf{\mu}}_{k}(t),\boldsymbol{\mathbf{\Sigma}}_{k}(t))}{\sum_{k}\pi_{k}\mathcal{N}(\boldsymbol{\mathbf{x}};\boldsymbol{\mathbf{\mu}}_{k}(t),\boldsymbol{\mathbf{\Sigma}}_{k}(t))}\
.$

This is essentially a softmax function, with each mode contributing based on
size and proximity to $\boldsymbol{\mathbf{x}}$. The sparseness of high-
dimensional Gaussians means that only the mode closest to
$\boldsymbol{\mathbf{x}}$ ($k^{\prime}$, say) will contribute appreciably to
the score, so that
$\boldsymbol{\mathbf{s}}(\boldsymbol{\mathbf{x}},t)\approx\boldsymbol{\mathbf{\Sigma}}_{k^{\prime}}(\boldsymbol{\mathbf{\mu}}_{k^{\prime}}-\boldsymbol{\mathbf{x}})$
to very good approximation. Hence, we can use the theory from the previous two
sections to accurately describe short time scale dynamics.

Figure 3: Conceptual view of mode-splitting.

Figure 4: Analytical predictions of convergence speed along dimensions of
different variances. A. The $\psi(t,\lambda_{k})$ function, with different
$\lambda_{k}=0.0,0.01,0.1,1,10,100$. It governs the dynamics of
$\boldsymbol{\mathbf{y}}_{t}$ along eigendimension
$\boldsymbol{\mathbf{u}}_{k}$ with variance $\lambda_{k}$. B. Time derivative
of the curves in A., i.e. the rate of change along each eigendimension. C. The
$\xi(t,\lambda_{k})$ function governing the dynamics of
$\hat{\boldsymbol{\mathbf{x}}}_{0}(\boldsymbol{\mathbf{x}}_{t})-\boldsymbol{\mathbf{\mu}}$
along eigendimension $\boldsymbol{\mathbf{u}}_{k}$, normalized by the standard
deviation $\sqrt{\lambda_{k}}$, showing the course of convergence for the
projected outcome. D. The time derivative of the curve in C., showing the
relative rates of change of
$\hat{\boldsymbol{\mathbf{x}}}_{0}(\boldsymbol{\mathbf{x}}_{t})-\boldsymbol{\mathbf{\mu}}$
along each eigendimension, highlighting that the projected outcome first moves
along directions of larger variance. The specific $\alpha_{t}$ schedule is
copied from ddpm-CIFAR-10 (see Appendix A). Other schedules yield similar
results.

###  4.5 Frequency and effects of mode-splitting

On longer time scales, the number of modes tends to increase exponentially: as
the noise scale gets smaller, the landscape reflects finer and finer details
of the training set, causing one mode to become two, two to become four, and
so on. From the particle’s point of view, the mean and covariance matrix of
the closest mode suddenly switches, then single mode dynamics runs until the
next splitting event.

To make our intuition more precise, consider a training set of $M$ closely
packed $D$-dimensional isotropic Gaussian modes. Given a coarse asymptotic
bound on high-dimensional sphere packing (Ball, 1992), we can estimate that a
given mode will have roughly $2^{D}$ nearest neighbors. At the beginning of
reverse diffusion, one mode splits into $2^{D}$ modes. Each of these splits
into $2^{D}$ more modes, each of those splits, and so on; after $k$
splittings, there are $2^{kD}$ modes. Hence, on the order of $\log M$
splittings are required to recover the resolution of the training set.

Splittings occur when the noise scale goes below the variances of the current
modes (equivalently: the distance between the underlying modes), uncovering
modes with variances that are a factor of roughly $2^{-D}$ smaller. After $k$
splittings, they are a factor of $2^{-kD}$ smaller. Because $\beta(t)$ is
usually chosen to depend linearly on time, the noise scale $\sigma_{t}$
changes exponentially with time, which means splittings should occur at
regular time intervals, and the characteristic amount of time between
splittings should scale as $1/\log M$.

According to this model, image generation should look hierarchical. High-level
features, which correspond to a smeared out collection of alike modes, are
committed to first. As these modes split, the particle ends up in the
attractor basin of one of the component modes, and the image commits to lower-
level features. This process continues until only high-frequency details
remain (Fig. 3).

This commitment to a particular collection of modes yields another reason that
later perturbations have less of an effect: they tend to push the particle
into nearby modes, which correspond to similar images. The typical distances
between particles that committed to different modes is expected to increase
exponentially with divergence time; for that reason, the strength of the
perturbation required to push them into nearby basins of attraction also
increases exponentially.

###  4.6 Numerical simulation

To sharpen our intuition about our analytical results (especially Eq. 15 and
18), we used a common $\alpha_{t}$ schedule (see Appendix A) and plotted the
function $\psi(t,\lambda_{k})$ and $\psi(t,0)$ which govern the
$\boldsymbol{\mathbf{x}}_{t}$ dynamics, for different $\lambda_{k}$, i.e. the
variance of the distribution in the direction $\boldsymbol{\mathbf{u}}_{k}$.
The projection coefficients of $\boldsymbol{\mathbf{x}}_{t}$ along high
variance directions increase at first, while the coefficients along the low
variance directions remain the same or shrink later (Fig. 4A).

By plotting the derivative of $\psi(t,\lambda)$ function, we can examine
$\dot{\boldsymbol{\mathbf{x}}}_{t}$ i.e. the velocity of particle movement
(Fig. 4B). It is initially dominated by contributions from the large variance
dimensions. Close to the end, it features comparable but opposite-sign
contributions from the high variance dimensions and the off-manifold (i.e.
noise) directions. This explains the observation that state differences
$\boldsymbol{\mathbf{x}}_{t-k}-\boldsymbol{\mathbf{x}}_{t}$ look interpretable
at first, and later look like noise-contaminated images (Fig. 1A, middle row).

Finally, we examined $\xi(t,\lambda_{k})$ which governs the dynamics of
$\hat{\boldsymbol{\mathbf{x}}}_{0}(\boldsymbol{\mathbf{x}}_{t})$ along each
principal direction $\boldsymbol{\mathbf{u}}_{k}$. These functions look
sigmoidal and converge to 1 when normalized by the standard deviation
$\sqrt{\lambda_{k}}$ (Fig. 4C). The functions converge earlier along large
variance dimensions, and later along small variance dimensions. Furthermore,
their time derivatives $\frac{d}{dt}\xi(t,\lambda_{k})$ are bump-like, and
peak earlier for high variance dimensions. In plain language, high variance
features are added to the endpoint estimate first, and low variance features
are added later. This partly explains our earlier observations about the order
of feature emergence (Fig. 1A, bottom row).

##  5 Validating theory on unconditional models

Even if approximating the training set as a high-dimensional Gaussian mixture
is valid, it is difficult to make detailed quantitative predictions about the
effects of mode-splitting on reverse diffusion dynamics. As a first step, we
can instead approximate the training set as a single Gaussian mode, and test
our earlier predictions. To do this, for a few accessible image data
sets—MNIST, CIFAR-10, and CelebA—we numerically computed the mean and
covariance of all images in pixel space. We then investigated three
unconditional pixel space diffusion models trained on these data sets (see
Appendix A). Although crude, this approximation is not unreasonable, since
e.g. face space may be well-approximated by a multivariate Gaussian (Turk &
Pentland, 1991).

##### Coarse agreement between theory and actual trajectory.

Under the single mode approximation, we can predict how
$\boldsymbol{\mathbf{x}}_{t}$ and
$\hat{\boldsymbol{\mathbf{x}}}_{0}(\boldsymbol{\mathbf{x}}_{t})$ will evolve.
In many cases, we can accurately describe the early part of the trajectory and
predict the low-frequency layout of the final image well (Fig. 5A-B). High-
frequency details like edges tend not to be well-predicted. When we compared
the theory and actual trajectory through time, we found the deviation blows up
around 20 steps (Fig. 5C); we think this is a signature of when the single
static mode assumption breaks down, and the trajectory starts to be guided by
the updated fine-grained mode. (See Fig. 12,13 for examples from MNIST and
CelebA.) On the other hand, this also showed that the single Gaussian mode
assumption is a good approximation in the initial phase. Since the high-
variance layout information is determined in this initial phase, our theory
illuminates the connection between the initial noise pattern
$\boldsymbol{\mathbf{x}}_{T}$ and the low frequency features in the final
sample $\boldsymbol{\mathbf{x}}_{0}$.

##### Off-manifold dynamics well-predicted by theory.

Eq. 15 indicates that off-manifold components of
$\boldsymbol{\mathbf{x}}_{t}-\alpha_{t}\boldsymbol{\mathbf{\mu}}$ should
shrink at a universal rate that depends only on $\alpha_{t}$. To test this, we
picked a few PCs along which our data sets vary the least, and found that
dynamics in these directions are well-predicted by our formula (Fig. 15). Even
though our theory does not account in detail for mode-splitting, its
predictions for off-manifold directions are still accurate.

##### Prediction errors happen along image manifold.

Next, we examined the deviation between our predictions and the real
trajectory. We found the deviation tends to lie within the subspace spanned by
the top PCs (Fig. 5D, Fig.12E, Fig.13E). One interpretation is that mode-
splitting majorly happens within the general image manifold: after mode-
splitting, the new effective mean and covariance change along the image
manifold, which induces error for the trajectory predicted by using the old
static mean and covariance.

Figure 5: Comparing analytical solution to actual trajectory for a CIFAR-10
diffusion model. A. Actual
$\hat{\boldsymbol{\mathbf{x}}}_{0}(\boldsymbol{\mathbf{x}}_{t})$ and
analytical solution using same initial conditions. B. Collections of actual
and theory-predicted $\boldsymbol{\mathbf{x}}_{0}$. C. Mean squared error
between the actual and analytical solution of trajectory
$\boldsymbol{\mathbf{x}}_{t}$. D. Fraction of squared error as a function of
PC, between the actual and analytical $\boldsymbol{\mathbf{x}}_{0}$.

##  6 Application: Stable Diffusion

In this section, we study Stable Diffusion, a powerful latent diffusion model
that can generate images from natural language text prompts (Rombach et al.,
2022). It is a conditional diffusion model: when we fix a text prompt $\tau$,
the model samples from the distribution $p(\boldsymbol{\mathbf{x}}|\tau)$.
Properties of this high-dimensional distribution like
$\boldsymbol{\mathbf{\mu}}$ and $\boldsymbol{\mathbf{\Sigma}}$ are extremely
hard to characterize (partly due to an inaccessible ground truth), but we can
still show that most of our qualitative predictions hold. We can also use our
theory to guide the search for ‘on-manifold’ directions, which can be used for
editing and creating variations.

Figure 6: Effects of perturbation during reverse diffusion. A. The PCs of the
state difference along the trajectory are on-manifold. B. The projected
outcome $G(\hat{\boldsymbol{\mathbf{x}}}_{0}(\boldsymbol{\mathbf{x}}_{t}))$
for the perturbed trajectory $\tilde{\boldsymbol{\mathbf{x}}_{t}}$ when the
same perturbation is injected at different steps $5,10,15,20,35$. Here the
perturbation is aligned with PC3 of the U-Net output along the original
trajectory. C. The LPIPS distance (Zhang et al., 2018) between the projected
outcome in the original trajectory and perturbed trajectory,
$D(G(\hat{\boldsymbol{\mathbf{x}}}_{0}(\boldsymbol{\mathbf{x}}_{t})),G(\hat{\boldsymbol{\mathbf{x}}}_{0}(\tilde{\boldsymbol{\mathbf{x}}}_{t})))$.

###  6.1 PCs of latent trajectory are on-manifold variations

We observed that Stable Diffusion trajectories are higher-dimensional than
other models’ trajectories (Table 2). What do these additional dimensions
represent? Our theory predicts that the dimensions outside the
$\boldsymbol{\mathbf{x}}_{0},\boldsymbol{\mathbf{x}}_{T}$ plane are lying
within the image manifold spanned by $\boldsymbol{\mathbf{U}}$ (Eq.4.2).
Indeed, when we visualized the PC directions through the decoder, they
appeared to be interpretable variations around the target image (Fig. 6A, 16).
This gave us an efficient way to find local on-manifold vectors using a single
sampling trajectory.

###  6.2 Effects of perturbation on image generation

Ho et al. (2020) use a noisy SDE solver, which implicitly includes
perturbations at every time step; they notice that earlier perturbations tend
to affect higher-level features, while later perturbations tend to affect
lower-level features. By perturbing the latent state at different time steps,
we notice the same phenomenon (Fig. 6A). Specifically, we perturb the latent
state once at time $t^{\prime}$ by an amount $K\delta\boldsymbol{\mathbf{x}}$,
then continue generation. At any time point, the effect on the final image
scales with the strength of perturbation (Fig. 18). As predicted by single
mode theory (Eq. 16), when strength and direction are fixed, the earlier the
perturbation, the larger the effect on the image (Fig. 6A,B). Plainly
speaking, the earlier perturbations have more time to be “amplified” to the
final result. We also noticed that the effect of perturbation is rapidly
reflected by the projected outcome
$G(\hat{\boldsymbol{\mathbf{x}}}_{0}(\boldsymbol{\mathbf{x}}_{t}))$, which
allows for fast testing and synthesis of perturbations (Fig. 6B).

At a fixed time point, not all perturbations are created equal. At the same
scale $K$, Gaussian noise perturbations often have much less perceptible
effect on the final image, while perturbations along trajectory PC directions
often have a significant effect. Quantitatively, the difference between the
perturbed and original trajectory along the perturbed grows through time when
the perturbation is along a PC direction; in contrast, the difference along a
noise direction will decrease over time (Fig.19,20). This observation aligns
with our theory (c.f. Eq. 16), which says off-manifold noise directions will
shrink to zero and signal directions will be amplified to the target standard
deviation $\sqrt{\lambda_{k}}$. Given a high-dimensional image space, it is
highly likely that a randomly chosen perturbation has little overlap with the
signal manifold. Can we design better perturbations?

###  6.3 Perturbation design guided by sampling trajectory

The PCs of a sampling trajectory correspond to meaningful and performant ‘on-
manifold’ perturbations. One nice feature of applying these perturbations is
that on-manifold directions need not be perfectly known, since the ‘noise’
directions will decay to zero by the end. Hence, as long as the perturbation
direction has enough signal direction in it, one tends to obtain good results.

Using this method, we can create manifolds of smoothly connected images,
answering a long-standing question about the existence of latent manifolds in
diffusion models (Fig.17 18). One nice application of the method is to create
local variations of an existing sample without changing the noise seed.

##  7 Discussion

To what extent do our main findings—outline-first and details-later image
generation, low-dimensional trajectories, and increasing commitment to image
elements—hold true for other diffusion model variants and other kinds of
generative models? Simulating reverse SDEs (Eq. 2) instead of ODEs should not
yield many qualitative differences; instead of linear ODEs, one has similarly-
behaved OU processes. The results of Antognini and Sohl-Dickstein (2018) show
OU trajectories are also very low-dimensional. Arguments similar to ones we
have made suggest other kinds of models may be e.g. hierarchical. Xu and Liu
et al.’s Poisson flow generative models (2022) share the qualitative splitting
behavior we described in Section 4.5, but feature smeared-out charge
distributions instead of smeared out Gaussian modes.

Our theory provides a conceptual bridge between diffusion models and GANs. In
the GAN literature, the idea of generating images by modeling progressively
greater levels of resolution is well-established, e.g. by successful
architectures like Progressive Growing GAN and StyleGAN (Karras et al., 2017,
2020). In this paradigm, early layers synthesize the rough layout of the image
from noise, while the last few layers add realistic details. At any layer, the
image can be visualized by a decoding head, with the decoded image gradually
adding details as it approaches the final result. We showed that when looking
at the projected outcome $\hat{\boldsymbol{\mathbf{x}}}_{0}$, diffusion models
have an intriguingly similar generative process: the first few steps map the
noise pattern to a rough outline, and then finer details are gradually added
in. The only significant difference is that GAN layers are replaced by reverse
diffusion sampling steps. This connection may inspire future designs of
diffusion models or GANs.

There may also be a sense in which the latent manifold geometry of GANs and
diffusion models is similar. Wang and Ponce (2021) found that perturbing the
large variance directions of a latent space metric has large and interpretable
effects on image generation in GANs. If we interpret the covariance matrix of
a Gaussian mode as inducing a metric on the latent space of diffusion models,
our perturbation-related observations can be cast in a similar light.

In psychology, the idea of generating mental images hierarchically is related
to the idea of non-commitment: not all details of a scene may be imagined.
Non-commitment is pervasive in human mental imagery (Bigelow et al., 2022). We
found that visualizing endpoint estimates also yields a sort of non-commitment
(e.g. to fine face details) early in reverse diffusion, and it would be
interesting to see to what extent this phenomenon resembles the psychological
one.

## References

  * Anderson (1982) Anderson, B. D.  Reverse-time diffusion equation models.  _Stochastic Processes and their Applications_ , 12(3):313–326, 1982.  ISSN 0304-4149.  doi: https://doi.org/10.1016/0304-4149(82)90051-5.  URL https://www.sciencedirect.com/science/article/pii/0304414982900515.
  * Antognini & Sohl-Dickstein (2018) Antognini, J. M. and Sohl-Dickstein, J.  PCA of high dimensional random walks with comparison to neural network training.  In _NeurIPS_ , pp. 10328–10337, 2018.
  * Ball (1992) Ball, K.  A lower bound for the optimal density of lattice packings.  _International Mathematics Research Notices_ , 1992(10):217–221, 05 1992.  ISSN 1073-7928.  doi: 10.1155/S1073792892000242.  URL https://doi.org/10.1155/S1073792892000242.
  * Bigelow et al. (2022) Bigelow, E. J., McCoy, J., and Ullman, T. D.  Non-commitment in mental imagery, Oct 2022.  URL psyarxiv.com/pn4zd.
  * Goodfellow et al. (2014) Goodfellow, I., Pouget-Abadie, J., Mirza, M., Xu, B., Warde-Farley, D., Ozair, S., Courville, A., and Bengio, Y.  Generative adversarial nets.  In Ghahramani, Z., Welling, M., Cortes, C., Lawrence, N., and Weinberger, K. (eds.), _Advances in Neural Information Processing Systems_ , volume 27. Curran Associates, Inc., 2014.  URL https://proceedings.neurips.cc/paper/2014/file/5ca3e9b122f61f8f06494c97b1afccf3-Paper.pdf.
  * Hertz et al. (2022) Hertz, A., Mokady, R., Tenenbaum, J., Aberman, K., Pritch, Y., and Cohen-Or, D.  Prompt-to-prompt image editing with cross attention control.  _arXiv preprint arXiv:2208.01626_ , 2022.
  * Ho & Salimans (2022) Ho, J. and Salimans, T.  Classifier-free diffusion guidance.  _arXiv preprint arXiv:2207.12598_ , 2022.
  * Ho et al. (2020) Ho, J., Jain, A., and Abbeel, P.  Denoising diffusion probabilistic models.  _Advances in Neural Information Processing Systems_ , 33:6840–6851, 2020.
  * Ho et al. (2022) Ho, J., Salimans, T., Gritsenko, A. A., Chan, W., Norouzi, M., and Fleet, D. J.  Video diffusion models.  In _ICLR Workshop on Deep Generative Models for Highly Structured Data_ , 2022.  URL https://openreview.net/forum?id=BBelR2NdDZ5.
  * Karras et al. (2017) Karras, T., Aila, T., Laine, S., and Lehtinen, J.  Progressive growing of GANs for improved quality, stability, and variation.  _arXiv preprint arXiv:1710.10196_ , 2017.
  * Karras et al. (2020) Karras, T., Laine, S., Aittala, M., Hellsten, J., Lehtinen, J., and Aila, T.  Analyzing and improving the image quality of Stylegan.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , pp. 8110–8119, 2020.
  * Kosslyn & Shwartz (1977) Kosslyn, S. M. and Shwartz, S. P.  A simulation of visual imagery.  _Cognitive Science_ , 1(3):265–295, 1977.  ISSN 0364-0213.  doi: https://doi.org/10.1016/S0364-0213(77)80020-7.  URL https://www.sciencedirect.com/science/article/pii/S0364021377800207.
  * Langlois & Roggman (1990) Langlois, J. H. and Roggman, L. A.  Attractive faces are only average.  _Psychological science_ , 1(2):115–121, 1990\.
  * Liu et al. (2022) Liu, L., Ren, Y., Lin, Z., and Zhao, Z.  Pseudo numerical methods for diffusion models on manifolds.  _arXiv preprint arXiv:2202.09778_ , 2022.
  * Rombach et al. (2022) Rombach, R., Blattmann, A., Lorenz, D., Esser, P., and Ommer, B.  High-resolution image synthesis with latent diffusion models.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pp. 10684–10695, 2022.
  * Ruderman (1994) Ruderman, D. L.  The statistics of natural images.  _Network: computation in neural systems_ , 5(4):517, 1994.
  * Sohl-Dickstein et al. (2015) Sohl-Dickstein, J., Weiss, E., Maheswaranathan, N., and Ganguli, S.  Deep unsupervised learning using nonequilibrium thermodynamics.  In Bach, F. and Blei, D. (eds.), _Proceedings of the 32nd International Conference on Machine Learning_ , volume 37 of _Proceedings of Machine Learning Research_ , pp. 2256–2265, Lille, France, 07–09 Jul 2015\. PMLR.  URL https://proceedings.mlr.press/v37/sohl-dickstein15.html.
  * Song et al. (2020) Song, J., Meng, C., and Ermon, S.  Denoising diffusion implicit models.  _arXiv preprint arXiv:2010.02502_ , 2020.
  * Song & Ermon (2019) Song, Y. and Ermon, S.  Generative modeling by estimating gradients of the data distribution.  In Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buc, F., Fox, E., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems_ , volume 32. Curran Associates, Inc., 2019.  URL https://proceedings.neurips.cc/paper/2019/file/3001ef257407d5a371a96dcd947c7d93-Paper.pdf.
  * Song et al. (2021) Song, Y., Sohl-Dickstein, J., Kingma, D. P., Kumar, A., Ermon, S., and Poole, B.  Score-based generative modeling through stochastic differential equations.  In _International Conference on Learning Representations_ , 2021.  URL https://openreview.net/forum?id=PxTIG12RRHS.
  * Sorscher et al. (2022) Sorscher, B., Ganguli, S., and Sompolinsky, H.  Neural representational geometry underlies few-shot concept learning.  _Proceedings of the National Academy of Sciences_ , 119(43):e2200800119, 2022.  doi: 10.1073/pnas.2200800119.  URL https://www.pnas.org/doi/abs/10.1073/pnas.2200800119.
  * Turk & Pentland (1991) Turk, M. and Pentland, A.  Eigenfaces for recognition.  _Journal of cognitive neuroscience_ , 3(1):71–86, 1991.
  * Wang & Ponce (2021) Wang, B. and Ponce, C. R.  A geometric analysis of deep generative image models and its applications.  In _International Conference on Learning Representations_ , 2021.  URL https://openreview.net/forum?id=GH7QRzUDdXG.
  * Wenliang & Moran (2022) Wenliang, L. K. and Moran, B.  Score-based generative model learn manifold-like structures with constrained mixing.  In _NeurIPS 2022 Workshop on Score-Based Methods_ , 2022.  URL https://openreview.net/forum?id=eSZqaIrDLZR.
  * Xu et al. (2022) Xu, Y., Liu, Z., Tegmark, M., and Jaakkola, T. S.  Poisson flow generative models.  In Oh, A. H., Agarwal, A., Belgrave, D., and Cho, K. (eds.), _Advances in Neural Information Processing Systems_ , 2022.  URL https://openreview.net/forum?id=voV_TRqcWh.
  * Yang et al. (2022) Yang, L., Zhang, Z., Song, Y., Hong, S., Xu, R., Zhao, Y., Shao, Y., Zhang, W., Cui, B., and Yang, M.-H.  Diffusion models: A comprehensive survey of methods and applications.  _arXiv preprint arXiv:2209.00796_ , 2022.
  * Zhang et al. (2018) Zhang, R., Isola, P., Efros, A. A., Shechtman, E., and Wang, O.  The unreasonable effectiveness of deep features as a perceptual metric.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pp. 586–595, 2018.

##  Appendix A Diffusion models used in numerical experiments

Table 1: Diffusion models used for this paper’s numerical experiments.
${\dagger}$: The MNIST diffusion model uses the upsampled $3\times 32\times 32$ RGB pixel space as sample space, while the original MNIST data set consists of $28\times 28$ single channel black and white images. Thus the effective dimensionality of these images is around $784$.  Data set | Hugging Face model_id | Dimensionality | Latents? | Conditional?
---|---|---|---|---
MNIST | dimpo/ddpm-mnist |  3x32x32=3072 ${\dagger}$ | $\times$ | $\times$
CIFAR-10 | google/ddpm-cifar10-32 | 3x32x32=3072 | $\times$ | $\times$
Lsun-Church | google/ddpm-church-256 | 3x256x256 = 196,608 | $\times$ | $\times$
CelebA-HQ | google/ddpm-celebahq-256 | 3x256x256 = 196,608 | $\times$ | $\times$
LAION-2B | runwayml/stable-diffusion-v1-5 | 4x64x64 = 16,384 | $\surd$ | $\surd$

##  Appendix B Supplementary Results

###  B.1 Visualizing image generation progress for other diffusion models

Figure 7: Visualizing the image generation process for the DDPM-CelebA model.
Same layout as Fig.1.

Figure 8: Visualizing the image generation process for the DDPM-Church model.
Same layout as Fig.1.

Figure 9: Visualizing the image generation process for the DDPM-MNIST model.
Same layout as Fig.1.

Figure 10: Visualizing the image generation process for the DDPM-CIFAR-10
model. Same layout as Fig.1.

###  B.2 Trajectory geometry statistics

Table 2: Trajectory geometry statistics for different diffusion models. Residual variance is here defined to be the squared norm of the error vector divided by the squared norm $\|\boldsymbol{\mathbf{x}}_{t}\|^{2}$. Residual variance is computed for three approximations: 1) projecting a trajectory onto the top 2 PCs, 2) projecting a trajectory onto the plane spanned by $\boldsymbol{\mathbf{x}}_{0}$ and $\boldsymbol{\mathbf{x}}_{T}$, 3) approximating the trajectory by a rotation $\alpha_{t}\boldsymbol{\mathbf{x}}_{0}+\sqrt{1-\alpha_{t}^{2}}\boldsymbol{\mathbf{x}}_{T}$. Effective dimensionality is defined as the number of PCs needed to account for 99.9% of the variance. Shown are the effective dimensionalities of the trajectory $\boldsymbol{\mathbf{x}}_{t}$, difference $\boldsymbol{\mathbf{x}}_{t-1}-\boldsymbol{\mathbf{x}}_{t}$, and the U-Net output $\boldsymbol{\mathbf{\epsilon}}_{\theta}(\boldsymbol{\mathbf{x}}_{t})$. All sampling used 51 time steps. |  | Residual Variance | Dim. for 99.9 Var.
---|---|---|---
| Sampler |  | top 2 PC
---
proj.
| $x_{0},x_{T}$
---
proj.
| $x_{0},x_{T}$
---
rotation
$x_{t}$ | $\Delta x_{t}$ | $\epsilon_{\theta}(x_{t})$
ddpm-mnist | DDIM | 0.08% | 0.26% | 1.70% | 2 | 5 | 8
ddpm-cifar10-32 | DDIM | 0.05% | 0.30% | 1.01% | 2 | 4 | 7
ddpm-church-256 | DDIM | 0.05% | 0.31% | 1.09% | 2 | 5 | 8
ddpm-celebahq-256 | DDIM | 0.02% | 0.10% | 1.07% | 2 | 4 | 7
stable-diffusion-v1-5 | PNDM | 0.46% | 0.81% | 2.49% | 5 | 34 | 28

Figure 11: Geometry of diffusion sampling trajectories for 5 different
diffusion models. Left panel: the error fraction of the projection onto the 2D
plane defined by $\boldsymbol{\mathbf{x}}_{0}$ and
$\boldsymbol{\mathbf{x}}_{T}$, and the error fraction under the rotation
approximation. Right panel: the cumulative explained variance of PCs for the
trajectory $\boldsymbol{\mathbf{x}}_{t}$, state differences
$\boldsymbol{\mathbf{x}}_{t-1}-\boldsymbol{\mathbf{x}}_{t}$, and the output
from the U-Net
$\boldsymbol{\mathbf{\epsilon}}_{\theta}(\boldsymbol{\mathbf{x}}_{t})$. 2 PCs
explained almost all variance for the $\boldsymbol{\mathbf{x}}_{t}$
trajectory, while the state difference
$\boldsymbol{\mathbf{x}}_{t-1}-\boldsymbol{\mathbf{x}}_{t}$ and U-Net outputs
are higher dimensional. For quantification see Tab. 2.

###  B.3 Validation of the single mode theory on CIFAR, MNIST, and CelebA

Figure 12: Comparing analytical solution and actual diffusion process for the
DDPM-MNIST model. A. True and predicted endpoint estimate
$\hat{\boldsymbol{\mathbf{x}}}_{0}(\boldsymbol{\mathbf{x}}_{t})$ throughout
reverse diffusion. B. Collection of samples of diffusion-generated images and
the corresponding images predicted by our analytical theory. C. Mean squared
error of the trajectory $\boldsymbol{\mathbf{x}}_{t}$. D. MSE of the endpoint
estimate during diffusion. E. $\boldsymbol{\mathbf{x}}_{0}$ prediction error
along each eigendimension.

Figure 13: Comparing analytical solution and actual diffusion process for the
DDPM-CelebA model. Same layout as Fig.12. Note that in A. the general layout
and shading around the face is consistent between the theory and actual
diffusion trajectory. Note in B, the zoomed-out version of the predicted and
actual sample look highly similar.

Figure 14: Comparing analytical solution and actual diffusion process for the
DDPM-CIFAR-10 model. Same layout as Fig.12.

Figure 15: Trajectory along off-manifold directions well-predicted by theory.
We computed the PC projection of actual trajectories and compared them to the
theoretical prediction given the same initial value; they aligned well. DDPM-
MNIST model.

###  B.4 Principal dimensions of Stable Diffusion sampling trajectories

Figure 16: Visualizing PCs of Stable Diffusion trajectories. We computed the
principal components of the trajectory $\boldsymbol{\mathbf{x}}_{t}$ and the
trajectory difference
$\boldsymbol{\mathbf{x}}_{t-1}-\boldsymbol{\mathbf{x}}_{t}$, and visualized
the scaled version of PC vectors through the decoder. We can see that they
represent an interpretable vector space around the actual sample
$\boldsymbol{\mathbf{x}}_{0}$.

###  B.5 Perturbation experiments

Figure 17: Example perturbation experiment result. Different rows represent
different perturbation times, from top to bottom $5,10,...50$. Different
columns represent different perturbation scales, left to right: negative to
positive, $-20,-15,...15,20$. Perturbations are along PC6 of U-Net output.
Stable Diffusion, PNDM sampler, seed 100. Prompt: “a portrait of an
aristocrat”.

Figure 18: Example perturbation experiment result. Different rows represent
different perturbation times, from top to bottom $5,10,...50$. Different
columns represent different perturbation scales, left to right: negative to
positive, $-20,-15,...15,20$. Perturbations are along PC5 of U-Net output.
Stable Diffusion, PNDM sampler, seed 101. Prompt: “a portrait of an
aristocrat”.

Figure 19: Example perturbation experiment result. Different rows represent
different perturbation times, from top to bottom $5,10,...50$. Different
columns represent different perturbation scales, left to right: negative to
positive, $-20,-15,...15,20$. Perturbations are along a random noise
direction, which is less effective than the previous PC perturbation. Stable
Diffusion, PNDM sampler, seed 101. Prompt: “a portrait of an aristocrat”.

Figure 20: Geometry of trajectory perturbation. A. Quantification of sample
deviation due to a perturbation at time $t^{\prime}$ and strength $K$, heatmap
encodes the perceptual image similarity per LPIPS. Left panel: perturbation
along PC02, Right panel: perturbation along a random vector. The perturbation
along the PC02 is much more effective than the random pattern in affecting the
sample. B. The underlying geometry, projecting the difference between the
perturbed trajectory and the original trajectory onto the unit perturbation
vector. As we predicted the perturbation along signal manifold gets amplified,
while the random perturbation decays over time. (16)

###  B.6 Impact of design choices on the geometry of diffusion models

The field of diffusion generative modeling has greatly advanced since the DDIM
paper (Song et al., 2020). Here, we consider how making slightly different
design choices (e.g. using different samplers) affects the geometry and
dynamics of sampling trajectories.

####  B.6.1 Differential equation solver

Stable Diffusion by default uses the deterministic PNDM sampler, but one can
use other solvers. Because the choice of solver slightly modifies the size and
direction of individual time steps, changing the solver is akin to doing an
early perturbation on the trajectory. Consistent with this expectation, we
observe a variety of effects on image generation with a fixed random seed and
different solvers: sometimes there is effectively no difference (e.g. DDIM and
PNDM), sometimes there is a slight difference, and sometimes there is a major
difference (e.g. PNDM vs LMSDiscrete Solver) (Fig.21). As for trajectory
geometry, usually, DDIM, PNDM, and DPMSolverMultistep create rotation-like 2D
trajectories well-predicted by our theory, while LMSDiscrete and EulerDiscrete
create more linear 1D trajectories.

####  B.6.2 Classifier-free guidance strength

Classifier-free guidance (Ho & Salimans, 2022) has been widely used in
conditional diffusion models as a method to generate samples highly aligned
with the conditional signal (e.g. prompt). We examined the effect of the
strength of classifier-free guidance on the geometry of trajectories. We found
that, generally, a higher guidance scale generates trajectories
$\boldsymbol{\mathbf{x}}_{t}$ and trajectory differences
$\Delta\boldsymbol{\mathbf{x}}_{t}$ with higher dimensionality (Fig.22 top).
Furthermore, when visualizing the top PC vectors, a larger number of
interpretable PC dimensions can be found for a higher guidance scale (Fig. 22
bottom). We also observed that with a smaller guidance value, the trajectory
is usually smooth; with a higher guidance value, it induced strong oscillatory
movement in the trajectory at the early phase, when combined with certain
higher-order schedulers e.g. the default solver PNDM (Liu et al., 2022) (Fig.
22 middle). This effect could be a feature for the sampler to explore the
landscape more. It could also be an artifact, which could be fixed by
modifying the sampler.

Figure 21: Effect of classifier-free guidance (cfg) strength and diffusion
sampler on the sample.

Figure 22: Effect of classifier-free guidance strength and sampler on
trajectory geometry. Top panel: Dimensionality of trajectory
$\boldsymbol{\mathbf{x}}_{t}$ and state difference
$\Delta\boldsymbol{\mathbf{x}}_{t}$, measured by explained variance of PCs.
Higher guidance induces higher dimensionality in $\boldsymbol{\mathbf{x}}_{t}$
and $\Delta\boldsymbol{\mathbf{x}}_{t}$. Middle panel: Trajectory difference
$\Delta\boldsymbol{\mathbf{x}}_{t}$ projected onto the
$\boldsymbol{\mathbf{x}}_{0},\boldsymbol{\mathbf{x}}_{T}$ plane. Higher
guidance induced oscillation in the search trajectory, esp. for PNDM sampler.
Bottom panel: Comparing the top 16 PCs of state difference
$\Delta\boldsymbol{\mathbf{x}}_{t}$, for high guidance (7.5) versus low
guidance (1.0) trajectory from the same noise seed. The samples are shown in
the middle. Trajectories sampled with higher guidance have more ‘on-manifold’
PC dimensions.

##  Appendix C Endpoint estimate trajectory examples

Figure 23: Example endpoint estimate trajectory. CelebA-HQ, DDIM sampler, seed
129.

Figure 24: Example endpoint estimate trajectory. CelebA-HQ, DDIM sampler, seed
152.

Figure 25: Example endpoint estimate trajectory. Stable Diffusion, PNDM
sampler, seed 101. Prompt: “a portrait of an aristocrat”.

Figure 26: Example endpoint estimate trajectory. Stable Diffusion, PNDM
sampler, seed 107. Prompt: “a portrait of an aristocrat”.

##  Appendix D Notation correspondence

Diffusion models usually have forward processes whose conditional
probabilities are

$p(\boldsymbol{\mathbf{x}}_{t}|\boldsymbol{\mathbf{x}}_{0})=\mathcal{N}(A_{t}\boldsymbol{\mathbf{x}}_{0},B_{t}\boldsymbol{\mathbf{I}})$
(20)

for all $t\in[0,T]$. In the limit of small time steps, the transition
probability distribution can be captured by the SDE

$\dot{\boldsymbol{\mathbf{x}}}_{t}=-C_{t}\boldsymbol{\mathbf{x}}_{t}+D_{t}\boldsymbol{\mathbf{\eta}}(t)$
(21)

where $\boldsymbol{\mathbf{\eta}}(t)$ is a vector of independent Gaussian
white noise terms.

Papers discussing these models may use slightly different notation. In the
table below, we briefly indicate how various choices of notation correspond to
one another. To make comparing discrete and continuous models easier, we
assume the time step size is $\Delta t=1$.

Table 3: Comparison of notation for diffusion model parameters. Paper | Citation | $A_{t}$ | $B_{t}$ | $C_{t}$ | $D_{t}$
---|---|---|---|---|---
DDPM | (Ho et al., 2020) | $\sqrt{\bar{\alpha}_{t}}$ | $1-\bar{\alpha}_{t}$ | $1-\sqrt{1-\beta_{t}}$ | $\sqrt{\beta_{t}}$
DDIM | (Song et al., 2020) | $\sqrt{\alpha_{t}}$ | $1-\alpha_{t}$ | $1-\sqrt{\alpha_{t}/\alpha_{t-1}}$ | $\sqrt{1-\alpha_{t}/\alpha_{t-1}}$
Stable Diff. | (Rombach et al., 2022) | $\alpha_{t}$ | $\sigma_{t}^{2}$ | $1-\alpha_{t}/\alpha_{t-1}$ | $\sqrt{\sigma_{t}^{2}-(\alpha_{t}/\alpha_{t-1})^{2}\sigma_{t-1}^{2}}$
VP SDE | (Song et al., 2021) | $\exp\left[-\frac{1}{2}\int_{0}^{t}\beta(s)ds\right]$ | $1-\exp\left[-\int_{0}^{t}\beta(s)ds\right]$ | $\beta(t)/2$ | $\sqrt{\beta(t)}$
Ours |  | $\alpha_{t}$ | $\sigma_{t}^{2}$ | $\beta(t)$ | $g(t)$

In the huggingface diffusers library implementation of reverse diffusion, the
function alphas_cumprod corresponds to our $\alpha_{t}^{2}$.

##  Appendix E High-dimensional Gaussians are thin shells

Understanding the counterintuitive properties of Gaussians in high dimensions
is important for understanding reverse diffusion dynamics in a landscape of
high-dimensional Gaussian modes. Although this counterintuitive behavior is
well-known (see e.g. (Song et al., 2021; Sorscher et al., 2022)), we present
the relevant details here for the reader’s convenience.

Consider an isotropic Gaussian centered at the origin in $D$ dimensions, where
$D\gg 1$. Its probability density function is

$p(\boldsymbol{\mathbf{x}})=\frac{1}{(2\pi\sigma^{2})^{D/2}}\exp\left\\{-\frac{\boldsymbol{\mathbf{x}}^{2}}{2\sigma^{2}}\right\\}\
.$ (22)

In spherical coordinates with $r:=\sqrt{\boldsymbol{\mathbf{x}}^{2}}$, the
probability of measuring a value near $r$ is

$p(r)=\frac{2}{\Gamma(D/2)(2\sigma^{2})^{D/2}}\
r^{D-1}\exp\left\\{-\frac{r^{2}}{2\sigma^{2}}\right\\}$ (23)

where the factor of $r^{D-1}$ comes from the Jacobian determinant associated
with the change of coordinates. The peak of this distribution occurs at
$r_{*}=\sqrt{D-1}\ \sigma$ and its mean occurs at $\langle
r\rangle=\sqrt{D}\sigma$. Its variance is $\text{var}(r)=2\sigma^{2}$, which
suggests that almost all of the probability mass is contained in a thin shell,
i.e.

$(\sqrt{D}-\sqrt{2})\sigma\leq r\leq(\sqrt{D}+\sqrt{2})\sigma\ .$ (24)

As $D$ gets large, fluctuations about this region become small compared to the
typical radius, and essentially all probability lies on the shell. By
appropriately modifying this argument, it is also true for a non-isotropic
Gaussian. In this more general case, the peak of $p(r)$ occurs on
$(D-1)$-dimensional ellipsoid defined by

$\boldsymbol{\mathbf{x}}^{T}\boldsymbol{\mathbf{\Sigma}}^{-1}\boldsymbol{\mathbf{x}}=D-1\
.$ (25)

Most probability is concentrated on the surface of this ellipsoid, i.e.

$\sqrt{D}-\sqrt{2}\leq\sqrt{\boldsymbol{\mathbf{x}}^{T}\boldsymbol{\mathbf{\Sigma}}^{-1}\boldsymbol{\mathbf{x}}}\leq\sqrt{D}+\sqrt{2}\
.$ (26)

The fact that $D$-dimensional Gaussians are essentially $(D-1)$-dimensional
ellipsoids has interesting consequences for mean-reversion dynamics (like what
we study in this paper) in high dimensions. As Antognini and Sohl-Dickstein
(2018) point out, OU processes in high dimensions almost never reach the mode,
or even get close to it, despite the defining characteristic of OU processes
being mean reversion. When particles reach the ellipsoid described above, the
driving force toward the mode is balanced by the noise in the process, causing
them to be overwhelmingly likely to remain on the ellipsoid.

##  Appendix F Theory of the projected outcome
$\hat{\boldsymbol{\mathbf{x}}}_{0}$

In this section, we study the endpoint estimate
$\hat{\boldsymbol{\mathbf{x}}}_{0}$ in the case of a Gaussian distribution. We
use this theory in the main text to understand the short time dynamics of
reverse diffusion.

###  F.1 Projected outcome for a Gaussian image distribution

If the image distribution is a Gaussian, i.e.
$p(\boldsymbol{\mathbf{x}})=\mathcal{N}(\boldsymbol{\mathbf{\mu}},\boldsymbol{\mathbf{\Sigma}})$,
then we have the explicit score function $\nabla_{\boldsymbol{\mathbf{x}}}\log
p(\boldsymbol{\mathbf{x}})=\boldsymbol{\mathbf{\Sigma}}^{-1}(\boldsymbol{\mathbf{\mu}}-\boldsymbol{\mathbf{x}})$.
The forward diffusion process yields the conditional distribution
$p(\boldsymbol{\mathbf{x}}_{t}\mid\boldsymbol{\mathbf{x}}_{0})=\mathcal{N}(\alpha_{t}\boldsymbol{\mathbf{x}}_{0},\sigma_{t}^{2}\boldsymbol{\mathbf{I}})$,
whose marginal distribution is
$p(\boldsymbol{\mathbf{x}}_{t})=\mathcal{N}(\alpha_{t}\boldsymbol{\mathbf{\mu}},\alpha_{t}^{2}\boldsymbol{\mathbf{\Sigma}}+\sigma_{t}^{2}\boldsymbol{\mathbf{I}})$.
When time $t$ increases, $\alpha_{t}$ decreases from 1 to 0, while
$\sigma_{t}$ increases.

The time-dependent score is

$\boldsymbol{\mathbf{s}}(\boldsymbol{\mathbf{x}},t)=\nabla_{\boldsymbol{\mathbf{x}}}\log
p(\boldsymbol{\mathbf{x}}_{t})=(\sigma_{t}^{2}\boldsymbol{\mathbf{I}}+\alpha_{t}^{2}\boldsymbol{\mathbf{\Sigma}})^{-1}(\alpha_{t}\boldsymbol{\mathbf{\mu}}-\boldsymbol{\mathbf{x}})\
.$

Note that in the most common implementation of diffusion models, people use a
neural network to parameterize the scaled negative score function
$\boldsymbol{\mathbf{\epsilon}}_{\theta}(\boldsymbol{\mathbf{x}},t)=-\sigma_{t}\boldsymbol{\mathbf{s}}(\boldsymbol{\mathbf{x}},t)$.
Then we can write the analytical form of the estimate
$\hat{\boldsymbol{\mathbf{x}}}_{0}(\boldsymbol{\mathbf{x}}_{t})$ as

$\displaystyle\hat{\boldsymbol{\mathbf{x}}}_{0}(\boldsymbol{\mathbf{x}}_{t})$
$\displaystyle=\frac{\boldsymbol{\mathbf{x}}_{t}-\sigma_{t}\boldsymbol{\mathbf{\epsilon}}_{\theta}(\boldsymbol{\mathbf{x}}_{t},t)}{\alpha_{t}}$
(27)
$\displaystyle=\frac{\boldsymbol{\mathbf{x}}_{t}+\sigma_{t}^{2}\boldsymbol{\mathbf{s}}(\boldsymbol{\mathbf{x}}_{t},t)}{\alpha_{t}}$
(28)
$\displaystyle=\frac{\boldsymbol{\mathbf{x}}_{t}+\sigma_{t}^{2}(\sigma_{t}^{2}\boldsymbol{\mathbf{I}}+\alpha_{t}^{2}\boldsymbol{\mathbf{\Sigma}})^{-1}(\alpha_{t}\boldsymbol{\mathbf{\mu}}-\boldsymbol{\mathbf{x}}_{t})}{\alpha_{t}}\
.$ (29)

Let’s examine the difference between
$\hat{\boldsymbol{\mathbf{x}}}_{0}(\boldsymbol{\mathbf{x}}_{t})$ and the mean
$\boldsymbol{\mathbf{\mu}}$. We find that

$\displaystyle\hat{\boldsymbol{\mathbf{x}}}_{0}(\boldsymbol{\mathbf{x}}_{t})-\boldsymbol{\mathbf{\mu}}$
$\displaystyle=\frac{1}{\alpha_{t}}(\boldsymbol{\mathbf{x}}_{t}-\alpha_{t}\boldsymbol{\mathbf{\mu}}+\sigma_{t}^{2}(\sigma_{t}^{2}\boldsymbol{\mathbf{I}}+\alpha_{t}^{2}\boldsymbol{\mathbf{\Sigma}})^{-1}(\alpha_{t}\boldsymbol{\mathbf{\mu}}-\boldsymbol{\mathbf{x}}_{t}))$
(30)
$\displaystyle=\frac{\sigma_{t}^{2}}{\alpha_{t}}(\frac{1}{\sigma_{t}^{2}}\boldsymbol{\mathbf{I}}-(\sigma_{t}^{2}\boldsymbol{\mathbf{I}}+\alpha_{t}^{2}\boldsymbol{\mathbf{\Sigma}})^{-1})(\boldsymbol{\mathbf{x}}_{t}-\alpha_{t}\boldsymbol{\mathbf{\mu}})\
.$

What can we learn from this formula? We explore several consequences of it
below.

##### The projected outcome is always exact for a delta function.

If the initial distribution $p(\boldsymbol{\mathbf{x}}_{0})$ is a delta
function, $\boldsymbol{\mathbf{\Sigma}}\to 0$. Then at any time $t$,
$p(\boldsymbol{\mathbf{x}}_{t})=\mathcal{N}(\alpha_{t}\boldsymbol{\mathbf{\mu}},\sigma_{t}^{2}\boldsymbol{\mathbf{I}})$,
and the projected outcome is
$\hat{\boldsymbol{\mathbf{x}}}_{0}(\boldsymbol{\mathbf{x}}_{t})-\boldsymbol{\mathbf{\mu}}\equiv
0,\forall\boldsymbol{\mathbf{x}}_{t},\forall t$. Hence, the projected outcome
is always exact, regardless of the position of $\boldsymbol{\mathbf{x}}_{t}$
and the time $t$. This can be regarded as one justification for this statistic
$\hat{\boldsymbol{\mathbf{x}}}_{0}(\boldsymbol{\mathbf{x}}_{t})$: it is exact
and invariant in the isotropic score field created by a point distribution.

This point may be relevant for understanding the very end of reverse diffusion
dynamics, when $\boldsymbol{\mathbf{x}}_{t}$ is likely to live within a score
field created by a single point. At such a time,
$\hat{\boldsymbol{\mathbf{x}}}_{0}(\boldsymbol{\mathbf{x}}_{t})$ scarcely
changes.

##### Projected outcome for an isotropic Gaussian.

When $\boldsymbol{\mathbf{\Sigma}}=\bar{\sigma}^{2}\boldsymbol{\mathbf{I}}$,
Eq. 30 takes a particularly simple form:

$\displaystyle\hat{\boldsymbol{\mathbf{x}}}_{0}(\boldsymbol{\mathbf{x}}_{t})-\boldsymbol{\mathbf{\mu}}$
$\displaystyle=\frac{1}{\alpha_{t}}\frac{\alpha_{t}^{2}\bar{\sigma}^{2}}{\sigma_{t}^{2}+\alpha_{t}^{2}\bar{\sigma}^{2}}(\boldsymbol{\mathbf{x}}_{t}-\alpha_{t}\boldsymbol{\mathbf{\mu}})\
.$ (31)

###  F.2 Low rank $\Sigma$

Images are commonly recognized as residing on a low-dimensional manifold. To
study analytically what this might imply for diffusion models that generate
images, we consider the specific case of a Gaussian whose covariance matrix
$\boldsymbol{\mathbf{\Sigma}}$ is degenerate and low rank. Let
$\boldsymbol{\mathbf{\Sigma}}=\boldsymbol{\mathbf{U}}\boldsymbol{\mathbf{\Lambda}}\boldsymbol{\mathbf{U}}^{T}$
be the eigendecomposition or compact SVD of the covariance matrix, where
$\boldsymbol{\mathbf{U}}$ is a $D\times r$ semi-orthogonal matrix whose
columns are normalized (i.e.
$\boldsymbol{\mathbf{U}}^{T}\boldsymbol{\mathbf{U}}=\boldsymbol{\mathbf{I}}_{r}$),
and $\boldsymbol{\mathbf{\Lambda}}$ is the $r\times r$ diagonal eigenvalue
matrix. Denote the $k$th column of $\boldsymbol{\mathbf{U}}$ by
$\boldsymbol{\mathbf{u}}_{k}$ and the $k$th diagonal element of
$\boldsymbol{\mathbf{\Lambda}}$ by $\lambda_{k}$.

We can use the Woodbury matrix inversion identity to write

$\displaystyle(\sigma_{t}^{2}\boldsymbol{\mathbf{I}}+\alpha_{t}^{2}\boldsymbol{\mathbf{\Sigma}})^{-1}$
$\displaystyle=(\sigma_{t}^{2}\boldsymbol{\mathbf{I}}+\alpha_{t}^{2}\boldsymbol{\mathbf{U}}\boldsymbol{\mathbf{\Lambda}}\boldsymbol{\mathbf{U}}^{T})^{-1}$
(32)
$\displaystyle=\frac{1}{\sigma_{t}^{2}}\boldsymbol{\mathbf{I}}-\frac{1}{\sigma_{t}^{4}}\boldsymbol{\mathbf{U}}(\frac{1}{\alpha_{t}^{2}}\boldsymbol{\mathbf{\Lambda}}^{-1}+\frac{1}{\sigma_{t}^{2}}\boldsymbol{\mathbf{U}}^{T}\boldsymbol{\mathbf{U}})^{-1}\boldsymbol{\mathbf{U}}^{T}$
(33)
$\displaystyle=\frac{1}{\sigma_{t}^{2}}\boldsymbol{\mathbf{I}}-\frac{1}{\sigma_{t}^{4}}\boldsymbol{\mathbf{U}}(\frac{1}{\alpha_{t}^{2}}\boldsymbol{\mathbf{\Lambda}}^{-1}+\frac{1}{\sigma_{t}^{2}}\boldsymbol{\mathbf{I}}_{r})^{-1}\boldsymbol{\mathbf{U}}^{T}\
.$ (34)

Define the diagonal matrix

$\displaystyle\tilde{\boldsymbol{\mathbf{\Lambda}}}_{t}$
$\displaystyle=\frac{1}{\sigma_{t}^{2}}(\frac{1}{\alpha_{t}^{2}}\boldsymbol{\mathbf{\Lambda}}^{-1}+\frac{1}{\sigma_{t}^{2}}\boldsymbol{\mathbf{I}}_{r})^{-1}$
(35)
$\displaystyle=\frac{1}{\sigma_{t}^{2}}\text{diag}[(\alpha_{t}^{2}\lambda_{k})^{-1}+(\sigma_{t}^{2})^{-1})^{-1}]$
(36)
$\displaystyle=\text{diag}\left[\frac{\alpha_{t}^{2}\lambda_{k}}{\alpha_{t}^{2}\lambda_{k}+\sigma_{t}^{2}}\right]\
.$ (37)

The diagonal elements of this matrix have the same form as the function

$\displaystyle\phi(\lambda,t):=\frac{\alpha_{t}^{2}\lambda}{\alpha_{t}^{2}\lambda+\sigma_{t}^{2}}\
.$ (38)

To improve our intuition about the temporal behavior of these diagonal
elements, we plotted this function (Fig.27).

We can use the matrix $\tilde{\boldsymbol{\mathbf{\Lambda}}}_{t}$ to write
$(\sigma_{t}^{2}\boldsymbol{\mathbf{I}}+\alpha_{t}^{2}\boldsymbol{\mathbf{\Sigma}})^{-1}$
as

$\displaystyle(\sigma_{t}^{2}\boldsymbol{\mathbf{I}}+\alpha_{t}^{2}\boldsymbol{\mathbf{\Sigma}})^{-1}$
$\displaystyle=\frac{1}{\sigma_{t}^{2}}(\boldsymbol{\mathbf{I}}-\boldsymbol{\mathbf{U}}\tilde{\boldsymbol{\mathbf{\Lambda}}}_{t}\boldsymbol{\mathbf{U}}^{T})\
.$ (39)

Using this result, we can write the endpoint estimate as

$\displaystyle\hat{\boldsymbol{\mathbf{x}}}_{0}(\boldsymbol{\mathbf{x}}_{t})-\boldsymbol{\mathbf{\mu}}$
$\displaystyle=\frac{\sigma_{t}^{2}}{\alpha_{t}}(\frac{1}{\sigma_{t}^{2}}\boldsymbol{\mathbf{I}}-(\sigma_{t}^{2}\boldsymbol{\mathbf{I}}+\alpha_{t}^{2}\boldsymbol{\mathbf{\Sigma}})^{-1})(\boldsymbol{\mathbf{x}}_{t}-\alpha_{t}\boldsymbol{\mathbf{\mu}})$
(40)
$\displaystyle=\frac{1}{\alpha_{t}}\boldsymbol{\mathbf{U}}\tilde{\boldsymbol{\mathbf{\Lambda}}}_{t}\boldsymbol{\mathbf{U}}^{T}(\boldsymbol{\mathbf{x}}_{t}-\alpha_{t}\boldsymbol{\mathbf{\mu}})\
.$ (41)

This formula has a series of interesting implications.

Figure 27: The function $\phi(\lambda,t)$ with different $\lambda$ values,
using a common $\alpha_{t}$ schedule.

##### Projected outcome stays on image manifold.

Note that the deviation of the endpoint estimate from the distribution mean
$\hat{\boldsymbol{\mathbf{x}}}_{0}(\boldsymbol{\mathbf{x}}_{t})-\boldsymbol{\mathbf{\mu}}$
always remains in the subspace spanned by the columns of
$\boldsymbol{\mathbf{U}}$, i.e.
$\hat{\boldsymbol{\mathbf{x}}}_{0}(\boldsymbol{\mathbf{x}}_{t})-\boldsymbol{\mathbf{\mu}}\in\text{span}(\boldsymbol{\mathbf{U}})$.
This an interesting result: if the data distribution is a low-dimensional
manifold (e.g. image manifold), the estimate
$\hat{\boldsymbol{\mathbf{x}}}_{0}$ will not deviate from this manifold. Even
if the projected outcome does not exactly reflect the true outcome, it will
not make errors out of the image manifold, i.e. orthogonal to high variance
directions. Visually, such directions would correspond to random noise, and
hence perturbations in those directions would be ‘nonsense’ perturbations to
images.

##### Projected outcome starts around the center of the distribution.

At the start of the reverse diffusion, $t\approx T$, $\alpha_{t}\approx 0$,
and $\sigma_{t}^{2}\gg\alpha_{t}^{2}$. This means
$\tilde{\boldsymbol{\mathbf{\Lambda}}}_{t}\approx\boldsymbol{\mathbf{0}}$ and
$\hat{\boldsymbol{\mathbf{x}}}_{0}(\boldsymbol{\mathbf{x}}_{t})-\boldsymbol{\mathbf{\mu}}\approx\boldsymbol{\mathbf{0}}$.
Hence, the projected outcome $\hat{\boldsymbol{\mathbf{x}}}_{0}$ initially
corresponds to the center of the distribution. For multi-class unconditional
generation, this center will be relatively class-ambiguous. For the
unconditional generation of faces, we observed that
$\hat{\boldsymbol{\mathbf{x}}}_{0}(\boldsymbol{\mathbf{x}}_{T})$ points to a
generic face, which is close to the ‘average face’.

##### Projected outcome ends in the real outcome.

At the end of reverse diffusion, $t=0$, $\sigma_{t}=0,$ and $\alpha_{t}=1$.
This means
$\hat{\boldsymbol{\mathbf{x}}}_{0}(\boldsymbol{\mathbf{x}}_{0})=\boldsymbol{\mathbf{x}}_{0}$,
i.e. the projected outcome corresponds to the real result of reverse
diffusion.

##### Features emerge in descending order of their variance.

Consider the projection of this vector
$\hat{\boldsymbol{\mathbf{x}}}_{0}(\boldsymbol{\mathbf{x}}_{t})-\boldsymbol{\mathbf{\mu}}$
on an eigenvector $\boldsymbol{\mathbf{u}}_{k}$ of
$\boldsymbol{\mathbf{\Sigma}}$:

$\displaystyle\boldsymbol{\mathbf{u}}_{k}^{T}(\hat{\boldsymbol{\mathbf{x}}}_{0}(\boldsymbol{\mathbf{x}}_{t})-\boldsymbol{\mathbf{\mu}})=\frac{1}{\alpha_{t}}\frac{\alpha_{t}^{2}\lambda_{k}}{\alpha_{t}^{2}\lambda_{k}+\sigma_{t}^{2}}\boldsymbol{\mathbf{u}}_{k}^{T}(\boldsymbol{\mathbf{x}}_{t}-\alpha_{t}\boldsymbol{\mathbf{\mu}})\
.$ (42)

When $\sigma_{t}^{2}\gg\alpha_{t}^{2}\lambda_{k}$, i.e. when the noise scale
is much larger than the signal scale, we have

$\boldsymbol{\mathbf{u}}_{k}^{T}(\hat{\boldsymbol{\mathbf{x}}}_{0}(\boldsymbol{\mathbf{x}}_{t})-\boldsymbol{\mathbf{\mu}})\ll\frac{\boldsymbol{\mathbf{u}}_{k}^{T}(\boldsymbol{\mathbf{x}}_{t}-\alpha_{t}\boldsymbol{\mathbf{\mu}})}{\alpha_{t}}\
,$ (43)

so the projected outcome is approximately $\boldsymbol{\mathbf{0}}$. At the
other extreme, when $\sigma_{t}^{2}\ll\alpha_{t}^{2}\lambda_{k}$, the signal
variance is much bigger than the noise variance, and we have

$\boldsymbol{\mathbf{u}}_{k}^{T}(\hat{\boldsymbol{\mathbf{x}}}_{0}(\boldsymbol{\mathbf{x}}_{t})-\boldsymbol{\mathbf{\mu}})\approx\frac{\boldsymbol{\mathbf{u}}_{k}^{T}(\boldsymbol{\mathbf{x}}_{t}-\alpha_{t}\boldsymbol{\mathbf{\mu}})}{\alpha_{t}}\
.$ (44)

To interpret this, note that if we regard $\boldsymbol{\mathbf{u}}_{k}$ as a
feature direction, then $\lambda_{k}$ is the variance along this feature
direction. The above equation tells us that this feature stays around the mean
value of the distribution when the noise variance is much larger than the
scaled variance of this feature. In plain language, when the signal scale
along a certain dimension is less than the noise scale, the projected outcome
along that dimension will remain undetermined.

Empirically, people have found that natural image spaces have spectra close to
$1/f$ (Ruderman, 1994), which means more image variance exists in low-
frequency features than in high-frequency features. Thus, we expect that in
the generating process, the projected outcome
$\hat{\boldsymbol{\mathbf{x}}}_{0}(x_{t})$ will specify the low-frequency
layout first, and then gradually add the high-frequency details.

##### Invariance of projected outcome to off-manifold perturbations.

Consider a perturbation of $\boldsymbol{\mathbf{x}}_{t}$ at time $t$ by
$\delta\boldsymbol{\mathbf{x}}$. The projected outcome will change by an
amount

$\displaystyle\hat{\boldsymbol{\mathbf{x}}}_{0}(\boldsymbol{\mathbf{x}}_{t}+\delta\boldsymbol{\mathbf{x}})-\hat{\boldsymbol{\mathbf{x}}}_{0}(\boldsymbol{\mathbf{x}}_{t})$
$\displaystyle=\frac{1}{\alpha_{t}}\boldsymbol{\mathbf{U}}\tilde{\boldsymbol{\mathbf{\Lambda}}}_{t}\boldsymbol{\mathbf{U}}^{T}\delta\boldsymbol{\mathbf{x}}\
.$ (45)

Notice that if $\boldsymbol{\mathbf{x}}_{t}$ is perturbed in a direction
orthogonal to the manifold spanned by the columns of
$\boldsymbol{\mathbf{U}}$, the projected outcome
$\hat{\boldsymbol{\mathbf{x}}}_{0}$ will not change. Thus, perturbation in a
random non-signal direction will not change the projected outcome of reverse
diffusion (and it will also, as we will see, not affect the real result
either). In contrast, if the perturbation is aligned with the image manifold,
the perturbation will affect the projected outcome.

##### Effect of a perturbation decreases over time.

As we can see from the formula above, the effect of the same perturbation
changes over time. After $\alpha_{t}^{2}\lambda_{k}\gg\sigma_{t}^{2}$, the
term
$\frac{\alpha_{t}^{2}\lambda_{k}}{\alpha_{t}^{2}\lambda_{k}+\sigma_{t}^{2}}$
is approximately $1$, so
$\frac{1}{\alpha_{t}}\frac{\alpha_{t}^{2}\lambda_{k}}{\alpha_{t}^{2}\lambda_{k}+\sigma_{t}^{2}}$
will decrease over time. We predict on this basis that after certain features
emerge, perturbations along those dimensions will have a limited effect.

##  Appendix G Theory of state space dynamics $\boldsymbol{\mathbf{x}}_{t}$

In this section, we study the dynamics of the latent state
$\boldsymbol{\mathbf{x}}_{t}$ and endpoint estimate
$\hat{\boldsymbol{\mathbf{x}}}_{0}$. These results are used in Section 4 of
the main text.

###  G.1 Alternative form of state dynamics

Using the form of the projected outcome (Eq. 30), we can write the probability
flow ODE (Eq. 6) as

$\displaystyle\dot{\boldsymbol{\mathbf{x}}}$
$\displaystyle=-\beta(t)\boldsymbol{\mathbf{x}}-\frac{1}{2}g^{2}(t)\nabla_{\boldsymbol{\mathbf{x}}}\log
p(\boldsymbol{\mathbf{x}},t)$ (46)
$\displaystyle=-\beta(t)\boldsymbol{\mathbf{x}}-\frac{1}{2}\frac{g^{2}(t)}{\sigma_{t}^{2}}(\alpha_{t}\hat{\boldsymbol{\mathbf{x}}}_{0}(\boldsymbol{\mathbf{x}})-\boldsymbol{\mathbf{x}})\
.$ (47)

Notice that

$\displaystyle\alpha(t)=\exp(-\int_{0}^{t}\beta(\tau)d\tau)\hskip
21.68121pt\frac{d}{dt}\alpha(t)=-\beta(t)\alpha(t)\ ,$ (48)

which allows us to write

$\displaystyle\frac{d}{dt}(\frac{\boldsymbol{\mathbf{x}}_{t}}{\alpha_{t}})$
$\displaystyle=\frac{\dot{\boldsymbol{\mathbf{x}}}_{t}}{\alpha_{t}}-\frac{\dot{\alpha}_{t}\boldsymbol{\mathbf{x}}_{t}}{\alpha_{t}^{2}}$
(49)
$\displaystyle=\frac{\dot{\boldsymbol{\mathbf{x}}}_{t}}{\alpha_{t}}+\frac{\beta_{t}\boldsymbol{\mathbf{x}}_{t}}{\alpha_{t}}$
(50)
$\displaystyle\frac{d}{dt}(\frac{\boldsymbol{\mathbf{x}}_{t}}{\alpha_{t}})$
$\displaystyle=-\frac{1}{2}\frac{g^{2}(t)}{\sigma_{t}^{2}}(\hat{\boldsymbol{\mathbf{x}}}_{0}(\boldsymbol{\mathbf{x}}_{t})-\frac{\boldsymbol{\mathbf{x}}_{t}}{\alpha_{t}})\
.$ (51)

From this, we can see that the quantity
$\boldsymbol{\mathbf{x}}_{t}/\alpha_{t}$ is being isotropically attracted
towards the moving target
$\hat{\boldsymbol{\mathbf{x}}}_{0}(\boldsymbol{\mathbf{x}}_{t})$ at a rate
determined by $\frac{1}{2}\frac{g^{2}(t)}{\sigma_{t}^{2}}$.

###  G.2 Analytical solution of reverse diffusion dynamics for a Gaussian
distribution

For a multivariate Gaussian distribution, the exact ODE governing the reverse
diffusion sampling process is

$\dot{\boldsymbol{\mathbf{x}}}=-\beta(t)\boldsymbol{\mathbf{x}}-\frac{1}{2}g^{2}(t)(\sigma_{t}^{2}\boldsymbol{\mathbf{I}}+\alpha_{t}^{2}\boldsymbol{\mathbf{\Sigma}})^{-1}(\alpha_{t}\boldsymbol{\mathbf{\mu}}-\boldsymbol{\mathbf{x}})\
.$ (52)

Using the eigendecomposition of $\Sigma$ described in the previous section,

$\displaystyle\dot{\boldsymbol{\mathbf{x}}}$
$\displaystyle=-\beta(t)\boldsymbol{\mathbf{x}}-\frac{1}{2}g^{2}(t)\frac{1}{\sigma_{t}^{2}}(\boldsymbol{\mathbf{I}}-\boldsymbol{\mathbf{U}}\tilde{\boldsymbol{\mathbf{\Lambda}}}_{t}\boldsymbol{\mathbf{U}}^{T})(\alpha_{t}\boldsymbol{\mathbf{\mu}}-\boldsymbol{\mathbf{x}})$
(53)

where we recall that $\tilde{\boldsymbol{\mathbf{\Lambda}}}$ is defined to be
the time-dependent diagonal matrix

$\displaystyle\tilde{\boldsymbol{\mathbf{\Lambda}}}_{t}$
$\displaystyle=\text{diag}\left[\frac{\alpha_{t}^{2}\lambda_{k}}{\alpha_{t}^{2}\lambda_{k}+\sigma_{t}^{2}}\right]\
.$ (54)

Consider the dynamics of the quantity
$\boldsymbol{\mathbf{x}}_{t}-\alpha_{t}\boldsymbol{\mathbf{\mu}}$. Using the
relationship between $\beta_{t}$ and $\alpha_{t}$, we have

$\displaystyle\frac{d}{dt}(\boldsymbol{\mathbf{x}}_{t}-\alpha_{t}\boldsymbol{\mathbf{\mu}})$
$\displaystyle=\dot{\boldsymbol{\mathbf{x}}}_{t}-\boldsymbol{\mathbf{\mu}}\dot{\alpha}_{t}$
(55)
$\displaystyle=\dot{\boldsymbol{\mathbf{x}}}_{t}+\beta_{t}\alpha_{t}\boldsymbol{\mathbf{\mu}}$
(56)
$\displaystyle=\beta_{t}(\alpha_{t}\boldsymbol{\mathbf{\mu}}-\boldsymbol{\mathbf{x}})-\frac{1}{2}g^{2}(t)\frac{1}{\sigma_{t}^{2}}(\boldsymbol{\mathbf{I}}-\boldsymbol{\mathbf{U}}\tilde{\boldsymbol{\mathbf{\Lambda}}}_{t}\boldsymbol{\mathbf{U}}^{T})(\alpha_{t}\boldsymbol{\mathbf{\mu}}-\boldsymbol{\mathbf{x}})$
(57)
$\displaystyle=(\frac{1}{2}g^{2}(t)\frac{1}{\sigma_{t}^{2}}(\boldsymbol{\mathbf{I}}-\boldsymbol{\mathbf{U}}\tilde{\boldsymbol{\mathbf{\Lambda}}}_{t}\boldsymbol{\mathbf{U}}^{T})-\beta_{t}\boldsymbol{\mathbf{I}})(\boldsymbol{\mathbf{x}}-\alpha_{t}\boldsymbol{\mathbf{\mu}})\
.$ (58)

If we assume that the forward process is a variance-preserving SDE, then
$\beta_{t}=\frac{1}{2}g^{2}(t)$, which implies
$\alpha_{t}^{2}=1-\sigma_{t}^{2}$. Using this, we obtain

$\displaystyle\frac{d}{dt}(\boldsymbol{\mathbf{x}}_{t}-\alpha_{t}\boldsymbol{\mathbf{\mu}})$
$\displaystyle=\beta_{t}(\frac{1}{\sigma_{t}^{2}}(\boldsymbol{\mathbf{I}}-\boldsymbol{\mathbf{U}}\tilde{\boldsymbol{\mathbf{\Lambda}}}_{t}\boldsymbol{\mathbf{U}}^{T})-\boldsymbol{\mathbf{I}})(\boldsymbol{\mathbf{x}}-\alpha_{t}\boldsymbol{\mathbf{\mu}})$
(59)
$\displaystyle=\beta_{t}((\frac{1}{\sigma_{t}^{2}}-1)\boldsymbol{\mathbf{I}}-\frac{1}{\sigma_{t}^{2}}\boldsymbol{\mathbf{U}}\tilde{\boldsymbol{\mathbf{\Lambda}}}_{t}\boldsymbol{\mathbf{U}}^{T}))(\boldsymbol{\mathbf{x}}-\alpha_{t}\boldsymbol{\mathbf{\mu}})\
.$ (60)

Define the variable
$\boldsymbol{\mathbf{y}}_{t}:=\boldsymbol{\mathbf{x}}_{t}-\alpha_{t}\boldsymbol{\mathbf{\mu}}$.
We have just shown that its dynamics are fairly ‘nice’, in the sense that the
above equation is well-behaved separable linear ODE. As we are about to show,
it is exactly solvable.

Write $\boldsymbol{\mathbf{y}}_{t}$ in terms of the orthonormal columns of
$\boldsymbol{\mathbf{U}}$ and a component that lies entirely in the orthogonal
space $\boldsymbol{\mathbf{U}}^{\perp}$:

$\boldsymbol{\mathbf{y}}_{t}=\boldsymbol{\mathbf{y}}^{\perp}(t)+\sum_{k=1}^{r}c_{k}(t)\boldsymbol{\mathbf{u}}_{k}\
,\ \boldsymbol{\mathbf{y}}^{\perp}(t)\in\boldsymbol{\mathbf{U}}^{\perp}\ .$
(61)

The dynamics of the coefficient $c_{k}(t)$ attached to the eigenvector
$\boldsymbol{\mathbf{u}}_{k}$ are

$\displaystyle\dot{c}_{k}(t)=\frac{d}{dt}(\boldsymbol{\mathbf{u}}_{k}^{T}\boldsymbol{\mathbf{y}}_{t})$
$\displaystyle=\beta_{t}((\frac{1}{\sigma_{t}^{2}}-1)-\frac{1}{\sigma_{t}^{2}}\frac{\alpha_{t}^{2}\lambda_{k}}{\alpha_{t}^{2}\lambda_{k}+\sigma_{t}^{2}})(\boldsymbol{\mathbf{u}}_{k}^{T}\boldsymbol{\mathbf{y}}_{t})$
(62)
$\displaystyle=\frac{\beta_{t}}{\sigma_{t}^{2}}(1-\sigma_{t}^{2}-\frac{\alpha_{t}^{2}\lambda_{k}}{\alpha_{t}^{2}\lambda_{k}+\sigma_{t}^{2}})c_{k}(t)$
(63)
$\displaystyle=\frac{\beta_{t}\alpha_{t}^{2}}{\sigma_{t}^{2}}(1-\frac{\lambda_{k}}{\alpha_{t}^{2}\lambda_{k}+\sigma_{t}^{2}})c_{k}(t)\
.$ (64)

Using the constraint that $\alpha_{t}^{2}+\sigma_{t}^{2}=1$, this becomes

$\displaystyle\dot{c}_{k}(t)$
$\displaystyle=\frac{\beta_{t}\alpha_{t}^{2}(1-\lambda_{k})}{\alpha_{t}^{2}\lambda_{k}+\sigma_{t}^{2}}c_{k}(t)\
.$ (65)

For the orthogonal space component $\boldsymbol{\mathbf{y}}^{\perp}(t)$, it
will stay in the orthogonal space $\boldsymbol{\mathbf{U}}^{\perp}$, and more
specifically the 1D space spanned by the initial
$\boldsymbol{\mathbf{y}}^{\perp}(t)$—so, when going backward in time, its
dynamics is simply a downscaling of $\boldsymbol{\mathbf{y}}^{\perp}(T)$.

$\displaystyle\dot{\boldsymbol{\mathbf{y}}}^{\perp}(t)$
$\displaystyle=\beta_{t}(\frac{1}{\sigma_{t}^{2}}-1)\boldsymbol{\mathbf{y}}^{\perp}(t)=\beta_{t}\frac{1-\sigma_{t}^{2}}{\sigma_{t}^{2}}\boldsymbol{\mathbf{y}}^{\perp}(t)=\frac{\beta_{t}\alpha_{t}^{2}}{\sigma_{t}^{2}}\boldsymbol{\mathbf{y}}^{\perp}(t)\
.$ (66)

Combining these two results and solving the ODEs in the usual way, we have the
trajectory solution

$\displaystyle\boldsymbol{\mathbf{y}}_{t}$
$\displaystyle=d(t)\boldsymbol{\mathbf{y}}^{\perp}(T)+\sum_{k=1}^{r}c_{k}(t)\boldsymbol{\mathbf{u}}_{k}$
(67) $\displaystyle d(t)$
$\displaystyle=\exp\left(\int_{T}^{t}d\tau\frac{\beta_{\tau}\alpha_{\tau}^{2}}{\sigma_{\tau}^{2}}\right)$
(68) $\displaystyle c_{k}(t)$
$\displaystyle=c_{k}(T)\exp\left(\int_{T}^{t}d\tau\frac{\beta_{\tau}\alpha_{\tau}^{2}(1-\lambda_{k})}{\alpha_{\tau}^{2}\lambda_{k}+\sigma_{\tau}^{2}})\right)\
.$ (69)

The initial conditions are

$\displaystyle c_{k}(T)$
$\displaystyle=\boldsymbol{\mathbf{u}}_{k}^{T}\boldsymbol{\mathbf{y}}_{T}$
(70) $\displaystyle\boldsymbol{\mathbf{y}}^{\perp}(T)$
$\displaystyle=\boldsymbol{\mathbf{y}}_{T}-\sum_{k=1}^{r}c_{k}(T)\boldsymbol{\mathbf{u}}_{k}\
,\ \boldsymbol{\mathbf{y}}^{\perp}(T)\in\boldsymbol{\mathbf{U}}^{\perp}\ .$
(71)

To solve the ODEs, it is helpful to use a particular reparameterization of
time. In particular, consider a reparameterization in terms of $\alpha_{t}$
using the relationship $-\beta_{t}\alpha_{t}dt=d\alpha_{t}$. The integral we
must do is

$\displaystyle\int_{T}^{t}d\tau\frac{\beta_{\tau}\alpha_{\tau}^{2}(1-\lambda_{k})}{\alpha_{\tau}^{2}\lambda_{k}+\sigma_{\tau}^{2}}$
$\displaystyle=\int_{T}^{t}d\tau\frac{\beta_{\tau}\alpha_{\tau}^{2}(1-\lambda_{k})}{1+\alpha_{\tau}^{2}(\lambda_{k}-1)}$
(72)
$\displaystyle=\int_{\alpha_{T}}^{\alpha_{t}}d\alpha_{\tau}\frac{\alpha_{\tau}(\lambda_{k}-1)}{1+\alpha_{\tau}^{2}(\lambda_{k}-1)}$
(73)
$\displaystyle=\frac{1}{2}\log(1+\alpha_{\tau}^{2}(\lambda_{k}-1))\Bigr{|}^{\alpha_{t}}_{\alpha_{T}}$
(74)
$\displaystyle=\frac{1}{2}\log\left(\frac{1+(\lambda_{k}-1)\alpha_{t}^{2}}{1+(\lambda_{k}-1)\alpha_{T}^{2}}\right)\
.$ (75)

Note that taking $\lambda_{k}=0$ gives us the solution to dynamics in the
directions orthogonal to the manifold. We have

$\displaystyle c_{k}(t)$
$\displaystyle=c_{k}(T)\sqrt{\frac{1+(\lambda_{k}-1)\alpha_{t}^{2}}{1+(\lambda_{k}-1)\alpha_{T}^{2}}}$
(76) $\displaystyle d(t)$
$\displaystyle=\sqrt{\frac{1-\alpha_{t}^{2}}{1-\alpha_{T}^{2}}}\ .$ (77)

The time derivatives of these coefficients are

$\displaystyle\dot{c}_{k}(t)$
$\displaystyle=c_{k}(T)\frac{-(\lambda_{k}-1)\alpha_{t}^{2}\beta_{t}}{\sqrt{(1+(\lambda_{k}-1)\alpha_{T}^{2})({1+(\lambda_{k}-1)\alpha_{t}^{2})}}}$
(78) $\displaystyle\dot{d}(t)$
$\displaystyle=\frac{\alpha_{t}^{2}\beta_{t}}{\sqrt{(1-\alpha_{T}^{2})(1-\alpha_{t}^{2})}}\
.$ (79)

Finally, we can write out the explicit solution for the trajectory
$\boldsymbol{\mathbf{x}}_{t}$:

$\displaystyle\boldsymbol{\mathbf{x}}_{t}=\alpha_{t}\boldsymbol{\mathbf{\mu}}+d(t)\boldsymbol{\mathbf{y}}^{\perp}(T)+\sum_{k=1}^{r}c_{k}(t)\boldsymbol{\mathbf{u}}_{k}\
.$ (80)

We can see that there are three terms: 1)
$\alpha_{t}\boldsymbol{\mathbf{\mu}}$, an increasing term that scales up to
the mean $\boldsymbol{\mathbf{\mu}}$ of the distribution; 2)
$d(t)\boldsymbol{\mathbf{y}}^{\perp}(T)$, a decaying term downscaling the
residual part of the initial noise vector, which is orthogonal to the data
manifold; and 3) the $c_{k}(t)\boldsymbol{\mathbf{u}}_{k}$ sum, each term of
which has independent dynamics.

We also now have the analytical solution for the projected outcome:

$\displaystyle\hat{\boldsymbol{\mathbf{x}}}_{0}(\boldsymbol{\mathbf{x}}_{t})-\boldsymbol{\mathbf{\mu}}$
$\displaystyle=\frac{1}{\alpha_{t}}\boldsymbol{\mathbf{U}}\tilde{\boldsymbol{\mathbf{\Lambda}}}_{t}\boldsymbol{\mathbf{U}}^{T}(\boldsymbol{\mathbf{x}}_{t}-\alpha_{t}\boldsymbol{\mathbf{\mu}})$
(81)
$\displaystyle=\sum_{k=1}^{r}c_{k}(t)\frac{\alpha_{t}\lambda_{k}}{\alpha_{t}^{2}\lambda_{k}+\sigma_{t}^{2}}\boldsymbol{\mathbf{u}}_{k}$
$\displaystyle=\sum_{k=1}^{r}c_{k}(T)\frac{\alpha_{t}\lambda_{k}}{\sqrt{(\alpha_{t}^{2}\lambda_{k}+\sigma_{t}^{2})(\alpha_{T}^{2}\lambda_{k}+\sigma_{T}^{2})}}\boldsymbol{\mathbf{u}}_{k}\
.$

Similarly, we can write out the exact formula for the trajectory’s tangent
vector $\dot{\boldsymbol{\mathbf{x}}}_{t}$:

$\displaystyle\dot{\boldsymbol{\mathbf{x}}}_{t}=$
$\displaystyle\dot{\alpha_{t}}\boldsymbol{\mathbf{\mu}}+\dot{d(t)}\boldsymbol{\mathbf{y}}^{\perp}(T)+\sum_{k=1}^{r}\dot{c_{k}}(t)\boldsymbol{\mathbf{u}}_{k}$
(82) $\displaystyle=$
$\displaystyle-\alpha_{t}\beta_{t}\boldsymbol{\mathbf{\mu}}+\frac{\alpha_{t}^{2}\beta_{t}}{\sqrt{(1-\alpha_{T}^{2})(1-\alpha_{t}^{2})}}\boldsymbol{\mathbf{y}}^{\perp}(T)-\sum_{k=1}^{r}c_{k}(T)\frac{(\lambda_{k}-1)\alpha_{t}^{2}\beta_{t}}{\sqrt{(1+(\lambda_{k}-1)\alpha_{T}^{2})({1+(\lambda_{k}-1)\alpha_{t}^{2})}}}\boldsymbol{\mathbf{u}}_{k}\
.$

To simply notation, we define two functions

$\displaystyle\psi(t,\lambda)$
$\displaystyle:=\sqrt{\frac{\sigma_{t}^{2}+\alpha_{t}^{2}\lambda}{\sigma_{T}^{2}+\alpha_{T}^{2}\lambda}}$
(83) $\displaystyle\xi(t,\lambda)$
$\displaystyle:=\frac{\alpha_{t}\lambda}{\sqrt{(\alpha_{t}^{2}\lambda+\sigma_{t}^{2})(\alpha_{T}^{2}\lambda+\sigma_{T}^{2})}}\
.$ (84)

In terms of these functions, our solutions can be written as

$\displaystyle\mathbf{x}_{t}=$
$\displaystyle\alpha_{t}\mathbf{\mu}+\psi(t,0)\mathbf{y}^{\perp}(T)+\sum_{k=1}^{r}\psi(t,\lambda_{k})c_{k}(T)\mathbf{u}_{k}$
(85) $\displaystyle\hat{\mathbf{x}}_{0}(\mathbf{x}_{t})=$
$\displaystyle\mathbf{\mu}+\sum_{k=1}^{r}\xi(t,\lambda_{k})c_{k}(T)\mathbf{u}_{k}\
.$ (86)

The dynamics of $\psi(t,\lambda_{k})$ and $\xi(t,\lambda_{k})$ are plotted in
main text (Fig.4). Finally, we note that there are various special
relationships between $\psi(t,\lambda)$, $\phi(t,\lambda)$, and
$\xi(t,\lambda)$:

$\displaystyle\lambda\frac{d}{dt}\psi(t,\lambda)=-(\lambda-1)\beta_{t}\alpha_{t}\xi(t,\lambda)$
(87)
$\displaystyle\xi(t,\lambda)=\frac{1}{\alpha_{t}}\psi(t,\lambda)\phi(t,\lambda)\
.$ (88)

##  Appendix H Reverse diffusion dynamics look approximately like a 2D
rotation

In this section, we derive Eq. 4.2 from the main text, which quantifies the
extent to which the dynamics of $\boldsymbol{\mathbf{x}}_{t}$ look like a
rotation within a 2D plane spanned by $\boldsymbol{\mathbf{x}}_{0}$ (the
reverse diffusion endpoint) and $\boldsymbol{\mathbf{x}}_{T}$ (the initial
noise). We will derive the formula by assuming that the training set consists
of a single Gaussian mode, but will comment later on the fact that this
assumption can be relaxed somewhat.

Assume that reverse diffusion begins at time $T$ with $\alpha_{T}\approx 0$,
and ends at time $t=0$ with $\alpha_{0}=1$. Using Eq. 15, at some intermediate
time $t$ we have that

$\boldsymbol{\mathbf{x}}_{t}=\alpha_{t}\boldsymbol{\mathbf{\mu}}+\sqrt{\frac{1-\alpha_{t}^{2}}{1-\alpha_{T}^{2}}}\boldsymbol{\mathbf{y}}^{\perp}(T)+\sum_{k=1}^{r}\sqrt{\frac{1+(\lambda_{k}-1)\alpha_{t}^{2}}{1+(\lambda_{k}-1)\alpha_{T}^{2}}}c_{k}(T)\boldsymbol{\mathbf{u}}_{k}\
.$ (89)

It is also true, by substituting $t=0$ and $t=T$, that

$\begin{split}\boldsymbol{\mathbf{\mu}}&=\boldsymbol{\mathbf{x}}_{0}-\sum_{k=1}^{r}\sqrt{\frac{\lambda_{k}}{1+(\lambda_{k}-1)\alpha_{T}^{2}}}c_{k}(T)\boldsymbol{\mathbf{u}}_{k}\\\
\boldsymbol{\mathbf{y}}^{\perp}(T)&=\boldsymbol{\mathbf{x}}_{T}-\alpha_{T}\boldsymbol{\mathbf{\mu}}-\sum_{k=1}^{r}c_{k}(T)\boldsymbol{\mathbf{u}}_{k}\
.\end{split}$ (90)

Using these two equations, we can rewrite Eq. 89 as

$\begin{split}\boldsymbol{\mathbf{x}}_{t}&=\alpha_{t}\boldsymbol{\mathbf{\mu}}+\sqrt{\frac{1-\alpha_{t}^{2}}{1-\alpha_{T}^{2}}}\left[\boldsymbol{\mathbf{x}}_{T}-\alpha_{T}\boldsymbol{\mathbf{\mu}}-\sum_{k=1}^{r}c_{k}(T)\boldsymbol{\mathbf{u}}_{k}\right]+\sum_{k=1}^{r}\sqrt{\frac{1+(\lambda_{k}-1)\alpha_{t}^{2}}{1+(\lambda_{k}-1)\alpha_{T}^{2}}}c_{k}(T)\boldsymbol{\mathbf{u}}_{k}\\\
&=\left[\alpha_{t}-\alpha_{T}\sqrt{\frac{1-\alpha_{t}^{2}}{1-\alpha_{T}^{2}}}\right]\boldsymbol{\mathbf{\mu}}+\sqrt{\frac{1-\alpha_{t}^{2}}{1-\alpha_{T}^{2}}}\boldsymbol{\mathbf{x}}_{T}+\sum_{k=1}^{r}\left\\{\sqrt{\frac{1+(\lambda_{k}-1)\alpha_{t}^{2}}{1+(\lambda_{k}-1)\alpha_{T}^{2}}}-\sqrt{\frac{1-\alpha_{t}^{2}}{1-\alpha_{T}^{2}}}\right\\}c_{k}(T)\boldsymbol{\mathbf{u}}_{k}\\\
&=\left[\alpha_{t}-\alpha_{T}\sqrt{\frac{1-\alpha_{t}^{2}}{1-\alpha_{T}^{2}}}\right]\boldsymbol{\mathbf{x}}_{0}+\sqrt{\frac{1-\alpha_{t}^{2}}{1-\alpha_{T}^{2}}}\boldsymbol{\mathbf{x}}_{T}+\boldsymbol{\mathbf{R}}(t)\end{split}$
(91)

where the remainder term $\boldsymbol{\mathbf{R}}(t)$ is equal to

$\boldsymbol{\mathbf{R}}(t)=\sum_{k=1}^{r}\left\\{\sqrt{\frac{1+(\lambda_{k}-1)\alpha_{t}^{2}}{1+(\lambda_{k}-1)\alpha_{T}^{2}}}-\sqrt{\frac{1-\alpha_{t}^{2}}{1-\alpha_{T}^{2}}}-\left[\alpha_{t}-\alpha_{T}\sqrt{\frac{1-\alpha_{t}^{2}}{1-\alpha_{T}^{2}}}\right]\sqrt{\frac{\lambda_{k}}{1+(\lambda_{k}-1)\alpha_{T}^{2}}}\right\\}c_{k}(T)\boldsymbol{\mathbf{u}}_{k}\
.$ (92)

The expression simplifies somewhat if we take $\alpha_{T}=0$. Doing so, we
obtain the equation seen in the main text:

$\begin{split}\boldsymbol{\mathbf{x}}_{t}\approx\
&\alpha_{t}\boldsymbol{\mathbf{x}}_{0}+\sqrt{1-\alpha_{t}^{2}}\
\boldsymbol{\mathbf{x}}_{T}+\sum_{k=1}^{r}\left\\{\sqrt{\sigma_{t}^{2}+\lambda_{k}\alpha_{t}^{2}}-\alpha_{t}\sqrt{\lambda_{k}}-\sigma_{t}\right\\}c_{k}(T)\boldsymbol{\mathbf{u}}_{k}\
.\end{split}$

If $\alpha_{T}$ is not very close to zero, there may be a small first-order
correction to the rotation-like dynamics, even if the sum terms are
approximately zero.

Although we derived this formula by assuming a single Gaussian mode, its form
does not actually depend on the details of the nearest mode. In the landscape
picture, although the mean and variance of the nearest mode regularly change
in a discontinuous way, the form of our rotation equation should stay the
same.