# The Hidden Linear Structure in
Score-Based Models and its Application

Binxu Wang
Department of Neurobiology
Harvard Medical School
Boston, MA
<EMAIL_ADDRESS>
&John J. Vastola
Department of Neurobiology
Harvard Medical School
Boston, MA
<EMAIL_ADDRESS>
They are currently affiliated with The Kempner Institute for the Study of
Natural and Artificial Intelligence at Harvard University.

###### Abstract

Score-based models have achieved remarkable results in the generative modeling
of many domains. By learning the gradient of smoothed data distribution, they
can iteratively generate samples from complex distribution e.g. natural
images. However, is there any universal structure in the gradient field that
will eventually be learned by any neural network? Here, we aim to find such
structures through a normative analysis of the score function. First, we
derived the closed-form solution to the scored-based model with a Gaussian
score. We claimed that for well-trained diffusion models, the learned score at
a high noise scale is well approximated by the linear score of Gaussian. We
demonstrated this through empirical validation of pre-trained images diffusion
model and theoretical analysis of the score function. This finding enabled us
to precisely predict the initial diffusion trajectory using the analytical
solution and to accelerate image sampling by 15-30% by skipping the initial
phase without sacrificing image quality. Our finding of the linear structure
in the score-based model has implications for better model design and data
pre-processing.

##  1 Motivation: Universal structures of score

The diffusion model or score-based model learns the gradient of the smoothed
data distribution (i.e. score), which is used to guide data sampling.
Mathematically, the score function is a dynamic vector field in the domain of
data, which changes with respect to time or noise scale in the diffusion
process. Obviously, the score function depends on the points in the dataset,
that’s why a neural network is trained to approximate it. However, governed by
the diffusion process, is there any universal structure of the score field
that we can leverage?

In this paper, we claimed that for diffusion models, the early phase of the
score vector field is dominated by a linear structure, i.e. the score of
Gaussian approximation of the data. We derived the closed-form solution to the
Gaussian score model. We found that by leveraging this linear structure and
the solution of diffusion trajectory, we can predict the early phase of the
diffusion sampling trajectory and accelerate diffusion sampling by skipping
this phase entirely, without affecting image quality. This finding has
implications for how to train or design diffusion models.

##  2 Formulation: Score-based model

We will use the unifying framework and notations in Karras et al., (2022)
throughout this paper. The score function $\nabla\log p(\mathbf{x},\sigma)$ is
defined as the gradient of the log data distribution convolved by a Gaussian
$\mathcal{N}(0,\sigma^{2}I)$. Here, we consider the deterministic sampler,
i.e. sampling by integrating the probability flow Ordinary Differential
Equation (ODE), with $\sigma_{t}$ a decreasing function of time $t$.

$\displaystyle d\mathbf{x}=-\dot{\sigma}_{t}\sigma_{t}\nabla_{\mathbf{x}}\log
p(\mathbf{x},\sigma_{t})dt$ (1)

The score function is learned by minimizing the denoising score-matching
objective in Eq. 2. Following Karras et al., (2022), we reparametrize the
score with the optimal ‘denoiser’ $D_{\theta}$ using a neural network.

$\displaystyle\mathcal{L}_{\theta}=\mathbb{E}_{\mathbf{y}\sim
p_{data}}\mathbb{E}_{\mathbf{n}\sim\mathcal{N}(0,\sigma^{2}I)}\|D_{\theta}(\mathbf{y}+\mathbf{n},\sigma)-\mathbf{y}\|^{2}_{2}.\;\Rightarrow\nabla_{\mathbf{x}}\log
p(\mathbf{x},\sigma)\approx(D_{\theta}(\mathbf{x},\sigma)-\mathbf{x})/\sigma^{2}$
(2)

Here, we will present results without time-dependent data scaling, but with
little modification, all the results applied to models with data scaling e.g.
DDPM (Ho et al.,, 2020), DDIM (Song et al.,, 2020), VP-ODE (Song et al.,,
2021).

Figure 1: Gaussian score function approximates the learned neural score A.
Conceptual plan of the paper. B. Approximation error of learned neural score
by different analytical scores: Gaussian, isotropic, Gaussian mixture, and
‘exact’ point cloud. Note the log scale of error on y axis. C. Deviation of
diffusion trajectory from the solution of analytical scores.

##  3 Theory: Closed-form solution to diffusion with Gaussian score

First, we showed that the diffusion trajectory is analytically solvable for
Gaussian training distribution. Let $\mathbf{x}_{t}\in\mathbb{R}^{D}$, and the
mean and covariance of the Gaussian be $\mathbf{\mu}$ and $\mathbf{\Sigma}$.
Assuming $\mathbf{\Sigma}$ has rank $r\leq D$, it has a compact singular value
decomposition (SVD)
$\mathbf{\Sigma}=\mathbf{U}\mathbf{\Lambda}\mathbf{U}^{T}$, where
$\mathbf{U}=[\mathbf{u}_{1},...,\mathbf{u}_{r}]$ is a $D\times r$ semi-
orthogonal matrix. The columns of $\mathbf{U}$ are the principal axes along
which the mode varies, and their span comprises the ‘data manifold’. At noise
scale $\sigma$, the score is that of
$\mathcal{N}(\mu,\Sigma+\sigma^{2}I)=(\sigma^{2}I+\Sigma)^{-1}(\mathbf{\mu}-\mathbf{x})$.
With Woodbury matrix inversion trick (Woodbury,, 1950), it simplifies as,

$\displaystyle\nabla\log
p(\mathbf{x};\sigma)=\frac{1}{\sigma^{2}}(I-\mathbf{U}\tilde{\Lambda}_{\sigma}\mathbf{U}^{T})(\mathbf{\mu}-\mathbf{x})\;,\;\;\tilde{\Lambda}_{\sigma}$
$\displaystyle=\mbox{diag}\big{[}\frac{\lambda_{k}}{\lambda_{k}+\sigma^{2}}\big{]}$
(3)

which is a linear function of state. It consists of an isotropic term pointing
to the mean $\mathbf{\mu}$, and a 2nd term depending on the PC spectrum
structure $\mathbf{U},\mathbf{\Lambda}$. With this linear score function,
probability flow (Eq. 1) becomes a linear ODE, which has the following closed-
form solution (Appendix 9.1).

$\displaystyle\mathbf{x}_{t}=$
$\displaystyle\mathbf{\mu}+\frac{\sigma_{t}}{\sigma_{T}}\
\mathbf{x}^{\perp}_{T}+\sum_{k=1}^{r}\psi(t,\lambda_{k})c_{k}(T)\mathbf{u}_{k}\;\;$
$\displaystyle\psi(t,\lambda_{k}):=\sqrt{\frac{\sigma_{t}^{2}+\lambda_{k}}{\sigma_{T}^{2}+\lambda_{k}}}$
(4) $\displaystyle\mathbf{x}^{\perp}_{T}:=$ $\displaystyle\
(\mathbf{I}-\mathbf{U}\mathbf{U}^{T})(\mathbf{x}_{T}-\mathbf{\mu})\;$
$\displaystyle c_{k}(T):=\ \mathbf{u}_{k}^{T}(\mathbf{x}_{T}-\mathbf{\mu})\ .$
(5)

The solution has 3 components, corresponding to the 3 parts of the initial
condition $\mathbf{x}_{T}$: 1) the distribution mean $\mathbf{\mu}$ stays; 2)
the off-manifold component $\mathbf{x}^{\perp}$ shrinks proportional to the
noise scale $\psi(t,0)=\frac{\sigma_{t}}{\sigma_{T}}$; and 3) the projections
along each eigenvector $\mathbf{u}_{k}$ moves governed by
$\psi(t,\lambda_{k})$. Similarly, the ideal denoising image along the
trajectory can be expressed as

$\displaystyle D(\mathbf{x}_{t},\sigma_{t})=\
\mathbf{\mu}+\sum_{k=1}^{r}\xi(t,\lambda_{k})c_{k}(T)\mathbf{u}_{k}\hskip
36.135pt\xi(t,\lambda):=\frac{\lambda}{\sqrt{(\lambda+\sigma_{t}^{2})(\lambda+\sigma_{T}^{2})}}\
.$ (6)

Note that even this simple closed-form solution can explain many phenomena of
diffusion models, including 1) the sampling trajectory of variance-preserving
diffusion models tend to be low dimensional, and more specifically, resembling
a two-dimensional rotation on the plane of initial and final state; 2) the
outlines or low-frequency aspect of images are determined earlier in the
generation process; 3) earlier and on manifold perturbations tend to have a
larger impact on the image generation process. Though not directly related to
the main claim of this paper, we encourage the readers to further explore
these connections in Wang and Vastola, (2023). This showed that qualitatively,
the sampling dynamics of diffusion models are largely dominated by the
Gaussian solution.

Figure 2: Validation and application of the Gaussian score approximation. A.
The denoiser images $D(\mathbf{x}_{t},t)$ along a diffusion sampling
trajectory compared with the Gaussian solution with the same initial condition
$\mathbf{x}_{T}$. B. Samples generated by the EDM model, Gaussian solution,
and the ‘exact’ scores from the same initial condition. C. Sampled image as a
function of skipped steps with the hybrid method combining Gaussian theory
with Heun’s method. D. Image quality (FID score) of the hybrid method as a
function of NFE, skipped steps, and skipped noise scale (see Appendix 7.5).

##  4 Validation: The Gaussian score structure in diffusion models

Most distributions useful enough to be fit by generative models are not
Gaussian. So how is the analytical solution of the Gaussian score model
related to practical diffusion models? Here we empirically tested our claim:
for pre-trained image diffusion models, the score in the early high-noise
phase is dominated by the linear score of the Gaussian approximation of data.
In Appendix 9.6, we also argued theoretically that, at a high enough noise
scale, the score function of bounded point cloud is indistinguishable from
that of the Gaussian with matching mean and covariance.

We picked several image datasets (CIFAR10, FFHQ64, AFHQv2-64 (Choi et al.,,
2020)) and the diffusion models pre-trained on them using the optimal
configuration in Karras et al., (2022). For each dataset, we computed several
analytical approximations of the score and compared that with the score
function of the neural network (Eq. 2): 1) Isotropic score (Iso),
$(\mathbf{\mu}-\mathbf{x})/\sigma^{2}$, which is the first order term of the
Gaussian score; 2) Gaussian score (Eq. 3), computed with the mean and
covariance of the training dataset; 3) Point cloud score (Exact), computed for
Gaussian mixture model with a Gaussian mode with negligible width at each data
point; 4) Per-class Gaussian Mixture Score (GMM), when the dataset has class
labels (e.g. CIFAR), we computed the score for Gaussian Mixture Model where
each mode has mean and covariance for each class (details in Sec.7.2). We
compared these score vector fields with the neural scores, and the diffusion
trajectories guided by them.

