Score Matching

Score matching fits a model of the score of a distribution, the gradient of its log-density with respect to the data, $\nabla_x \log p(x)$, instead of the density itself. For images the score is a vector field that points, at every image, toward more probable images. A gradient with respect to $x$ ignores constant factors, so the score never involves the normalizing constant that makes flexible image models intractable. Hyvärinen (2005) introduced score matching; its denoising form (Vincent, 2011) trains score-based generative models and, after a change of variables, the diffusion models whose notation this article follows.

Definition

Score matching estimates a model score $s_\theta(x) \approx \nabla_x \log p(x)$ by minimizing the expected squared distance between the two vector fields under the data distribution $p$, the Fisher divergence (Song, Garg, Shi & Ermon, 2019). Since the data score is unknown, usable objectives remove it by integration by parts (implicit score matching), by noising the data (denoising score matching) or by random projections (sliced score matching). The classical score is a gradient with respect to parameters; Hyvärinen notes that his is the score for a hypothetical location parameter.

Intuition

A model $p_\theta(x) = \tilde p_\theta(x) / Z_\theta$ with $Z_\theta = \int \tilde p_\theta(x)\, dx$ is easy to write down and hard to normalize: for $8 \times 8$ patches $Z_\theta$ is a 64-dimensional integral, and maximum likelihood must resort to Markov chain Monte Carlo or approximations (Hyvärinen, 2005). The score sidesteps it, because $\log Z_\theta$ does not depend on $x$:

$$\nabla_x \log p_\theta(x) = \nabla_x \log \tilde p_\theta(x).$$

A score also suffices to generate, by Langevin dynamics, so score-based generative modeling learns the score and then samples with it (Song & Ermon, 2019).

Formal Definition

Explicit and implicit score matching

The explicit objective is

$$J_{\mathrm{ESM}}(\theta) = \tfrac{1}{2}\, \mathbb{E}_{p(x)} \big\lVert s_\theta(x) - \nabla_x \log p(x) \big\rVert^2 .$$

Hyvärinen (2005, Theorem 1) proved that, if $p$ and $s_\theta$ are differentiable, the relevant expectations are finite and $p(x)\, s_\theta(x) \to 0$ as $\lVert x \rVert \to \infty$,

$$J_{\mathrm{ESM}}(\theta) = \mathbb{E}_{p(x)} \Big[ \operatorname{tr}\!\big(\nabla_x s_\theta(x)\big) + \tfrac{1}{2} \lVert s_\theta(x) \rVert^2 \Big] + C = J_{\mathrm{ISM}}(\theta) + C,$$

with $C$ independent of $\theta$; Vincent (2011) named the two forms explicit (ESM) and implicit (ISM). Expanding the square, only the cross term $-\mathbb{E}[s_\theta \cdot \nabla \log p]$ involves the unknown score, and integration by parts along each coordinate turns it into $\mathbb{E}[\operatorname{tr} \nabla_x s_\theta]$. In one dimension,

$$\int p\, (\log p)'\, f\, dx = \int p'\, f\, dx = -\int p\, f'\, dx .$$

If the model family contains the data distribution, is identifiable and positive, the estimator is consistent, provided the optimizer finds the global minimum (Theorem 2, Corollary 3). For a Gaussian model it equals the maximum-likelihood estimate for any sample (§3.1), and Hyvärinen fitted an overcomplete model of 200 filters to $8 \times 8$ natural-image patches, whose normalizing constant is unknown (§3.3). For neural networks the trace is the obstacle: computed naively it costs $D$ times more backward passes than the score itself for a $D$-pixel image (Song et al., 2019, §2.1).

Denoising score matching

Vincent (2011) perturbed each training image, $\tilde x = x + \sigma z$ with $z \sim \mathcal{N}(0, I)$, and targeted the score of the noisy density $q_\sigma(\tilde x) = \int p(x)\, \mathcal{N}(\tilde x; x, \sigma^2 I)\, dx$ with the denoising score matching (DSM) objective

$$J_{\mathrm{DSM}}(\theta) = \tfrac{1}{2}\, \mathbb{E}_{x, \tilde x} \big\lVert s_\theta(\tilde x) - \nabla_{\tilde x} \log q_\sigma(\tilde x \mid x) \big\rVert^2, \qquad \nabla_{\tilde x} \log q_\sigma(\tilde x \mid x) = -\frac{\tilde x - x}{\sigma^2} = -\frac{z}{\sigma}.$$