#### Gaussian score predicts the learned score at high noise

First, we validated the score vector prediction. For each noise scale
$\sigma$, we evaluated the neural and analytical scores at random points
sampled from $\mathcal{N}(0,\sigma^{2}I)$. We computed the fraction of
unexplained variance (residual) as a function of the noise scale for each
analytical score approximator (Fig. 3). For all three spaces, at most noise
levels ($\sigma>1.08$), the Gaussian score explains almost all variance (>
99%) of the neural score. As a control, at all noise levels, the Gaussian
score explains more variance than the isotropic score without information of
the covariance. Further, for CIFAR10, at all levels, the 10-class Gaussian
mixture score predicts neural score slightly better than the Gaussian score,
which shows that adding more modes is indeed helpful to capture the details of
the score field.

#### Neural score deviates from ‘exact’ score at low noise

However, this trend breaks down when taken to an extreme. In all three
datasets, the exact point cloud score shows a slight improvement over the
neural score in high-noise conditions (Fig. 3); however, it diverges
significantly from the neural score under low-noise conditions. This dichotomy
implies that in high-noise regimes, the neural score captures more nuanced
features of the data points than a single Gaussian does. Conversely, in low-
noise regimes, the neural score diverges from the ‘exact’ score, resembling
more closely a Gaussian or Gaussian mixture score. This indicates that the
neural network learned a smoother score field, particularly in low-noise
conditions (Kadkhodaie et al.,, 2023; Scarvelis et al.,, 2023), thereby
mitigating overfitting and preventing memorization, since learning the exact
score would lead to generating samples exclusively from the training set.

#### Gaussian solution predicts early diffusion trajectory

Next, we tested how well the solutions to ODE (Eq. 1) with the closed-form
scores predicted the sampling trajectories of actual diffusion models. We
found that the early phase of reverse diffusion is well-predicted by the
Gaussian solution (Fig. 4, 5). The predicted and actual trajectory stayed
close (MSE < 0.01) for around 9 sampling steps (noise scale $\sigma_{t}=1.92$)
for CIFAR10 and 17 steps ($\sigma_{t}=4.37$) for FFHQ and AFHQ until
diverging. This roughly matches the scale where the Gaussian approximation of
the score field breaks down, and the score needs to be approximated by that of
a more complicated distribution (Fig.3). Visually, as the low-frequency
information is determined early on Ho et al., (2020), the ‘layout’ or shading
of the final image is predicted by Gaussian solution, while high-frequency
details such as edges are less well-predicted (Fig. 2 A,B).

#### Neural diffusion trajectory deviates from the ‘exact’ score at low noise

In line with our score field observations, the early diffusion trajectory
aligns better with the exact score trajectory than with the Gaussian score
trajectory. However, the late trajectory based on neural score diverges from
the exact score trajectory, yielding notably different samples in the end
(Fig. 2 B).

In summary, these results show that at the high noise regime, the neural score
function is dominated by the Gaussian approximation, but with details better
predicted by the more elaborate Gaussian mixture score or ‘exact’ score. But
at low noise regime, the actual neural score function learned is different
from the ‘exact’ score, and it’s more similar to the ‘smoother’ score of
Gaussian or Gaussian mixture, possibly due to the regularizing effect of the
neural network function approximator.

##  5 Application: Accelerating Diffusion with analytical teleportation.

The Gaussian analytical solution provides a surprisingly good approximation to
the early sampling trajectory, thus a direct application is to accelerate
diffusion by ‘teleportation’. Namely, instead of evaluating the neural score
function and integrating the probability flow ODE, we can directly evaluate
$\mathbf{x}_{t}$ using the Gaussian solution with $\mathbf{\mu}$ and
$\mathbf{\Sigma}$ of the training set. In principle, this speedup can be
combined with any deterministic sampler. We showcase its effectiveness with
the optimized 2nd order Heun sampler from Karras et al., (2022), which has
demonstrated near state-of-the-art image quality and efficiency. For each
initial state, we skipped to different times with Eq.4 and continued with the
Heun sampler from that time (see Sec.7.4, Alg.1). We tested this hybrid
sampler on unconditional diffusion models of CIFAR-10, FFHQ-64, and AFHQv2-64
Karras et al., (2022). We found that, consistently, we can save 15-30% of
neural function evaluation (NFE) with less than 3% increase in the FID score
(Fig.2 C,D). Intriguingly, for CIFAR10 and AFHQ, skipping steps can even
reduce the FID score by 1% (full results in Sec.8.2 Fig.6, Tab.2,3,4). We also
noticed that the number of skippable steps roughly corresponds to the noise
scale where the neural score deviates from the Gaussian proxy.

##  6 Discussion

In summary, we showed that even for natural images, the learned score at a
high noise regime is well predicted by the linear score of Gaussian; at the
low noise regime, it is better approximated by a smoother Gaussian (Mixture)
score instead of an ‘exact’ score. We leveraged this finding to accelerate
diffusion sampling. This result bears further implications for training and
designs of diffusion models.

Implications on Model Design Since the score function to be learned is
dominated by the Gaussian score, we can we can directly build this linear
structure into the neural network to assist learning, e.g. to add a by-pass
pathway based on the Gaussian approximation score of training distribution,
and let the neural network learn the nonlinear residual to the Gaussian terms
only.

Implications on Training Distribution Further, we may be able to reshape the
training distribution to assist learning. We hypothesize that if we pre-
condition the target distribution by whitening its spectrum, then the target
score become more isotropic, and the score network may converge faster. This
might explain the training efficiency of latent diffusion models Rombach et
al., (2022): with the KL regularization, the auto-encoder not only compressed
the state space but also ‘whiten’ the image distribution by morphing it closer
to isotropic Gaussian distribution. Similar whitening tricks is worth
exploring for general diffusion models.

## Acknowledgments and Disclosure of Funding

B.W. was funded by Victoria Quan fellowship at Harvard Medical School, and
Kempner research fellowship. We thank Carlos R. Ponce lab, and O2 Cluster at
Harvard Medical School for providing computational resources for this work. We
thank Zhengdao Chen, Jacob Zavatone-Veth, Cengiz Pehlevan, Yao Fu, Hao Sun for
providing feedback to an early version of this work.

## References

  * Choi et al., (2020) Choi, Y., Uh, Y., Yoo, J., and Ha, J.-W. (2020).  Stargan v2: Diverse image synthesis for multiple domains.  In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pages 8188–8197.
  * Griffiths, (2005) Griffiths, D. J. (2005).  Introduction to electrodynamics.
  * Ho et al., (2020) Ho, J., Jain, A., and Abbeel, P. (2020).  Denoising diffusion probabilistic models.  Advances in Neural Information Processing Systems, 33:6840–6851.
  * Kadkhodaie et al., (2023) Kadkhodaie, Z., Guth, F., Simoncelli, E. P., and Mallat, S. (2023).  Generalization in diffusion models arises from geometry-adaptive harmonic representation.  arXiv preprint arXiv:2310.02557.
  * Karras et al., (2022) Karras, T., Aittala, M., Aila, T., and Laine, S. (2022).  Elucidating the design space of diffusion-based generative models.  arXiv preprint arXiv:2206.00364.
  * Liu et al., (2022) Liu, L., Ren, Y., Lin, Z., and Zhao, Z. (2022).  Pseudo numerical methods for diffusion models on manifolds.  arXiv preprint arXiv:2202.09778.
  * Rombach et al., (2022) Rombach, R., Blattmann, A., Lorenz, D., Esser, P., and Ommer, B. (2022).  High-resolution image synthesis with latent diffusion models.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 10684–10695.
  * Scarvelis et al., (2023) Scarvelis, C., Borde, H. S. d. O., and Solomon, J. (2023).  Closed-form diffusion models.  arXiv preprint arXiv:2310.12395.
  * Song et al., (2020) Song, J., Meng, C., and Ermon, S. (2020).  Denoising diffusion implicit models.  arXiv preprint arXiv:2010.02502.
  * Song et al., (2021) Song, Y., Sohl-Dickstein, J., Kingma, D. P., Kumar, A., Ermon, S., and Poole, B. (2021).  Score-based generative modeling through stochastic differential equations.  In International Conference on Learning Representations.
  * Wang and Vastola, (2023) Wang, B. and Vastola, J. J. (2023).  Diffusion models generate images like painters: an analytical theory of outline first, details later.  arXiv preprint arXiv:2303.02490.
  * Woodbury, (1950) Woodbury, M. A. (1950).  Inverting modified matrices.  Department of Statistics, Princeton University.

##  7 Detailed Methods

###  7.1 Image datasets and Models

Table 1: Specifications of the image datasets and diffusion models Dataset name | Num Samples | Resolution | Pre-trained Model Spec
---|---|---|---
CIFAR10 | 50000 | 32 | edm-cifar10-32x32-uncond-vp
FFHQ64 | 70000 | 64 | edm-ffhq-64x64-uncond-vp
AFHQv2-64 | 15803 | 64 | edm-afhqv2-64x64-uncond-vp

###  7.2 Idealized Score Approximations

In the paper, we compared several analytical approximations of the score. We
listed their formula below.

Isotropic score, only depends on the mean of data, isotropically pointing
towards $\mu$.

$\frac{\mu-\mathbf{x}}{\sigma^{2}}$ (7)

This is equivalent to approximating the whole training dataset with its mean.

Gaussian score, for a Gaussian distribution of
$\mathcal{N}(\mathbf{x};\mathbf{\mu},\mathbf{\Sigma})$, its score is

$\displaystyle\nabla_{\mathbf{x}}\log\mathcal{N}(\mathbf{x};\mathbf{\mu},\sigma^{2}\mathbf{I}+\mathbf{\Sigma})$
(8) $\displaystyle=(\sigma^{2}I+\Sigma)^{-1}(\mu-\mathbf{x})$ (9)
$\displaystyle=\frac{1}{\sigma^{2}}(I-U\tilde{\Lambda}_{\sigma}U^{T})(\mathbf{\mu}-\mathbf{x})$
(10)
$\displaystyle=\frac{\mu-\mathbf{x}}{\sigma^{2}}-\frac{1}{\sigma^{2}}U\tilde{\Lambda}_{\sigma}U^{T}(\mathbf{\mu}-\mathbf{x})$
(11) $\displaystyle\tilde{\Lambda}_{\sigma}$
$\displaystyle=diag\big{[}\frac{\lambda_{k}}{\lambda_{k}+\sigma^{2}}\big{]}$
(12)

it depends on the mean and covariance of the data. It added a non-isotropic
correction term based on covariance to the isotropic score.

Gaussian mixture score. For a general Gaussian mixture with $k$ component
$q(\mathbf{x})=\sum^{k}_{i}\pi_{i}\mathcal{N}(\mathbf{x};\mathbf{\mu}_{i},\mathbf{\Sigma}_{i})$,
its score function at noise scale $\sigma$ is the following,

$\begin{split}&\nabla_{\mathbf{x}}\log\Big{(}\sum^{k}_{i}\pi_{i}\mathcal{N}(\mathbf{x};\mathbf{\mu}_{i},\sigma^{2}\mathbf{I}+\mathbf{\Sigma}_{i})\Big{)}\\\
=&\sum_{i}-(\sigma^{2}\mathbf{I}+\mathbf{\Sigma}_{i})^{-1}(\mathbf{x}-\mathbf{\mu}_{i})\frac{\pi_{i}\mathcal{N}(\mathbf{x};\mathbf{\mu}_{i},\sigma^{2}\mathbf{I}+\mathbf{\Sigma}_{i})}{\sum^{k}_{j}\pi_{j}\mathcal{N}(\mathbf{x};\mathbf{\mu}_{j},\sigma^{2}\mathbf{I}+\mathbf{\Sigma}_{j}))}\\\
=&\sum_{i}-(\sigma^{2}\mathbf{I}+\mathbf{\Sigma}_{i})^{-1}(\mathbf{x}-\mathbf{\mu}_{i})w_{i}(\mathbf{x},\sigma)\
.\end{split}$ (13)

Where the covariance matrix for each Gaussian component can be eigen-
decomposed and inverted efficiently. It’s a weighted average of the score of
each Gaussian mode, with a softmax-like weighting function
$w_{i}(\mathbf{x},\sigma)$.

Exact point cloud score. For a set of data points $\\{\mathbf{y}_{i}\\}$, the
score of $\mathbf{x}$ at noise scale $\sigma$ is