The target needs no derivatives of the network. Vincent proved that $J_{\mathrm{DSM}}$ equals explicit score matching against $q_\sigma$ up to a constant, so its minimizer is the score of the noisy data, not of $p$. In our reading, the proof amounts to $\nabla \log q_\sigma(\tilde x) = \mathbb{E}[-(\tilde x - x)/\sigma^2 \mid \tilde x]$: the noisy score is the per-pair target averaged over the clean images that could have produced $\tilde x$. Hence the best denoiser is $\mathbb{E}[x \mid \tilde x] = \tilde x + \sigma^2 \nabla \log q_\sigma(\tilde x)$; denoising and score estimation are one problem.

In the diffusion notation, $x_t = \sqrt{\bar\alpha_t}\, x_0 + \sqrt{1 - \bar\alpha_t}\, \epsilon$ and the target is $\nabla_{x_t} \log q(x_t \mid x_0) = -\epsilon / \sqrt{1 - \bar\alpha_t}$. Weighting the squared error of the DSM loss at step $t$ by $1 - \bar\alpha_t$ and writing $s_\theta(x_t, t) = -\epsilon_\theta(x_t, t) / \sqrt{1 - \bar\alpha_t}$ gives

$$(1 - \bar\alpha_t) \Big\lVert s_\theta(x_t, t) + \frac{\epsilon}{\sqrt{1 - \bar\alpha_t}} \Big\rVert^2 = \big\lVert \epsilon - \epsilon_\theta(x_t, t) \big\rVert^2 ,$$

the noise-prediction loss of DDPM. Ho, Jain and Abbeel (2020, §3.2) observed that their loss resembles denoising score matching over noise levels indexed by $t$, and Song et al. (2021, Eq. 3) wrote it in this weighted form. In our reading, the $\epsilon$–score relation of the diffusion model article is thus a property of the training objective.

Sliced score matching

Song, Garg, Shi and Ermon (2019) kept the clean-data score but compared projections onto random directions $v$, such as standard normal vectors:

$$J_{\mathrm{SSM}}(\theta) = \mathbb{E}_{v}\, \mathbb{E}_{p(x)} \Big[ v^\top \nabla_x s_\theta(x)\, v + \tfrac{1}{2} \big(v^\top s_\theta(x)\big)^2 \Big].$$

$v^\top \nabla_x s_\theta\, v$ costs one extra gradient computation per direction rather than one per pixel; they proved the estimator consistent and asymptotically normal. Song and Ermon (2019) report that sliced score matching needs around four times the computation of denoising score matching.

Langevin Dynamics and Annealing

Langevin dynamics iterates

$$x_k = x_{k-1} + \frac{\eta}{2}\, \nabla_x \log p(x_{k-1}) + \sqrt{\eta}\, z_k, \qquad z_k \sim \mathcal{N}(0, I),$$

with step size $\eta$ (Song & Ermon, 2019, Eq. 4, who write $\epsilon$). As $\eta \to 0$ and the number of steps grows, $x_k$ becomes an exact sample under regularity conditions; for finite steps they neglect the Metropolis–Hastings correction.

Song and Ermon (2019, §3) identified why a single learned score fails on images. Under the manifold hypothesis, images concentrate near a low-dimensional set where the score is undefined and score matching inconsistent; sliced score matching on raw CIFAR-10 gave a fluctuating loss that converged once nearly imperceptible noise $\mathcal{N}(0, 0.0001)$ was added (Fig. 1). In low-density regions, where few training images fall, the loss barely constrains the learned score, yet Langevin chains often start there. And even an exact score misleads: inside modes with disjoint supports the score ignores the mixture weights, so Langevin samples do too, and with nearly disjoint modes mixing needs tiny steps and very many of them (§3.2.2).

Their noise-conditional score network (NCSN) $s_\theta(x, \sigma)$ is one network trained by DSM at a geometric sequence of levels $\sigma_1 > \dots > \sigma_L$,

$$\frac{1}{L} \sum_{i=1}^{L} \sigma_i^2\, \tfrac{1}{2}\, \mathbb{E} \Big\lVert s_\theta(\tilde x, \sigma_i) + \frac{\tilde x - x}{\sigma_i^2} \Big\rVert^2 ,$$

the weight $\sigma_i^2$ balancing the terms because they observed empirically that $\lVert s_\theta \rVert \propto 1/\sigma$ near the optimum. Large noise fills empty regions and merges modes; small noise leaves images almost unchanged. Annealed Langevin dynamics (Algorithm 1) runs $T$ steps at each level from $\sigma_1$ down to $\sigma_L$, with step $\eta_i = \eta\, \sigma_i^2 / \sigma_L^2$, each level starting from the previous one’s samples. For pixels in $[0, 1]$ they used $L = 10$, $\sigma_1 = 1$, $\sigma_{10} = 0.01$, $T = 100$, $\eta = 2 \times 10^{-5}$ and a uniform-noise start, and reported an Inception score of 8.87 and a Fréchet Inception Distance (FID) of 25.32 on unconditional CIFAR-10 (Table 1). A baseline with the single level $\sigma = 0.01$ and plain Langevin dynamics failed to produce reasonable images (§5; Appendix C.1 of arXiv v3). For their two-Gaussian toy example the paper disagrees with itself: the main text gives $\sigma_1 = 10$, $\sigma_{10} = 0.1$ (§4.3), Appendix B.1 gives $\sigma_1 = 20$, $\sigma_{10} = 1$.

Song et al. (2021) took the number of levels to infinity. A forward stochastic differential equation diffuses the data, and generation solves its time reversal, whose drift contains $\nabla_x \log p_t(x)$, learned by a continuous weighted DSM loss (Eq. 7); the diffusion model article gives both SDEs. Langevin dynamics survives as the corrector of their predictor–corrector samplers, which alternate a reverse-SDE step with score-based Markov chain Monte Carlo (§4.2), and the deterministic probability-flow ODE (Eq. 13), also given there, shares the SDE’s marginals, permitting exact likelihoods and adaptive solvers. The diffusion sampling article covers these solvers.

Connection to Denoising Autoencoders

Vincent’s paper was about autoencoders. For a denoising autoencoder with tied weights, a sigmoid hidden layer, a linear decoder $r(\tilde x) = W^\top \operatorname{sigmoid}(W \tilde x + b) + c$ and Gaussian corruption of variance $\sigma^2$, he defined an energy whose score is $(r(x) - x)/\sigma^2$. Substituted into $J_{\mathrm{DSM}}$, it gives the squared reconstruction error divided by $2\sigma^4$: training the autoencoder is score matching of that energy model against a Parzen (Gaussian-kernel) estimate of the data density. The denoiser’s residual $r(\tilde x) - \tilde x$ is a scaled score, pointing from a corrupted image back toward the data. Vincent noted that this gives the autoencoder a proper energy, so inputs can be ranked and, in principle, sampled. In our reading, a diffusion denoiser $\epsilon_\theta$ is the same idea without the energy, with one network for all noise levels; NCSN likewise outputs the vector field directly rather than differentiating an energy, to avoid higher-order gradients (Song & Ermon, 2019, §2.1).

Examples

The three objectives on image patches. The data are three $8 \times 8$ patterns, a vertical edge, a horizontal edge and a centered square (values $\pm 1$, mutually orthogonal), with weights 0.6, 0.3 and 0.1 and pixel noise of standard deviation 0.1, so every noisy marginal is a mixture of Gaussians with an exact score. At a step with $\bar\alpha_t = 0.5$ the code averages the DSM target over the posterior of $x_0$ at four noisy patches, then evaluates $J_{\mathrm{ESM}}$, $J_{\mathrm{ISM}}$ (with an exact Jacobian trace) and $J_{\mathrm{DSM}}$ for four candidate score models on 200,000 training pairs.

import numpy as np

# Three 8x8 patterns in [-1, 1]: vertical edge, horizontal edge, centered square.
# They have equal norm (8) and are mutually orthogonal.
yy, xx = np.mgrid[:8, :8]
mu = np.array([np.where(xx >= 4, 1.0, -1.0),                          # vertical edge
               np.where(yy >= 4, 1.0, -1.0),                          # horizontal edge
               np.where((abs(yy - 3.5) < 2) & (abs(xx - 3.5) < 2), 1.0, -1.0)]
              ).reshape(3, 64)
w = np.array([0.6, 0.3, 0.1])          # mixture weights
tau2 = 0.1 ** 2                        # within-pattern pixel noise: x0 = mu_k + 0.1 * xi
D = 64