$\displaystyle\nabla_{\mathbf{x}}\log\Big{(}\sum^{N}_{i}\frac{1}{N}\mathcal{N}(\mathbf{x};\mathbf{y}_{i},\sigma^{2}\mathbf{I}\Big{)}$
(14) $\displaystyle=$
$\displaystyle\frac{1}{\sigma^{2}}\left[-\mathbf{x}+\sum_{i}w_{i}(\mathbf{x},\sigma)\mathbf{y}_{i}\right]$
(15) $\displaystyle=$
$\displaystyle\frac{1}{\sigma^{2}}\left[-\mathbf{x}+\sum_{i}\frac{\exp\big{(}-\frac{1}{2\sigma^{2}}\|\mathbf{y}_{i}-\mathbf{x}\|^{2}\big{)}}{\sum_{j}\exp\big{(}-\frac{1}{2\sigma^{2}}\|\mathbf{y}_{j}-\mathbf{x}\|^{2}\big{)}}\mathbf{y}_{i}\right]\
.$ (16) $\displaystyle w_{i}(\mathbf{x},\sigma):=$
$\displaystyle\frac{\exp\big{(}-\frac{1}{2\sigma^{2}}\|\mathbf{y}_{i}-\mathbf{x}\|^{2}\big{)}}{\sum_{j}\exp\big{(}-\frac{1}{2\sigma^{2}}\|\mathbf{y}_{j}-\mathbf{x}\|^{2}\big{)}}=\mbox{softmax}\Big{(}-\frac{1}{2\sigma^{2}}\|\mathbf{y}_{i}-\mathbf{x}\|^{2}\Big{)}$
(17)

###  7.3 Sampling Diffusion trajectory with Gaussian Mixture scores

For the Gaussian score, we are able to compute the whole sampling trajectory
analytically. But for a Gaussian mixture with more than one mode, we need to
use numerical integration. We evaluated the numerical score with the Gaussian
mixture model defined as above and integrated Eq. 1 with off-the-shelf Runge-
Kutta 4 integrator (solve_ivp from scipy). We chose $\sigma(t)=t$, and
integrate $t$ backward from sigma $80.0$ to $0.0$. To compare the trajectory
with the one sampled with the original Heun method, we evaluated the
trajectory at the same discrete time steps as the Karras et al., (2022) paper.
The $i$th noise level is the following,

$\sigma_{i}=\Big{(}\sigma_{max}^{1/\rho}+\frac{i}{n_{step}-1}(\sigma_{min}^{1/\rho}-\sigma_{max}^{1/\rho})\Big{)}^{\rho}$
(18)

we chose the same hyper parameter $\sigma_{min},\sigma_{max},\rho$ and
$n_{step}$ as the original EDM paper.

###  7.4 Hybrid sampling method

As we stated in the paper, the Hybrid sampling scheme can be combined with any
deterministic sampler (e.g. DDIM Song et al., (2020), PNDM Liu et al.,
(2022)). Here we detailed its implementation with Heun’s method.

0: Data: mean and covariance of dataset $\mu$, $\Sigma$ and its
eigendecomposition $\mathbf{\Sigma}=\mathbf{U}\mathbf{\Lambda}\mathbf{U}^{T}$

0: Original parameter: $\sigma_{min},\sigma_{max},\rho$

0: Parameter: Skip time/noise scale $t^{\prime}=\sigma_{skip}$

0: Input: $\mathbf{x}_{T}$

1: $\mathbf{x}_{t^{\prime}}\leftarrow$ Gaussian solution at $t^{\prime}$ using
$\mu$ and $\Sigma$:

$\mathbf{x}_{t^{\prime}}=\mathbf{\mu}+\frac{\sigma_{t^{\prime}}}{\sigma_{max}}(\mathbf{I}-\mathbf{U}\mathbf{U}^{T})(\mathbf{x}_{T}-\mathbf{\mu})+\sum_{k=1}^{r}\sqrt{\frac{\sigma_{t^{\prime}}^{2}+\lambda_{k}}{\sigma_{max}^{2}+\lambda_{k}}}\mathbf{u}_{k}\mathbf{u}_{k}^{T}(\mathbf{x}_{T}-\mathbf{\mu})$

2: Continue sampling using Heun method with $\sigma_{\text{max}}=t^{\prime}$
and initial state $\mathbf{x}_{t^{\prime}}$, get the final sample
$\mathbf{x}_{0}$.

3: return $\mathbf{x}_{0}$

Algorithm 1 Hybrid Sampling Method with Heun’s method

For a closer comparison with the baseline Heun method, we used the following
strategy for choosing $\sigma_{skip}$. The Heun method will sample a sequence
of $n_{step}$ noise levels
$[\sigma_{0},\sigma_{1},\sigma_{2},...\sigma_{n_{step}}]$, where
$\sigma_{0}=\sigma_{max}$ and $\sigma_{n_{step}}=\sigma_{min}$. For each
initial condition $\mathbf{x}_{T}$, we use the analytical solution to evaluate
$\mathbf{x}_{t^{\prime}}$ or integrate the probability flow ODE to time
$t^{\prime}$, where we chose $t^{\prime}$ as the $i$-th noise level
$\sigma_{i}$ (Eq. 18). Then we will skip the first $i$ step in the Heun method
and start at initial state $\mathbf{x}_{t^{\prime}}$. Note that Heun’s method
is a 2nd order method, which evaluates the neural score function twice for
each step in $\mathbf{x}_{t}$. Because of this, skipping $i$ time steps will
save $2i$ neural function evaluations (note x-axes in Fig. 6).

###  7.5 FID score computation and baseline

We used the same code for FID score computation as in Karras et al., (2022).
For each sampler, we sampled the same initial noise state $x_{T}$ with random
seeds 0-49999. We computed the FID score based on 50,000 samples.

For our baseline, we picked the same model configurations as reported in Tab.
2, specifically, Variance preserving (VP) config F, the unconditional model
for CIFAR10, FFHQ 64 and AFHQv2 64. The default sampling steps are 18 time
steps ($NFE=35$) for CIFAR10, and 40 time steps ($NFE=79$) for FFHQ and AFHQ.

##  8 Extended Results

###  8.1 Detailed Validation Results

Here we presented the full results for all datasets: approximating neural
scores with analytical scores (Fig.3), deviations between sampling trajectory
(Fig.4) and denoiser trajectory (Fig.5) guided by neural score and analytical
score.

Figure 3: Approximation error of score learned by score neural network with
various analytical approximations. Note the log scale on y. Figure 4:
Deviation between the state trajectory $\mathbf{x}_{t}$ of EDM neural network
and the analytical score. The thick line denotes the mean over initial
conditions; the shaded area denotes 25%, 75% quantile range over the initial
conditions. Figure 5: Deviation between the denoiser $D(\mathbf{x}_{t},t)$ of
EDM neural network and the analytical score. The thick line denotes the mean
over initial conditions; the shaded area denotes 25%, 75% quantile range over
the initial conditions. Figure 6: Image quality as a function of skipping
steps for hybrid sampling approach. Note the main x-axes are the number of
Neural Function Evaluation (NFE); the secondary x-axes are the number of
skipping steps from the Heun sampler; the tertiary-axes are the time or noise
level $\sigma_{skip}$ at which we evaluate the Gaussian solution. See
Tab.4,2,3 for numbers.

###  8.2 Detailed Gaussian teleportation FID scores results

The effects of applying the Gaussian analytical teleportation on the FID
scores are presented below (Fig.6, Tab.3,2,4).

For CIFAR-10 model, we found teleportation not only reduces the number of
required neural function calls, but also lowers the FID score. Admittedly, the
FID score effect is quite small on the CIFAR-10 model: it goes from 1.958 (no
teleportation) to 1.933 (5 skipping steps), which is around 1% decrease.
Visually, the difference in sample images is mostly negligible; there are some
changes, but they are plausible changes in image details. Hence, teleportation
saves 29% of NFE and improves the sample quality, without changing the model
or the sampler. Further, skipping 6 steps will save 34% NFE but increase the
FID score by less than .2%.

For AFHQv2-64, the teleportation can save 25-30% of NFE (10-12 steps skipped)
or reduce the FID by 2% (6 steps skipped). For FFHQ64, the teleportation saves
10-15% of NFE (4-6 steps skipped) with around 3% increase in FID.

The full table of FID scores as a function of the number of skipped steps is
shown below.

The fact that replacing neural network evaluations with the analytical
solution can even improve the already low FID score (albeit very slightly in
the tests we just ran) is intriguing, and it is not immediately obvious why
this is true. One possibility is that, though using the Inception feature
space, the FID score is still measuring the similarity of image distribution
with a Gaussian approximation. Thus using an explicit Gaussian approximation
of data to guide diffusion may help this metric.

Table 2: FFHQ64 FID with analytical teleportation Nskip | NFE | time/noise scale | FID
---|---|---|---
0 | 79 | 80.0 | 2.464
2 | 75 | 60.1 | 2.489
4 | 71 | 44.6 | 2.523
6 | 67 | 32.7 | 2.561
8 | 63 | 23.6 | 2.617
10 | 59 | 16.8 | 2.681
12 | 55 | 11.7 | 2.841
14 | 51 | 8.0 | 3.243
16 | 47 | 5.4 | 4.402
20 | 39 | 2.2 | 15.451
Table 3: AFHQV264 FID with analytical teleportation Nskip | NFE | time/noise scale | FID
---|---|---|---
0 | 79 | 80.0 | 2.043
2 | 75 | 60.1 | 2.029
4 | 71 | 44.6 | 2.016
6 | 67 | 32.7 | 2.003
8 | 63 | 23.6 | 2.005
10 | 59 | 16.8 | 2.026
12 | 55 | 11.7 | 2.102
14 | 51 | 8.0 | 2.359
16 | 47 | 5.4 | 3.206
20 | 39 | 2.2 | 14.442
Table 4: CIFAR10 FID with analytical teleportation Nskip | NFE | time/noise scale | FID
---|---|---|---
0 | 35 | 80.0 | 1.958
1 | 33 | 57.6 | 1.955
2 | 31 | 40.8 | 1.949
3 | 29 | 28.4 | 1.940
4 | 27 | 19.4 | 1.932
5 | 25 | 12.9 | 1.934
6 | 23 | 8.4 | 1.963
7 | 21 | 5.3 | 2.123
8 | 19 | 3.3 | 3.213
10 | 15 | 1.1 | 23.947
12 | 11 | 0.3 | 109.178

###  8.3 Additional Validation experiments with Variance Preserving Diffusion

Here we compared the Gaussian solutions for variance-preserving ODEs
(Sec.9.2), with the actual diffusion trajectories sampled by DDIM sampler. We
used pre-trained DDPM / DDIM models on MNIST, CIFAR10, LSUN church, CelebA
256. The dataset details are listed below (Tab. 5).

Table 5: Diffusion models used for 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
---|---|---
MNIST | dimpo/ddpm-mnist |  3x32x32=3072 ${\dagger}$
CIFAR-10 | google/ddpm-cifar10-32 | 3x32x32=3072
Lsun-Church | google/ddpm-church-256 | 3x256x256 = 196,608
CelebA-HQ | google/ddpm-celebahq-256 | 3x256x256 = 196,608

#### Teleportation results

For the MNIST and CIFAR-10 models, we can easily skip 40% of the initial steps
with the Gaussian solution without much of a perceptible change in the final
sample (Fig.7E, Fig.8D). Quantitatively, for CIFAR10 model, we found skipping
up to 40% of the initial steps can even slightly decrease the Frechet
Inception Distance score (FID), and hence improve the quality of generated
samples (Fig.7F). For models of higher resolution data sets like CelebA-HQ, we
need to be more careful; skipping more than 20% of the initial steps will
induce some perceptible distortions in the generated images (Fig.7E bottom),
which suggests that the Gaussian approximation is less effective for larger
images. The reason may have to do with a low-quality covariance matrix
estimate, which could arise from the number of training images being small
compared to the effective dimensionality of the image manifold.

Figure 7: Comparing analytical solutions of sampling trajectory with DDIM
diffusion model for CIFAR-10. A. $\hat{\mathbf{x}}_{0}(\mathbf{x}_{t})$ of a
DDIM trajectory and the Gaussian solution with the same initial condition
$\mathbf{x}_{T}$. B. Samples generated by DDIM and the analytical theories
from the same initial condition. C. Mean squared error between the
$\mathbf{x}_{t}$ trajectory of DDIM and Gaussian solution. D. Comparing the
state trajectory and final sample of three normative models (Gaussian, GMM,
exact) with DDIM. E. Hybrid sampling method combines Gaussian theory
prediction with DDIM. F. Image quality of the hybrid method (FID score) as
function of different numbers of skipped steps (see Appendix 7.5).

Figure 8:  Comparing diffusion dynamics guided by Gaussian and Gaussian
mixture score with the ones learned by a neural network (MNIST). A.
Unconditional model. In each row, it shows the image generated from single-
mode Gaussian, 10-mode Gaussian mixture, Delta GMM defined on the whole MNIST
dataset, and the unconditional diffusion model.
B. Conditional model of the 10 classes. In each row, from left to right, it
shows the image generated via single Gaussian, delta GMM defined on all
training data of that class, and the conditional diffusion model with
different classifier free guidance strength 0.0 - 5.0.
C. Unconditional model, Deviation between trajectory predicted by unimodal
Gaussian theory and 10-mode GMM and exact delta GMM. It shows that Gaussian
mixture the same effect on the diffusion trajectory in the first phase as a
matched unimodal Gaussian.
D. Unconditional model, Deviation between the trajectory predicted by
different GMM theories and that sampled by DDIM. It shows that the trajectory
sampled by DDIM is actually closer to that predicted by unimodal or 10 mode
GMM, than the exact delta gmm model. This means the model didn’t learn the
exact score, but some coarse-grained approximation to it.
E. Conditional model, Deviation between unimodal Gaussian theory and exact
delta GMM for different digits class. It shows some digits are better
approximated by Gaussian than others.
F. Conditional model, Deviation between the trajectory predicted by exact
delta GMM and that sampled by DDIM with different guidance scales. It shows
that a larger guidance scale push the sampling result from the exact training
distribution further.

##  9 Detailed Derivations

###  9.1 Detailed derivation to the solution of the Gaussian Diffusion

In this section, we derive the analytic solution $\mathbf{x}_{t}$ for the
reverse diffusion trajectory of the Gaussian score model. We will assume that
the score function corresponds to a Gaussian distribution whose mean is
$\mathbf{\mu}$ and covariance matrix is $\mathbf{\Sigma}$. We will also
assume, given that images are usually regarded as residing on a low-
dimensional manifold within pixel space, that the rank $r$ of the covariance
matrix may be less than the dimensionality $D$ of state space. Let
$\mathbf{\Sigma}=\mathbf{U}\mathbf{\Lambda}\mathbf{U}^{T}$ be the
eigendecomposition or compact SVD of the covariance matrix, where $\mathbf{U}$
is a $D\times r$ semi-orthogonal matrix whose columns are normalized (i.e.
$\mathbf{U}^{T}\mathbf{U}=\mathbf{I}_{r}$), and $\mathbf{\Lambda}$ is the
$r\times r$ diagonal eigenvalue matrix. Denote the $k$th column of
$\mathbf{U}$ by $\mathbf{u}_{k}$ and the $k$th diagonal element of
$\mathbf{\Lambda}$ by $\lambda_{k}$.

The probability flow ODE used for the EDM model reads,

$\displaystyle d\mathbf{x}$
$\displaystyle=-\sigma(\sigma^{2}I+\Sigma)^{-1}(\mu-\mathbf{x})d\sigma$ (19)
$\displaystyle d\mathbf{x}$
$\displaystyle=-\frac{1}{\sigma}(I-\mathbf{U}\tilde{\Lambda}_{\sigma}\mathbf{U}^{T})(\mathbf{\mu}-\mathbf{x})d\sigma\;,\;\;$
(20) $\displaystyle\tilde{\Lambda}_{\sigma}$
$\displaystyle=diag\big{[}\frac{\lambda_{k}}{\lambda_{k}+\sigma^{2}}\big{]}$
(21)

We can see it’s a linear ODE with respect to $\mathbf{x}-\mu$, with
disentangled dynamics along each eigen mode of $\mathbf{u}_{k}$. Choosing our
dynamic variable to be projection coefficient on each axes
$c_{k}(\sigma)=\mathbf{u}_{k}^{T}(\mathbf{x}-\mu)$. Then its dynamics could be
written as

$\displaystyle
dc_{k}(\sigma)=\frac{\sigma}{\lambda_{k}+\sigma^{2}}c_{k}(\sigma)d\sigma$ (22)

Integrating this, we get

$\displaystyle\frac{dc_{k}(\sigma)}{c_{k}(\sigma)}$
$\displaystyle=d\log\sqrt{\lambda_{k}+\sigma^{2}}$ (23) $\displaystyle d\log
c_{k}(\sigma)$ $\displaystyle=d\log\sqrt{\lambda_{k}+\sigma^{2}}$ (24)
$\displaystyle c_{k}(\sigma)$ $\displaystyle=\sqrt{\lambda_{k}+\sigma^{2}}K$
(25)

With integral constant $K$. Using the initial condition
$c_{k}(T)=\mathbf{u}_{k}^{T}(\mathbf{x}_{T}-\mu)$, then

$\displaystyle\frac{c_{k}(\sigma)}{c_{k}(T)}=\sqrt{\frac{\lambda_{k}+\sigma^{2}}{\lambda_{k}+\sigma_{T}^{2}}}=:\psi(\lambda_{k},\sigma)$
(26)

Thus, we can arrive at the solution Eq.4.

$\displaystyle\mathbf{x}_{\sigma}-\mu$
$\displaystyle=\sum_{k}c_{k}(\sigma)\mathbf{u}_{k}$ (27)
$\displaystyle\mathbf{x}_{\sigma}$
$\displaystyle=\mu+\sum_{k}\psi(\lambda_{k},\sigma)c_{k}(T)\mathbf{u}_{k}$
(28) $\displaystyle\mathbf{x}_{\sigma}$
$\displaystyle=\mu+\sum_{k}\sqrt{\frac{\lambda_{k}+\sigma^{2}}{\lambda_{k}+\sigma_{T}^{2}}}\mathbf{u}_{k}\mathbf{u}_{k}^{T}(\mathbf{x}_{T}-\mu)$
(29)

When the $\mathbf{U}$ is rank deficient, we will have the off-manifold term as
in Eq. 4.

For the denoising outcome, we have Eq. 2,

$\displaystyle D(\mathbf{x},\sigma)$
$\displaystyle=\mathbf{x}+\sigma^{2}\nabla\log p(\mathbf{x},\sigma)$ (30)
$\displaystyle=\mathbf{x}+(I-\mathbf{U}\tilde{\Lambda}_{\sigma}\mathbf{U}^{T})(\mathbf{\mu}-\mathbf{x})$
(31)
$\displaystyle=\mathbf{\mu}+\mathbf{U}\tilde{\Lambda}_{\sigma}\mathbf{U}^{T}(\mathbf{x}-\mathbf{\mu})$
(32)

Using the solution above, we have

$\displaystyle D(\mathbf{x}_{\sigma},\sigma)$
$\displaystyle=\mathbf{\mu}+\sum_{k}\frac{\lambda_{k}}{\sqrt{(\lambda_{k}+\sigma^{2})(\lambda_{k}+\sigma_{T}^{2})}}\mathbf{u}_{k}\mathbf{u}_{k}^{T}(\mathbf{x}_{T}-\mu)$
(33)

If we define

$\displaystyle\xi(\lambda_{k},\sigma):=\frac{\lambda_{k}}{\sqrt{(\lambda_{k}+\sigma^{2})(\lambda_{k}+\sigma_{T}^{2})}}$
(34)

and use the definition of $c_{k}(T)$, then we get the trajectory for denoiser
image as Eq.6.

$\displaystyle D(\mathbf{x}_{\sigma},\sigma)$
$\displaystyle=\mathbf{\mu}+\sum_{k}\xi(\lambda_{k},\sigma)c_{k}(T)\mathbf{u}_{k}$
(35)

###  9.2 Detailed derivation to the solution of the Gaussian Diffusion with
VP-SDE

As an extension, we present the solution for the VP-SDE forward process, which
has a data scaling term $\alpha_{t}$ for $\mathbf{x}$. As you will see, the
solutions are largely the same, just with an additional $\alpha_{t}$ term.

#### Forward/reverse process.

We consider forward processes defined by the stochastic differential equation
(SDE)

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

where $\beta(t)$ controls the decay of signal, $g(t)$ is a time-dependent
noise amplitude, $\mathbf{\eta}(t)$ is a vector of independent Gaussian white
noise terms, and time runs from $t=0$ to $T$. Its reverse process is

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

where $\mathbf{s}(\mathbf{x},t):=\nabla_{\mathbf{x}}\log p(\mathbf{x},t)$ is
the score function, and where we use the standard convention that time runs
backward, i.e. from $t=T$ to $0$. In this section, we focus on one popular
forward process: the variance-preserving SDE, which enforces the constraint
$\beta(t)=\frac{1}{2}g^{2}(t)$. The marginal probabilities of this process are

$\displaystyle
p(\mathbf{x}_{t}|\mathbf{x}_{0})=\mathcal{N}(\mathbf{x}_{t}|\alpha_{t}\mathbf{x}_{0},\sigma_{t}^{2}\mathbf{I})\hskip
21.68121pt$ (38)
$\displaystyle\alpha_{t}:=e^{-\int_{0}^{t}\beta(t^{\prime})dt^{\prime}}\hskip
21.68121pt\sigma_{t}^{2}:=1-e^{-2\int_{0}^{t}\beta(t^{\prime})dt^{\prime}}$
(39)

where $\alpha_{t}$ and $\sigma_{t}$ represent the signal and noise scale,
satisfying $\alpha_{t}^{2}+\sigma_{t}^{2}=1$. Normally, as $t$ goes from $0\to
T$, signal scale $\alpha_{t}$ monotonically decreases from $1\to 0$ and
$\sigma$ increases from $0\to 1$. Note that there exists a deterministic
probability flow ODE with the same marginal probabilities Song et al., (2021)
as the reverse SDE:

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

where time again runs backward from $t=1$ to $t=0$. In practice, instead of
the SDE (Eq.36), this deterministic process is often used to sample from the
distribution Karras et al., (2022). The behavior of the probability flow ODE
will be our main focus.

#### Gaussian Score model

For a Gaussian score model, the probability flow ODE that reverses a VP-SDE
forward process Song et al., (2021) is

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

Using the decomposition of $\mathbf{\Sigma}$ described above Woodbury, (1950),

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

where $\tilde{\mathbf{\Lambda}}_{t}$ is defined to be the time-dependent
diagonal matrix

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

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

$\displaystyle\frac{d}{dt}(\mathbf{x}_{t}-\alpha_{t}\mathbf{\mu})$
$\displaystyle=\dot{\mathbf{x}}_{t}-\mathbf{\mu}\dot{\alpha}_{t}$ (44)
$\displaystyle=\dot{\mathbf{x}}_{t}+\beta_{t}\alpha_{t}\mathbf{\mu}$ (45)
$\displaystyle=\beta_{t}(\alpha_{t}\mathbf{\mu}-\mathbf{x})-\frac{1}{2}g^{2}(t)\frac{1}{\sigma_{t}^{2}}(\mathbf{I}-\mathbf{U}\tilde{\mathbf{\Lambda}}_{t}\mathbf{U}^{T})(\alpha_{t}\mathbf{\mu}-\mathbf{x})$
(46)
$\displaystyle=\left[\frac{1}{2}g^{2}(t)\frac{1}{\sigma_{t}^{2}}(\mathbf{I}-\mathbf{U}\tilde{\mathbf{\Lambda}}_{t}\mathbf{U}^{T})-\beta_{t}\mathbf{I}\right](\mathbf{x}-\alpha_{t}\mathbf{\mu})\
.$ (47)

Assuming 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}(\mathbf{x}_{t}-\alpha_{t}\mathbf{\mu})$
$\displaystyle=\beta_{t}\left[\frac{1}{\sigma_{t}^{2}}(\mathbf{I}-\mathbf{U}\tilde{\mathbf{\Lambda}}_{t}\mathbf{U}^{T})-\mathbf{I}\right](\mathbf{x}-\alpha_{t}\mathbf{\mu})$
(48)
$\displaystyle=\beta_{t}\left[(\frac{1}{\sigma_{t}^{2}}-1)\mathbf{I}-\frac{1}{\sigma_{t}^{2}}\mathbf{U}\tilde{\mathbf{\Lambda}}_{t}\mathbf{U}^{T})\right](\mathbf{x}-\alpha_{t}\mathbf{\mu})\
.$ (49)

Define the variable $\mathbf{y}_{t}:=\mathbf{x}_{t}-\alpha_{t}\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 $\mathbf{y}_{t}$ in terms of the orthonormal columns of $\mathbf{U}$ and
a component that lies entirely in the orthogonal space $\mathbf{U}^{\perp}$:

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

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

$\displaystyle\dot{c}_{k}(t)=\frac{d}{dt}(\mathbf{u}_{k}^{T}\mathbf{y}_{t})$
$\displaystyle=\beta_{t}\left[(\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}}\right](\mathbf{u}_{k}^{T}\mathbf{y}_{t})$
(51)
$\displaystyle=\frac{\beta_{t}}{\sigma_{t}^{2}}\left(1-\sigma_{t}^{2}-\frac{\alpha_{t}^{2}\lambda_{k}}{\alpha_{t}^{2}\lambda_{k}+\sigma_{t}^{2}}\right)c_{k}(t)$
(52)
$\displaystyle=\frac{\beta_{t}\alpha_{t}^{2}}{\sigma_{t}^{2}}\left(1-\frac{\lambda_{k}}{\alpha_{t}^{2}\lambda_{k}+\sigma_{t}^{2}}\right)c_{k}(t)\
.$ (53)

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)\
.$ (54)

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

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

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

$\displaystyle\mathbf{y}_{t}$
$\displaystyle=d(t)\mathbf{y}^{\perp}(T)+\sum_{k=1}^{r}c_{k}(t)\mathbf{u}_{k}$
(56) $\displaystyle d(t)$
$\displaystyle=\exp\left(\int_{T}^{t}d\tau\frac{\beta_{\tau}\alpha_{\tau}^{2}}{\sigma_{\tau}^{2}}\right)$
(57) $\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)\
.$ (58)

The initial conditions are

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

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)}$
(61)
$\displaystyle=\int_{\alpha_{T}}^{\alpha_{t}}d\alpha_{\tau}\frac{\alpha_{\tau}(\lambda_{k}-1)}{1+\alpha_{\tau}^{2}(\lambda_{k}-1)}$
(62)
$\displaystyle=\frac{1}{2}\log(1+\alpha_{\tau}^{2}(\lambda_{k}-1))\Bigr{|}^{\alpha_{t}}_{\alpha_{T}}$
(63)
$\displaystyle=\frac{1}{2}\log\left(\frac{1+(\lambda_{k}-1)\alpha_{t}^{2}}{1+(\lambda_{k}-1)\alpha_{T}^{2}}\right)\
.$ (64)

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}}}=c_{k}(T)\sqrt{\frac{\sigma_{t}^{2}+\lambda_{k}\alpha_{t}^{2}}{\sigma_{T}^{2}+\lambda_{k}\alpha_{T}^{2}}}$
(65) $\displaystyle d(t)$
$\displaystyle=\sqrt{\frac{1-\alpha_{t}^{2}}{1-\alpha_{T}^{2}}}=\frac{\sigma_{t}}{\sigma_{T}}\
.$ (66)

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})}}}$
(67) $\displaystyle\dot{d}(t)$
$\displaystyle=\frac{\alpha_{t}^{2}\beta_{t}}{\sqrt{(1-\alpha_{T}^{2})(1-\alpha_{t}^{2})}}\
.$ (68)

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

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

The solution is a sum of on- and off-manifold components:

$\begin{split}\mathbf{x}_{t}=\alpha_{t}\mathbf{\mu}+\frac{\sigma_{t}}{\sigma_{T}}\
\mathbf{y}^{\perp}_{T}+\sum_{k=1}^{r}\psi(t,\lambda_{k})c_{k}(T)\mathbf{u}_{k}\hskip
36.135pt\psi(t,\lambda_{k})=\sqrt{\frac{\sigma_{t}^{2}+\lambda_{k}\alpha_{t}^{2}}{\sigma_{T}^{2}+\lambda_{k}\alpha_{T}^{2}}}\\\
\mathbf{y}^{\perp}_{T}=\
(\mathbf{I}-\mathbf{U}^{T}\mathbf{U})(\mathbf{x}_{T}-\alpha_{T}\mathbf{\mu})\hskip
36.135ptc_{k}(T)=\ \mathbf{u}_{k}^{T}(\mathbf{x}_{T}-\alpha_{T}\mathbf{\mu})\
.\end{split}$ (70)

The initial condition $\mathbf{x}_{T}$ can be decomposed into contributions
along each principal direction, an off-manifold contribution, and a
$\alpha_{T}\mathbf{\mu}$ contribution. There are three terms in the dynamics
for these three parts: 1) the scaling up of the distribution mean; 2) the
scaling down of the off-manifold component $\mathbf{y}^{\perp}_{T}$,
proportional to the noise scale $\psi(t,0)=\frac{\sigma_{t}}{\sigma_{T}}$; and
3) the on-manifold movement along each eigenvector governed by
$\psi(t,\lambda_{k})$.

We also now have the analytical solution for the projected outcome or the
ideal denoiser, which is defined as:

$D(\mathbf{x}_{t},\sigma_{t}):=\frac{\mathbf{x}_{t}+\sigma_{t}^{2}\mathbf{s}(\mathbf{x}_{t},t)}{\alpha_{t}}$
(71)

$\displaystyle D(\mathbf{x}_{t},\sigma_{t})-\mathbf{\mu}$
$\displaystyle=\frac{1}{\alpha_{t}}\mathbf{U}\tilde{\mathbf{\Lambda}}_{t}\mathbf{U}^{T}(\mathbf{x}_{t}-\alpha_{t}\mathbf{\mu})$
(72)
$\displaystyle=\sum_{k=1}^{r}c_{k}(t)\frac{\alpha_{t}\lambda_{k}}{\alpha_{t}^{2}\lambda_{k}+\sigma_{t}^{2}}\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})}}\mathbf{u}_{k}\
.$

Thus,