# One diffusion step with alpha_bar = 0.5: x_t = sqrt(ab) x0 + sqrt(1 - ab) eps.
ab = 0.5
v = ab * tau2 + 1 - ab                 # p_t is a mixture of N(sqrt(ab) mu_k, v I)

def resp(xt, weights=w):
    d2 = ((xt[:, None, :] - np.sqrt(ab) * mu[None]) ** 2).sum(-1)
    lg = np.log(weights) - d2 / (2 * v)
    r = np.exp(lg - lg.max(1, keepdims=True))
    return r / r.sum(1, keepdims=True)

def true_score(xt, weights=w):         # exact grad log p_t(x_t)
    return -(xt - np.sqrt(ab) * resp(xt, weights) @ mu) / v

def divergence(xt, weights=w):         # trace of the Jacobian of true_score
    r = resp(xt, weights)
    g = -(xt[:, None, :] - np.sqrt(ab) * mu[None]) / v          # grad of each log-component
    dr = r[:, :, None] * (g - (r[:, :, None] * g).sum(1, keepdims=True))
    return -D / v + np.sqrt(ab) / v * np.einsum("nkd,kd->n", dr, mu)

rng = np.random.default_rng(0)
def sample_x0(n):
    return mu[rng.choice(3, n, p=w)] + np.sqrt(tau2) * rng.standard_normal((n, D))

# 1) The DSM target -eps / sqrt(1 - ab), averaged over the clean images that could
#    have produced a given noisy x_t, equals the true score at x_t.
xt_test = np.sqrt(ab) * sample_x0(4) + np.sqrt(1 - ab) * rng.standard_normal((4, D))
pv = 1 / (1 / tau2 + ab / (1 - ab))    # posterior variance of x0 within a component
print("x_t   |true score|   |MC mean of DSM target - true score|   (M = 100,000)")
for i, xt in enumerate(xt_test):
    r = resp(xt[None])[0]
    M = 100_000
    k = rng.choice(3, M, p=r)                                    # component | x_t
    m = pv * (mu[k] / tau2 + np.sqrt(ab) * xt / (1 - ab))        # posterior mean
    x0 = m + np.sqrt(pv) * rng.standard_normal((M, D))           # x0 ~ q(x0 | x_t)
    eps = (xt - np.sqrt(ab) * x0) / np.sqrt(1 - ab)
    target = -eps / np.sqrt(1 - ab)                              # grad log q(x_t | x0)
    s = true_score(xt[None])[0]
    print(f"{i:3d} {np.linalg.norm(s):12.3f} {np.linalg.norm(target.mean(0) - s):18.4f}")

# 2) Explicit, implicit and denoising score matching differ by constants: compare candidate
#    score models s(x_t) on the same 200,000 training pairs (x0, eps).
N = 200_000
x0 = sample_x0(N)
eps = rng.standard_normal((N, D))
xt = np.sqrt(ab) * x0 + np.sqrt(1 - ab) * eps
s_true = true_score(xt)
models = {"true score": (true_score, divergence),
          "0.5 x true score": (lambda x: 0.5 * true_score(x), lambda x: 0.5 * divergence(x)),
          "uniform weights": (lambda x: true_score(x, np.ones(3) / 3),
                              lambda x: divergence(x, np.ones(3) / 3)),
          "N(0, I) score": (lambda x: -x, lambda x: -D * np.ones(len(x)))}
print(f"\n{'model':18s}{'J_ESM':>10s}{'J_ISM':>10s}{'J_DSM':>10s}{'ISM - ESM':>11s}{'DSM - ESM':>11s}")
for name, (f, div) in models.items():
    s = f(xt)
    j_esm = 0.5 * ((s - s_true) ** 2).sum(1).mean()                 # needs the true score
    j_ism = (0.5 * (s ** 2).sum(1) + div(xt)).mean()                # Hyvarinen (2005)
    j_dsm = 0.5 * ((s + eps / np.sqrt(1 - ab)) ** 2).sum(1).mean()  # Vincent (2011)
    print(f"{name:18s}{j_esm:10.3g}{j_ism:10.3f}{j_dsm:10.3f}{j_ism - j_esm:11.3f}{j_dsm - j_esm:11.3f}")

Output:

x_t   |true score|   |MC mean of DSM target - true score|   (M = 100,000)
  0       11.402             0.0032
  1       10.411             0.0031
  2       11.559             0.0036
  3       11.792             0.0033