$\displaystyle
D(\mathbf{x}_{t},\sigma_{t})=\mathbf{\mu}+\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})}}\mathbf{u}_{k}$
(73)

###  9.3 Score function of general Gaussian mixture models

Let

$q(\mathbf{x})=\sum_{i}\pi_{i}\mathcal{N}(\mathbf{x};\mathbf{\mu}_{i},\mathbf{\Sigma}_{i})$
(74)

be a Gaussian mixture distribution, where the $\pi_{i}$ are mixture weights,
$\mathbf{\mu}_{i}$ is the $i$-th mean, and $\mathbf{\Sigma}_{i}$ is the $i$-th
covariance matrix. The score function for this distribution is

$\begin{split}\nabla_{\mathbf{x}}\log
q(\mathbf{x})=&\frac{\sum_{i}\pi_{i}\nabla_{\mathbf{x}}\mathcal{N}(\mathbf{x};\mathbf{\mu}_{i},\mathbf{\Sigma}_{i})}{q(\mathbf{x})}\\\
=&\sum_{i}-\Sigma_{i}^{-1}(\mathbf{x}-\mathbf{\mu}_{i})\frac{\pi_{i}\mathcal{N}(\mathbf{x};\mathbf{\mu}_{i},\mathbf{\Sigma}_{i})}{q(\mathbf{x})}\\\
=&\sum_{i}\frac{\pi_{i}\mathcal{N}(\mathbf{x};\mathbf{\mu}_{i},\mathbf{\Sigma}_{i})}{q(\mathbf{x})}\nabla_{\mathbf{x}}\log\mathcal{N}(\mathbf{x};\mathbf{\mu}_{i},\mathbf{\Sigma}_{i})\\\
=&\sum_{i}w_{i}(\mathbf{x})\nabla_{\mathbf{x}}\log\mathcal{N}(\mathbf{x};\mathbf{\mu}_{i},\mathbf{\Sigma}_{i})\end{split}$
(75)

where we have defined the mixing weights

$\displaystyle
w_{i}(\mathbf{x}):=\frac{\pi_{i}\mathcal{N}(\mathbf{x};\mathbf{\mu}_{i},\mathbf{\Sigma}_{i})}{q(\mathbf{x})}\
.$ (76)

Thus, the score of the Gaussian mixture is a weighted mixture of the score
fields of each of the individual Gaussians.

In the context of diffusion, we are interested in the time-dependent score
function. Given a Gaussian mixture initial condition, the end result of the
VP-SDE forward process will also be a Gaussian mixture:

$\displaystyle
p_{t}(\mathbf{x})=\sum_{i}\pi_{i}\mathcal{N}(\mathbf{x};\mathbf{\mu}_{i},\sigma_{t}^{2}\mathbf{I}+\mathbf{\Sigma}_{i})\
.$ (77)

The corresponding time-dependent score is

$\begin{split}\mathbf{s}(\mathbf{x},t)&=\nabla_{\mathbf{x}}\log
p_{t}(\mathbf{x})\\\
=&\sum_{i}-(\sigma_{t}^{2}\mathbf{I}+\mathbf{\Sigma}_{i})^{-1}(\mathbf{x}-\mathbf{\mu}_{i})\frac{\pi_{i}\mathcal{N}(\mathbf{x};\mathbf{\mu}_{i},\sigma_{t}^{2}\mathbf{I}+\mathbf{\Sigma}_{i})}{p_{t}(\mathbf{x})}\\\
=&\sum_{i}-(\sigma_{t}^{2}\mathbf{I}+\mathbf{\Sigma}_{i})^{-1}(\mathbf{x}-\mathbf{\mu}_{i})w_{i}(\mathbf{x},t)\
.\end{split}$ (78)

Note that we have a formula for
$(\sigma_{t}^{2}\mathbf{I}+\mathbf{\Sigma}_{i})^{-1}$ (derived using the
Woodbury matrix inversion identity; see Eq.3) in terms of the (compact) SVD of
$\mathbf{\Sigma}_{i}$. We can use it to write

$\displaystyle\mathbf{s}(\mathbf{x},t)=\frac{1}{\sigma_{t}^{2}}\sum_{i}-(\mathbf{I}-\mathbf{U}_{i}\tilde{\mathbf{\Lambda}_{i}}_{t}\mathbf{U}_{i}^{T})(\mathbf{x}-\mathbf{\mu}_{i})w_{i}(\mathbf{x},t)$
(79)

where

$\mathbf{\Sigma}_{i}=\mathbf{U}_{i}\mathbf{\Lambda}_{i}\mathbf{U}_{i}^{T}\hskip
72.26999pt\tilde{\mathbf{\Lambda}}_{t}:=\text{diag}\left[\frac{\lambda_{k}}{\lambda_{k}+\sigma_{t}^{2}}\right]\
.$ (80)

This representation of the score function is numerically convenient, since
(once the SVDs of each covariance matrix have been obtained), it can be
evaluated using a relatively small number of matrix multiplications, which are
cheaper than the covariance matrix inversions that a naive implementation of
the Gaussian mixture score function would require.

We used this formula (Eq.79) and an off-the-shelf ODE solver to simulate the
reverse diffusion trajectory of a 10-mode Gaussian mixture score model (Fig.
7).

###  9.4 Score function of Gaussian mixture with identical and isotropic
covariance

Assume each Gaussian mode has covariance
$\mathbf{\Sigma}_{i}=\sigma^{2}\mathbf{I}$ and that every mixture weight is
the same (i.e. $\pi_{i}=\pi_{j},\forall i,j$). Then the score for this kind of
specific Gaussian mixture is

$\begin{split}\nabla_{\mathbf{x}}\log
q(\mathbf{x})=&\frac{\sum_{i}\pi_{i}\nabla_{\mathbf{x}}\mathcal{N}(\mathbf{x};\mathbf{\mu}_{i},\sigma^{2}\mathbf{I})}{\sum_{i}\pi_{i}\mathcal{N}(\mathbf{x};\mathbf{\mu}_{i},\sigma^{2}\mathbf{I})}\\\
=&\frac{\sum_{i}-\frac{1}{\sigma^{2}}(\mathbf{x}-\mathbf{\mu}_{i})\exp\big{(}-\frac{1}{2\sigma^{2}}\|\mathbf{x}-\mathbf{\mu}_{i}\|_{2}^{2}\big{)}}{\sum_{i}\exp\big{(}-\frac{1}{2\sigma^{2}}\|\mathbf{x}-\mathbf{\mu}_{i}\|_{2}^{2}\big{)}}\\\
=&\frac{1}{\sigma^{2}}\sum_{i}w_{i}(\mathbf{x})(\mathbf{\mu}_{i}-\mathbf{x})\
,\end{split}$ (81)

where the weight $w_{i}(\mathbf{x})$ is a softmax of the negative squared
distance to all the means, with $\sigma^{2}$ functioning as a temperature
parameter:

$\displaystyle w_{i}(\mathbf{x})$
$\displaystyle=Softmax(\big{\\{}-\frac{1}{2\sigma^{2}}\|\mathbf{x}-\mathbf{\mu}_{i}\|_{2}^{2}\big{\\}})=\frac{\exp\big{(}-\frac{1}{2\sigma^{2}}\|\mathbf{x}-\mathbf{\mu}_{i}\|_{2}^{2}\big{)}}{\sum_{i}\exp\big{(}-\frac{1}{2\sigma^{2}}\|\mathbf{x}-\mathbf{\mu}_{i}\|_{2}^{2}\big{)}}\
.$ (82)

Since the weights $w_{i}(\mathbf{x})$ sum to $1$, we can also write the score
function in the suggestive form

$\displaystyle\nabla_{\mathbf{x}}\log
q(\mathbf{x})=\frac{(\sum_{i}w_{i}(\mathbf{x})\mathbf{\mu}_{i})-\mathbf{x}}{\sigma^{2}}\
.$ (83)

This has a form analogous to the score of a single Gaussian mode—but instead
of $\mathbf{x}$ being ‘attracted’ towards a single mean $\mathbf{\mu}$, it is
attracted towards a weighted combination of all of the means, with modes
closer to the state $\mathbf{x}$ being more highly weighted.

###  9.5 Score of exact (delta mixture) score model

A particularly interesting special case of the Gaussian mixture model is the
delta mixture model used in the main text, whose components are vanishing-
width Gaussians centered on the training images. In particular, consider a
data set $\\{\mathbf{y}_{i}\\}$ with $i=1,...,N$, so that the starting
distribution is

$p(\mathbf{x})=\frac{1}{N}\sum_{i}\delta(\mathbf{x}-\mathbf{y}_{i})\ .$ (84)

At time $t$, the marginal distribution will be a Gaussian mixture

$\displaystyle
p_{t}(\mathbf{x}_{t})=\frac{1}{N}\sum_{i}\mathcal{N}(\mathbf{x}_{t};\mathbf{y}_{i},\sigma_{t}^{2}\mathbf{I})\
.$ (85)

Then using the Eq. 81 above we have

$\begin{split}s(\mathbf{x}_{t},t)=\nabla\log
p_{t}(\mathbf{x}_{t})&=\frac{1}{\sigma_{t}^{2}}\sum_{i}w_{i}(\mathbf{x}_{t})(\mathbf{y}_{i}-\mathbf{x}_{t})\\\
&=\frac{1}{\sigma_{t}^{2}}\left[-\mathbf{x}_{t}+\sum_{i}w_{i}(\mathbf{x}_{t})\mathbf{y}_{i}\right]\\\
&=\frac{1}{\sigma_{t}^{2}}\left[-\mathbf{x}_{t}+\sum_{i}Softmax(\big{\\{}-\frac{1}{2\sigma_{t}^{2}}\|\mathbf{y}_{i}-\mathbf{x}_{t}\|^{2}\big{\\}})\mathbf{y}_{i}\right]\
.\end{split}$ (86)

The endpoint estimate of the distribution is

$\begin{split}\hat{\mathbf{x}}_{0}(\mathbf{x}_{t})&=\frac{\mathbf{x}_{t}+\sigma_{t}^{2}\nabla\log
p(\mathbf{x}_{t})}{}\\\
&=\sum_{i}Softmax(\big{\\{}-\frac{1}{2\sigma_{t}^{2}}\|\mathbf{y}_{i}-\mathbf{x}_{t}\|^{2}\big{\\}})\mathbf{y}_{i}\\\
&=\sum_{i}w_{i}(\mathbf{x}_{t})\mathbf{y}_{i}\ .\end{split}$ (87)

Thus, the endpoint estimate is a weighted average of training data, with the
softmax of negative squared distance as weights and $\sigma_{t}^{2}$ as a
temperature parameter.

###  9.6 Arguments of the approximation of score field of Gaussian and point
cloud

In this section, we will argue from a theoretical perspective, the score field
of a bounded point cloud is equivalent to a Gaussian with matching mean and
covariance, when the noise level is high and when the query point is far from
the bounded point cloud. We are going to use a technique inspired by multi-
pole expansion in electrodynamics (Griffiths,, 2005).

#### Set up

Assume we have a set of data points $\\{\mathbf{y}_{i}\\},i=1...N$, their mean
and covariance are defined as

$\displaystyle\mu$ $\displaystyle=\frac{1}{N}\sum_{i}\mathbf{y}_{i}$ (88)
$\displaystyle\Sigma$
$\displaystyle=\frac{1}{N}\sum_{i}\mathbf{y}_{i}\mathbf{y}_{i}^{T}-\mu\mu^{T}$
(89)

these are the first moment and the second central moment of the distribution.
The Gaussian distribution with matching first two moments is
$\mathcal{N}(\mu,\Sigma)$. We are going to show that when the noise level
$\sigma$ and $\|\mathbf{x}-\mu\|$ are both much larger than the standard
deviation of the distribution, the score at noise level $\sigma$ of the
original dataset (point cloud) is equivalent to the score of this Gaussian
distribution.

#### Score expansion for Gaussian distribution

Consider the Gaussian distribution, at noise level $\sigma$, the distribution
is $\mathcal{N}(\mu,\Sigma+\sigma^{2}I)$. Then its score can be written as

$\displaystyle\log p(\mathbf{x};\sigma)$
$\displaystyle=(\sigma^{2}I+\Sigma)^{-1}(\mathbf{\mu}-\mathbf{x})$ (90)
$\displaystyle=\frac{1}{\sigma^{2}}(I-U\tilde{\Lambda}_{\sigma}U^{T})(\mathbf{\mu}-\mathbf{x})$
(91)

in which the

$\displaystyle\tilde{\Lambda}_{\sigma}$
$\displaystyle=diag\big{[}\frac{\lambda_{k}}{\lambda_{k}+\sigma^{2}}\big{]}$
(92)
$\displaystyle=diag\big{[}\frac{\lambda_{k}}{\sigma^{2}}-\frac{\lambda_{k}^{2}}{\sigma^{2}(\sigma^{2}+\lambda_{k})}\big{]}$
(93)
$\displaystyle=diag\big{[}\frac{\lambda_{k}}{\sigma^{2}}-\frac{\lambda_{k}^{2}}{\sigma^{4}}+\frac{\lambda_{k}^{3}}{\sigma^{4}(\sigma^{2}+\lambda_{k})}\big{]}$
(94)
$\displaystyle=\frac{1}{\sigma^{2}}\Lambda-\frac{1}{\sigma^{4}}\Lambda^{2}+\frac{1}{\sigma^{6}}\Lambda^{3}-\frac{1}{\sigma^{8}}\Lambda^{4}+...$
(95) $\displaystyle=\frac{1}{\sigma^{2}}\Lambda-\Delta$ (96)

where the residual term $\Delta$ is the 2nd order term of
$\lambda/\sigma^{2}$, $\Delta\sim(\lambda/\sigma^{2})^{2}$.

Using this expansion, the score field can be expressed as

$\displaystyle\log p(\mathbf{x};\sigma)$
$\displaystyle=\frac{1}{\sigma^{2}}(I-\frac{1}{\sigma^{2}}U\Lambda
U^{T}+U\Delta U^{T})(\mathbf{\mu}-\mathbf{x})$ (97)
$\displaystyle=\frac{1}{\sigma^{2}}(\mathbf{\mu}-\mathbf{x})-\frac{1}{\sigma^{4}}\Sigma(\mathbf{\mu}-\mathbf{x})+\frac{1}{\sigma^{2}}U\Delta
U^{T}(\mathbf{\mu}-\mathbf{x})$ (98)
$\displaystyle=\frac{1}{\sigma^{2}}(\mathbf{\mu}-\mathbf{x})-\frac{1}{\sigma^{4}}\Sigma(\mathbf{\mu}-\mathbf{x})+\frac{1}{\sigma^{6}}U\Lambda^{2}U^{T}(\mathbf{\mu}-\mathbf{x})-\frac{1}{\sigma^{8}}U\Lambda^{3}U^{T}(\mathbf{\mu}-\mathbf{x})+...$
(99)
$\displaystyle=\frac{1}{\sigma^{2}}(\mathbf{\mu}-\mathbf{x})-\frac{1}{\sigma^{4}}\Sigma(\mathbf{\mu}-\mathbf{x})+\frac{1}{\sigma^{6}}\Sigma^{2}(\mathbf{\mu}-\mathbf{x})-\frac{1}{\sigma^{8}}\Sigma^{3}(\mathbf{\mu}-\mathbf{x})+...$
(100)

From this result, we can express the score field of a Gaussian distribution by
a series of terms: the first-order term is the isotropic attraction to the
mean, the second-order term is the anisotropic term conditioned by the
covariance matrix; the higher-order terms are increasingly anisotropic with
higher power of the covariance matrix, i.e. larger condition number.

Note that, this series is always linear with respect to
$\mathbf{\mu}-\mathbf{x}$, and it doesn’t contain any quadratic term of
$\mathbf{x}$.

#### Score expansion for delta point distribution

Consider the data points, at noise level $\sigma$, the noised distribution is
a Gaussian mixture with $N$ components,
$q(\mathbf{x};\sigma)=\frac{1}{N}\sum_{i}\mathcal{N}(\mathbf{x};\mathbf{y}_{i},\sigma^{2}I)$.
As we have derived above, the score of a Gaussian mixture with identical
variance at each mode can be expressed as

$\displaystyle\nabla_{\mathbf{x}}\log q(\mathbf{x};\sigma)$
$\displaystyle=\frac{1}{\sigma^{2}}\sum_{i}w_{i}(\mathbf{x})(\mathbf{y}_{i}-\mathbf{x})$
(101)
$\displaystyle=\frac{1}{\sigma^{2}}\bigg{(}\sum_{i}w_{i}(\mathbf{x})\mathbf{y}_{i}-\mathbf{x}\bigg{)}$
(102)
$\displaystyle=\frac{1}{\sigma^{2}}(\mathbf{\mu}-\mathbf{x})+\frac{1}{\sigma^{2}}\sum_{i}w_{i}(\mathbf{x})(\mathbf{y}_{i}-\mathbf{\mu})$
(103)

Where the weighting function is

$\displaystyle
w_{i}(\mathbf{x})=\frac{\exp\big{(}-\frac{1}{2\sigma^{2}}\|\mathbf{x}-\mathbf{y}_{i}\|^{2}_{2}\big{)}}{\sum^{N}_{j}\exp\big{(}-\frac{1}{2\sigma^{2}}\|\mathbf{x}-\mathbf{y}_{j}\|^{2}_{2}\big{)}}$
(104)

Here we rewrite all the distances using the distributional mean as a reference

$\displaystyle\|\mathbf{x}-\mathbf{y}_{i}\|^{2}_{2}$
$\displaystyle=\|\mathbf{x}-\mu+\mu-\mathbf{y}_{i}\|^{2}_{2}$ (105)
$\displaystyle=\|\mathbf{x}-\mu\|^{2}_{2}+\|\mu-\mathbf{y}_{i}\|^{2}_{2}+2(\mathbf{x}-\mu)^{T}(\mu-\mathbf{y}_{i})$
(106)

Multiply the same term
$\exp\big{(}\frac{1}{2\sigma^{2}}\|\mathbf{x}-\mu\|^{2}_{2}\big{)}$ at
numerator and denominator, we get

$\displaystyle
w_{i}(\mathbf{x})=\frac{\exp\big{(}-\frac{1}{2\sigma^{2}}(\|\mu-\mathbf{y}_{i}\|^{2}_{2}+2(\mathbf{x}-\mu)^{T}(\mu-\mathbf{y}_{i}))\big{)}}{\sum^{N}_{j}\exp\big{(}-\frac{1}{2\sigma^{2}}(\|\mu-\mathbf{y}_{j}\|^{2}_{2}+2(\mathbf{x}-\mu)^{T}(\mu-\mathbf{y}_{j}))\big{)}}$
(107)

Note the $(\mathbf{x}-\mu)^{T}(\mu-\mathbf{y}_{i})$ is on the order of
$\|\mathbf{x}-\mu\|\|\mu-\mathbf{y}_{i}\|\propto\sigma\sqrt{d}\sqrt{tr\Lambda}$.
$\|\mu-\mathbf{y}_{j}\|^{2}_{2}$ is on the order of $tr\Lambda$.
$\|\mathbf{x}-\mu\|^{2}$ is on the order of $d\sigma^{2}$.

We can expand the exponential functions to the first-order of
$\sqrt{d}\sqrt{tr\Lambda}/\sigma$

$\displaystyle w_{i}(\mathbf{x})\approx$ $\displaystyle\
\frac{1-\frac{1}{\sigma^{2}}(\mathbf{x}-\mu)^{T}(\mu-\mathbf{y}_{i})-\frac{1}{2\sigma^{2}}\|\mu-\mathbf{y}_{i}\|^{2}_{2}}{N+\sum^{N}_{j}\big{[}-\frac{1}{\sigma^{2}}(\mathbf{x}-\mu)^{T}(\mu-\mathbf{y}_{j})-\frac{1}{2\sigma^{2}}\|\mu-\mathbf{y}_{j}\|^{2}_{2}\big{]}}$
(108) $\displaystyle=$ $\displaystyle\
\frac{1-\frac{1}{\sigma^{2}}(\mathbf{x}-\mu)^{T}(\mu-\mathbf{y}_{i})-\frac{1}{2\sigma^{2}}\|\mu-\mathbf{y}_{i}\|^{2}_{2}}{N+\sum^{N}_{j}\big{[}-\frac{1}{2\sigma^{2}}\|\mu-\mathbf{y}_{j}\|^{2}_{2}\big{]}}$
(109) $\displaystyle=$ $\displaystyle\
\frac{1-\frac{1}{\sigma^{2}}(\mathbf{x}-\mu)^{T}(\mu-\mathbf{y}_{i})-\frac{1}{2\sigma^{2}}\|\mu-\mathbf{y}_{i}\|^{2}_{2}}{N-N\frac{tr\Sigma}{2\sigma^{2}}}$
(110)

Here we used the relationship

$\displaystyle\sum^{N}_{j}(\mu-\mathbf{y}_{j})=0$ (111)
$\displaystyle\sum^{N}_{j}\|\mu-\mathbf{y}_{j}\|^{2}_{2}=Ntr\Sigma$ (112)

Now consider $\sum_{i}w_{i}(\mathbf{x})\mathbf{y}_{i}$