model                  J_ESM     J_ISM     J_DSM  ISM - ESM  DSM - ESM
true score                 0   -63.395     0.634    -63.395      0.634
0.5 x true score        15.8   -47.532    16.466    -63.366      0.632
uniform weights     4.46e-10   -63.395     0.634    -63.395      0.634
N(0, I) score           31.5   -31.846    32.147    -63.365      0.629

The averaged target matches the true score, of norm 10.4 to 11.8, to within 0.0031 to 0.0036, about the Monte Carlo error of 100,000 draws by our calculation. $J_{\mathrm{ISM}} - J_{\mathrm{ESM}}$ stays between $-63.395$ and $-63.365$ and $J_{\mathrm{DSM}} - J_{\mathrm{ESM}}$ between 0.629 and 0.634 across candidates, differences of Monte Carlo size, so the three objectives rank models identically while only $J_{\mathrm{ESM}}$ needs the true score. The score built with uniform weights does as well as the true one: at this noise level the modes are so far apart that the weights barely change the score where the data lie. It runs in about 2 seconds.

One noise level versus annealing. With Song and Ermon’s perturbation $\tilde x = x + \sigma z$ (a rescaling of the diffusion step; $\bar\alpha_t = 0.5$ is $\sigma = 1$), the code measures how much uniform weights change the score at ten levels from $\sigma = 20$ to $0.1$, then samples 1,000 patches per seed from uniform noise with the exact score: plain Langevin dynamics at $\sigma = 0.1$, with the true and with reversed weights, and annealed Langevin dynamics with 100 and 300 steps per level, and assigns each sample to its nearest pattern.

import numpy as np

# Same data as before: three orthogonal 8x8 patterns plus pixel noise of std 0.1.
yy, xx = np.mgrid[:8, :8]
mu = np.array([np.where(xx >= 4, 1.0, -1.0),                          # vertical edge
               np.where(yy >= 4, 1.0, -1.0),                          # horizontal edge
               np.where((abs(yy - 3.5) < 2) & (abs(xx - 3.5) < 2), 1.0, -1.0)]
              ).reshape(3, 64)
w = np.array([0.6, 0.3, 0.1])
tau2 = 0.1 ** 2

# Song & Ermon perturbation: x~ = x + sigma * z, so q_sigma is a mixture of
# N(mu_k, (tau^2 + sigma^2) I) and its score is exact.
def score(x, sigma, weights=w):
    v = tau2 + sigma ** 2
    d2 = (x ** 2).sum(1, keepdims=True) - 2 * x @ mu.T + (mu ** 2).sum(1)
    lg = np.log(weights) - d2 / (2 * v)
    r = np.exp(lg - lg.max(1, keepdims=True))
    r /= r.sum(1, keepdims=True)
    return -(x - r @ mu) / v

sigmas = np.geomspace(20.0, 0.1, 10)    # sigma_1 = 20 ... sigma_L = 0.1, L = 10
eps_step = 0.1 * sigmas[-1] ** 2       # so that alpha_i = 0.1 sigma_i^2

# How much do the weights matter to the score at each noise level? Fisher divergence
# between the score with the true weights and with uniform weights, under q_sigma.
rng = np.random.default_rng(0)
x0 = mu[rng.choice(3, 20_000, p=w)] + 0.1 * rng.standard_normal((20_000, 64))
print("sigma   0.5 E||s(true w) - s(uniform w)||^2")
for s in sigmas:
    x = x0 + s * rng.standard_normal(x0.shape)
    print(f"{s:6.2f}  {0.5 * ((score(x, s) - score(x, s, np.ones(3) / 3)) ** 2).sum(1).mean():.2e}")

def langevin(x, levels, steps, rng, weights=w):
    # Song & Ermon (2019), Algorithm 1: step alpha_i = eps * sigma_i^2 / sigma_L^2.
    for s in levels:
        a = eps_step * s ** 2 / sigmas[-1] ** 2
        for _ in range(steps):
            x = x + 0.5 * a * score(x, s, weights) + np.sqrt(a) * rng.standard_normal(x.shape)
    return x

def shares(x):
    d2 = ((x[:, None, :] - mu[None]) ** 2).mean(-1)
    k = d2.argmin(1)
    return np.bincount(k, minlength=3) / len(k), np.sqrt(d2.min(1).mean())