$\displaystyle\sum^{N}_{i}w_{i}(\mathbf{x})\mathbf{y}_{i}$
$\displaystyle\approx\frac{1}{N(1-tr\Sigma/2\sigma^{2})}\sum^{N}_{i}\bigg{[}1-\frac{1}{\sigma^{2}}(\mathbf{x}-\mu)^{T}(\mu-\mathbf{y}_{i})-\frac{1}{2\sigma^{2}}\|\mu-\mathbf{y}_{i}\|^{2}_{2}\bigg{]}\mathbf{y}_{i}$
(113)
$\displaystyle=\frac{1}{N(1-tr\Sigma/2\sigma^{2})}\bigg{[}N\mu-\frac{1}{\sigma^{2}}\sum^{N}_{i}(\mathbf{x}-\mu)^{T}(\mu-\mathbf{y}_{i})\mathbf{y}_{i}-\frac{1}{2\sigma^{2}}\sum^{N}_{i}\|\mu-\mathbf{y}_{i}\|^{2}_{2}\mathbf{y}_{i}\bigg{]}$
(114)

Here we used the relationship

$\displaystyle\sum^{N}_{i}(\mathbf{x}-\mu)^{T}(\mu-\mathbf{y}_{i})\mathbf{y}_{i}$
$\displaystyle=\sum^{N}_{i}(\mu-\mathbf{y}_{i})^{T}(\mathbf{x}-\mu)\mathbf{y}_{i}$
(115)
$\displaystyle=\sum^{N}_{i}\mathbf{y}_{i}(\mu-\mathbf{y}_{i})^{T}(\mathbf{x}-\mu)$
(116)
$\displaystyle=\sum^{N}_{i}(\mathbf{y}_{i}-\mu)(\mu-\mathbf{y}_{i})^{T}(\mathbf{x}-\mu)$
(117) $\displaystyle=-N\Sigma(\mathbf{x}-\mu)$ (118)

and

$\displaystyle\sum^{N}_{i}\|\mu-\mathbf{y}_{i}\|^{2}_{2}\mathbf{y}_{i}$
$\displaystyle=\sum^{N}_{i}\|\mu-\mathbf{y}_{i}\|^{2}_{2}(\mathbf{y}_{i}-\mu)+\mu\sum^{N}_{i}\|\mu-\mathbf{y}_{i}\|^{2}_{2}$
(119)
$\displaystyle=\sum^{N}_{i}\|\mathbf{y}_{i}-\mu\|^{2}_{2}(\mathbf{y}_{i}-\mu)+\mu
Ntr\Sigma$ (120) $\displaystyle=N\mathbf{\gamma}+Ntr\Sigma\mu$ (121)

where
$\mathbf{\gamma}:=1/N\sum^{N}_{i}\|\mathbf{y}_{i}-\mu\|^{2}_{2}(\mathbf{y}_{i}-\mu)$
is a form of third centered moment, which measures the assymetry of the
distribution with respect to the mean.

Take these Eqs. 118,121 into Eq. 114, we have

$\displaystyle\sum^{N}_{i}w_{i}(\mathbf{x})\mathbf{y}_{i}$
$\displaystyle\approx\frac{1}{N(1-tr\Sigma/2\sigma^{2})}\bigg{[}N\mu+\frac{1}{\sigma^{2}}N\Sigma(\mathbf{x}-\mu)-\frac{1}{2\sigma^{2}}\big{(}N\mathbf{\gamma}+Ntr\Sigma\mu\big{)}\bigg{]}$
(122)
$\displaystyle=\mu+\frac{1}{\sigma^{2}(1-tr\Sigma/2\sigma^{2})}\Sigma(\mathbf{x}-\mu)-\frac{1}{2\sigma^{2}(1-tr\Sigma/2\sigma^{2})}\mathbf{\gamma}$
(123)

Taking together, the score function reads,

$\displaystyle\nabla_{\mathbf{x}}\log q(\mathbf{x};\sigma)$
$\displaystyle=\frac{1}{\sigma^{2}}(\sum_{i}w_{i}(\mathbf{x})\mathbf{y}_{i}-\mathbf{x})$
(124)
$\displaystyle\approx\frac{1}{\sigma^{2}}(\mu-\mathbf{x})+\frac{1}{\sigma^{4}}\Sigma(\mathbf{x}-\mu)-\frac{1}{2\sigma^{4}}\gamma+...$
(125)

We can see the expansion of the score for point cloud (Eq. 125) share the
first two terms with the expansion of the Gaussian score (Eq. 100): the
isotropic term of $(\mu-\mathbf{x})/\sigma^{2}$ and the first order
anisotropic term determined by the covariance.

#### Alternative score expansion for delta point distribution

Another way of expansion is the following

$\displaystyle
w_{i}(\mathbf{x})=\frac{\exp\big{(}-\frac{1}{2\sigma^{2}}(\|\mu-\mathbf{y}_{i}\|^{2}_{2}+2(\mathbf{x}-\mu)^{T}(\mu-\mathbf{y}_{i}))-tr\Sigma\big{)}}{\sum^{N}_{j}\exp\big{(}-\frac{1}{2\sigma^{2}}(\|\mu-\mathbf{y}_{j}\|^{2}_{2}+2(\mathbf{x}-\mu)^{T}(\mu-\mathbf{y}_{j})-tr\Sigma)\big{)}}$
(126)

Note the $(\mathbf{x}-\mu)^{T}(\mu-\mathbf{y}_{i})$ is on the order of
$\|\mathbf{x}-\mu\|\|\mu-\mathbf{y}_{i}\|\propto\sigma\sqrt{d}\sqrt{tr\Lambda}$.
$\|\mu-\mathbf{y}_{j}\|^{2}_{2}$ is on the order of $tr\Lambda$.
$\|\mathbf{x}-\mu\|^{2}$ is on the order of $d\sigma^{2}$.

We can expand the exponential functions to the first-order of
$\sqrt{d}\sqrt{tr\Lambda}/\sigma$

$\displaystyle w_{i}(\mathbf{x})\approx$ $\displaystyle\
\frac{1-\frac{1}{\sigma^{2}}(\mathbf{x}-\mu)^{T}(\mu-\mathbf{y}_{i})-\frac{1}{2\sigma^{2}}\|\mu-\mathbf{y}_{i}\|^{2}_{2}+\frac{1}{2\sigma^{2}}tr\Sigma}{N+\sum^{N}_{j}\big{[}-\frac{1}{\sigma^{2}}(\mathbf{x}-\mu)^{T}(\mu-\mathbf{y}_{j})-\frac{1}{2\sigma^{2}}\|\mu-\mathbf{y}_{j}\|^{2}_{2}+\frac{1}{2\sigma^{2}}tr\Sigma\big{]}}$
(127) $\displaystyle=$ $\displaystyle\
\frac{1-\frac{1}{\sigma^{2}}(\mathbf{x}-\mu)^{T}(\mu-\mathbf{y}_{i})-\frac{1}{2\sigma^{2}}\|\mu-\mathbf{y}_{i}\|^{2}_{2}+\frac{1}{2\sigma^{2}}tr\Sigma}{N+\sum^{N}_{j}\big{[}-\frac{1}{2\sigma^{2}}\|\mu-\mathbf{y}_{j}\|^{2}_{2}+\frac{N}{2\sigma^{2}}tr\Sigma\big{]}}$
(128) $\displaystyle=$ $\displaystyle\
\frac{1-\frac{1}{\sigma^{2}}(\mathbf{x}-\mu)^{T}(\mu-\mathbf{y}_{i})-\frac{1}{2\sigma^{2}}\|\mu-\mathbf{y}_{i}\|^{2}_{2}+\frac{1}{2\sigma^{2}}tr\Sigma}{N}$
(129) $\displaystyle=$
$\displaystyle\frac{1}{N\sigma^{2}}\big{[}\sigma^{2}-(\mathbf{x}-\mu)^{T}(\mu-\mathbf{y}_{i})-\frac{1}{2}\|\mu-\mathbf{y}_{i}\|^{2}_{2}+\frac{1}{2}tr\Sigma\big{]}$
(130)

Here we used the identity Eq.112 to simplify the denominator.

$\displaystyle\sum_{i}w_{i}(\mathbf{x})(\mathbf{y}_{i}-\mathbf{\mu})$
$\displaystyle\approx\frac{1}{N\sigma^{2}}\big{[}N\Sigma(\mathbf{x}-\mu)-\sum_{i}\frac{1}{2}\|\mu-\mathbf{y}_{i}\|^{2}(\mathbf{y}_{i}-\mathbf{\mu})\big{]}$
(131)
$\displaystyle=\frac{1}{\sigma^{2}}\Sigma(\mathbf{x}-\mu)-\frac{1}{2\sigma^{2}}\sum_{i}\|\mu-\mathbf{y}_{i}\|^{2}(\mathbf{y}_{i}-\mathbf{\mu})$
(132)

Thus

$\displaystyle\nabla_{\mathbf{x}}\log q(\mathbf{x};\sigma)$
$\displaystyle=\frac{1}{\sigma^{2}}(\mathbf{\mu}-\mathbf{x})+\frac{1}{\sigma^{2}}\sum_{i}w_{i}(\mathbf{x})(\mathbf{y}_{i}-\mathbf{\mu})$
(133)
$\displaystyle\approx\frac{1}{\sigma^{2}}(\mathbf{\mu}-\mathbf{x})+\frac{1}{\sigma^{4}}\Sigma(\mathbf{x}-\mu)-\frac{1}{2\sigma^{4}}\sum_{i}\|\mu-\mathbf{y}_{i}\|^{2}(\mathbf{y}_{i}-\mathbf{\mu})+...$
(134)

Similarly, using this expansion (Eq. 134), we can see the expansion of the
score for point cloud share the first two terms with the expansion of the
Gaussian score (Eq. 100): the isotropic term of $(\mu-\mathbf{x})/\sigma^{2}$
and the first order anisotropic term determined by the covariance.

#### Validity of Assumption

For a $d$ dimensional Gaussian distribution
$\mathbf{x}\sim\mathcal{N}(0,\sigma^{2}I)$, the squared norm
$\|\mathbf{x}/\sigma\|^{2}\sim\chi^{2}(d)$,
$\mathbb{E}\|\mathbf{x}\|^{2}=d\sigma^{2}$. As
$\sum_{i}^{N}\|\mathbf{y}_{i}-\mathbf{\mu}\|^{2}=Ntr\Lambda$,
$\|\mathbf{y}_{i}-\mathbf{\mu}\|^{2}$ is on the order of
$tr\Lambda=\sum_{k}\lambda_{k}$. For a bounded point cloud where
$\|\mu\|\ll\|\mathbf{x}\|$ and $\|\mu-\mathbf{y}_{i}\|\ll\|\mathbf{x}\|$.
Then, $\|\mathbf{x}-\mathbf{y}_{i}\|^{2}$ and $\|\mathbf{x}-\mu\|^{2}$ are
also on the order of $d\sigma^{2}$. Further works are needed to show the more
precise condition at which the Gaussian approximation holds well.