n = 1000
runs = [("single sigma = 0.1", sigmas[-1:], 3000, w),
        ("single, weights reversed", sigmas[-1:], 3000, w[::-1]),
        ("annealed, T = 100", sigmas, 100, w),
        ("annealed, T = 300", sigmas, 300, w)]
print(f"\nshare of samples nearest each pattern: mean over seeds 0-4 [min, max], {n} samples per seed")
print(f"{'sampler':25s}{'evals':>6s}{'v-edge':>21s}{'h-edge':>21s}{'square':>21s}{'RMS':>7s}")
for name, levels, steps, weights in runs:
    res = []
    for seed in range(5):
        rng = np.random.default_rng(seed)
        x = rng.uniform(-1, 1, (n, 64))            # uniform noise, as in Song & Ermon
        f, rms = shares(langevin(x, levels, steps, rng, weights))
        res.append(np.append(f, rms))
    m, lo, hi = np.mean(res, 0), np.min(res, 0), np.max(res, 0)
    print(f"{name:25s}{len(levels) * steps:6d}"
          + "".join(f"{a:7.3f} [{b:.3f}, {c:.3f}]" for a, b, c in zip(m[:3], lo[:3], hi[:3]))
          + f"{m[3]:7.3f}")
print(f"{'true weights':31s}" + "".join(f"{v:7.3f}" + " " * 14 for v in w))

Output:

sigma   0.5 E||s(true w) - s(uniform w)||^2
 20.00  2.41e-05
 11.10  2.26e-04
  6.16  1.59e-03
  3.42  5.16e-03
  1.90  1.71e-03
  1.05  3.29e-14
  0.58  0.00e+00
  0.32  0.00e+00
  0.18  0.00e+00
  0.10  0.00e+00

share of samples nearest each pattern: mean over seeds 0-4 [min, max], 1000 samples per seed
sampler                   evals               v-edge               h-edge               square    RMS
single sigma = 0.1         3000  0.329 [0.308, 0.343]  0.331 [0.302, 0.352]  0.340 [0.328, 0.355]  0.142
single, weights reversed   3000  0.327 [0.307, 0.338]  0.330 [0.304, 0.351]  0.343 [0.332, 0.360]  0.142
annealed, T = 100          1000  0.578 [0.569, 0.589]  0.300 [0.294, 0.307]  0.123 [0.117, 0.129]  0.143
annealed, T = 300          3000  0.583 [0.562, 0.608]  0.301 [0.272, 0.338]  0.116 [0.098, 0.127]  0.142
true weights                     0.600                0.300                0.100              

The weights matter most between $\sigma \approx 2$ and $6$, peaking at $\sigma = 3.42$, less at the two largest levels, where the modes merge, and not at all, to rounding, at $\sigma \le 1.05$. Langevin dynamics at $\sigma = 0.1$ therefore ignores them, producing about one third of each pattern whether the vertical edge has weight 0.6 or 0.1; in our reading each chain settles into the pattern nearest its start, and our three patterns are symmetric. Annealing recovers the weights, averaging 0.578, 0.300 and 0.123 with 100 steps per level and 0.583, 0.301 and 0.116 with 300 over seeds 0 to 4. The excess of the rarest pattern shrinks with more steps, suggesting incomplete mixing at the levels where the weights are decided. All samplers end at an RMS distance of 0.142 to 0.143 from the nearest pattern, against $\sqrt{0.1^2 + 0.1^2} \approx 0.141$ for data noise plus the final level. It runs in under 30 seconds.

Common Misconceptions

• “Denoising score matching learns the data score.” It learns the score of the noisy data, which approximates the data score only for small noise (Song & Ermon, 2019), and a single small level leaves low-density regions unlearned.
• “An exact score gives correct samples.” Plain Langevin dynamics with the exact score misallocates well-separated modes; annealing is part of the method.

Where It Is Used

• Generative models. The losses of diffusion models, NCSN and the score SDE are weighted denoising score matching.
• Restoration and editing. Song et al. (2021) performed class-conditional generation, inpainting and colorization with one unconditional score model, and Song and Ermon (2019) inpainted with annealed Langevin dynamics. A learned score is an image prior, and a score model is a denoiser.
• Unnormalized models. Hyvärinen’s original setting, including natural-image patch models.
• Implicit distributions. Song et al. (2019) used sliced score matching to train VAEs with implicit encoders and Wasserstein autoencoders.