# Tweedie Moment Projected Diffusions for
Inverse Problems

Benjamin Boys & Mark Girolami
University of Cambridge & The Alan Turing Institute
{bb515<EMAIL_ADDRESS>Jakiw Pidstrigach & Sebastian Reich
Universitat Potsdam
{pidstrigach<EMAIL_ADDRESS>Alan Mosca
nplan.io
<EMAIL_ADDRESS>O. Deniz Akyildiz
Imperial College London
<EMAIL_ADDRESS>
###### Abstract

Diffusion generative models unlock new possibilities for inverse problems as
they allow for the incorporation of strong empirical priors into the process
of scientific inference. Recently, diffusion models received significant
attention for solving inverse problems by posterior sampling, but many
challenges remain open due to the intractability of this sampling process.
Prior work resorted to Gaussian approximations to conditional densities of the
reverse process, leveraging Tweedie’s formula to parameterise its mean,
complemented with various heuristics. In this work, we leverage higher order
information using Tweedie’s formula and obtain a finer approximation with a
principled covariance estimate. This novel approximation removes any time-
dependent step-size hyperparameters required by earlier methods, and enables
higher quality approximations of the posterior density which results in better
samples. Specifically, we tackle noisy linear inverse problems and obtain a
novel approximation to the gradient of the likelihood. We then plug this
gradient estimate into various diffusion models and show that this method is
optimal for a Gaussian data distribution. We illustrate the empirical
effectiveness of our approach for general linear inverse problems on toy
synthetic examples as well as image restoration using pretrained diffusion
models as the prior. We show that our method improves the sample quality by
providing statistically principled approximations to diffusion posterior
sampling problem.

##  1 Introduction

Solving inverse problems is one of the most crucial goals in scientific
inference (Tarantola, 2005; Stuart, 2010). Diffusion models (Song et al.,
2020; Ho et al., 2020) recently emerged as a strong alternative to classical
Bayesian approaches to solve inverse problems to build samplers for
complicated and implicit conditional distributions. With their flexibility,
these models replace handcrafted priors on the latent signal with pretrained
and strong empirical priors. For example, given a latent signal
${\mathbf{x}}_{0}$ (say a face image), we can train a diffusion model to
sample from the prior $p({\mathbf{x}}_{0})$. The idea for solving inverse
problems is to leverage extraordinary modelling power of diffusion models to
learn samplers for priors and couple this with a given likelihood
$p({\mathbf{y}}|{\mathbf{x}}_{0})$ for a given data ${\mathbf{y}}$ to sample
from posterior densities for inverse problems. However, designing a diffusion
model for the posterior comes with its own challenges due to intractability.
Despite its challenges, this approach has been recently taking off with lots
of activity in the field, e.g. for compressed sensing (Bora et al., 2017;
Kadkhodaie & Simoncelli, 2021), projecting score-based stochastic differential
equations (SDEs) (Song et al., 2021b), gradient-based approaches (Daras et
al., 2022; Chung et al., 2022b), magnetic reasonance imaging (MRI) by
approximating annealed Langevin dynamics with approximate scores (Jalal et
al., 2021), image restoration (Kawar et al., 2022), score-based models as
priors but with a normalising flow approach (Feng et al., 2023), variational
approaches (Mardani et al., 2023; Feng & Bouman, 2023). Most relevant ideas to
us, which we will review in detail in Section 5, use Tweedie’s formula (Efron,
2011; Laumont et al., 2022) to approximate the smoothed likelihood, e.g.
Diffusion posterior sampling (DPS) (Chung et al., 2022a; 2023) and pseudo-
guided diffusion models ($\Pi$GDM) (Song et al., 2023). Similar approaches are
also exploited using singular-value decomposition (SVD) based approaches
(Kawar et al., 2021). In this work, we develop an approach which builds a
tighter approximation to optimal formulae for approximating the scores of the
posterior diffusion model.

Figure 1: (Plot) error to target posterior for a Gaussian random field. (Top
row) visualisation of the empirical mean and variance of the 1500 samples that
were used to compute this error against the analytical moments. (Bottom)
Wasserstein distances of different methods w.r.t. sample size. For details,
see Appendix E.2.

This paper is devoted developing novel methods to solving inverse problems,
given a latent (target) signal ${\mathbf{x}}_{0}\in\mathbb{R}^{d_{x}}$, noisy
observed data ${\mathbf{y}}\in\mathbb{R}^{d_{y}}$, a known linear observation
map ${\mathbf{H}}$, and a pretrained diffusion prior. The main tool we use is
Tweedie’s formula to obtain both the mean and the covariance for approximating
diffused likelihoods, to be used for building the final posterior score
approximation. This is as opposed to previous works which only utilised first
moment approximations using Tweedie’s formula (Chung et al., 2022a; Song et
al., 2023). We show that utilising covariance approximation, with moment
projections, provides a principled scheme with improved performance, which we
term Tweedie Moment Projected Diffusions (TMPD).

To demonstrate our method briefly, Figure 1 demonstrates a sampling scenario
of a Gaussian random field (GRF) whose mean and variance entries are plotted
under “Analytic” column111We only plot the variances for visualisation while
the GRF has a full covariance.. We demonstrate the approximations under this
setting provided by our method (TMPD) and its diagonal (cheaper) version
(DTMPD), compared with $\Pi$GDM (Song et al., 2023), and DPS (Chung et al.,
2022a). The figure demonstrates the optimality our method: Our first and
second moment approximations become exact in this case. This results in a
drastic performance improvement stemming from the statistical optimality of
our method for near-Gaussian settings and also unlocks a possible line for
theoretical research for understanding similar diffusion models for inverse
problems.

In what follows, we will first introduce the technical background in Section 2
and then describe TMPD in detail in Section 3. We will then provide some
theoretical results about our method in Section 4 and provide a discussion to
closely related work in literature in Section 5. Finally, Section 6 will
present experiments on Gaussian mixtures, image inpainting, and super-
resolution, demonstrating quantitative and qualitative improvements provided
by TMPD.

##  2 Technical Background

There are two main diffusion modelling paradigms, namely, score-based models
(SGMs) (Song et al., 2020) and denoising diffusion probabilistic models (DDPM)
(Ho et al., 2020). In both, the goal is to sample from a target distribution
$p_{0}:=p_{\text{data}}$. To that end, an artificial path $p_{t}$ is
introduced, with the property that $p_{t}$ will approach $\mathcal{N}(0,I)$
for large $t$, i.e., $p_{t}\to\mathcal{N}(0,I_{d})$ as $t\to\infty$. Then, one
learns to reverse this process, so that one can sample from a standard normal
distribution and transform these samples to samples from $p_{\text{data}}$. In
the SDE paradigm, the interpolation parameter $t$ will take continuous values
in $t\in[0,T]$, whereas in the DDPM setting, $t$ is discrete. However, the
DDPM setting can be seen as discretization of the SDE (Song et al., 2020).

In SGMs, a stochastic differential equation (SDE) is used to noise the data. A
very common choice is to use the time-rescaled Ornstein-Uhlenbeck process:

$\mathrm{d}{\mathbf{x}}_{t}=-\frac{1}{2}\beta(t){\mathbf{x}}_{t}\mathrm{d}t+\sqrt{\beta(t)}\mathrm{d}{\mathbf{w}}_{t},\quad{\mathbf{x}}_{0}\sim
p_{0}=p_{\text{data}}.$ (1)

The corresponding reverse SDE is then given by

$\mathrm{d}{\mathbf{z}}_{t}=\frac{1}{2}\beta(T-t){\mathbf{z}}_{t}\mathrm{d}t+\beta(T-t)\nabla_{{\mathbf{z}}_{t}}\log
p_{T-t}({\mathbf{z}}_{t})\mathrm{d}t+\sqrt{\beta(T-t)}\mathrm{d}\bar{{\mathbf{w}}}_{t},\quad{\mathbf{z}}_{0}\sim
p_{T}.$ (2)

The transition kernels for the time-changed SDE can be derived from the
transition kernels of a normal OU process
$p_{t}({\mathbf{x}}_{t}|{\mathbf{x}}_{0})=\mathcal{N}({\mathbf{x}}_{t};\sqrt{\alpha_{t}}{\mathbf{x}}_{0},v_{t}{\mathbf{I}}_{d_{{\textnormal{x}}}})$
where $\alpha_{t}:=\exp\left(-\int_{0}^{t}\beta(s)\mathrm{d}s\right)$ and
$v_{t}:=1-\alpha_{t}$. A parameterisation that performs well in practice is
$\beta(t)=\beta_{\text{min}}+t(\beta_{\text{max}}-\beta_{\text{min}})$. In the
diffusion modelling literature, the OU process is also sometimes called
variance preserving SDE. However, this is not the only SDE that is suitable
for the forward process. See Appendix C for details on a time rescaled
Brownian motion (Variance Exploding (Song et al., 2020) SDE).

There are two usual approximations to solve equation 2. First, we do not know
$p_{T}$, since it is a noised version of the distribution $p_{\text{data}}$.
However for $T$ large enough, we can approximate $p_{T}\approx
p_{\text{ref}}=\mathcal{N}(0,I_{d})$. We also do not have $\nabla\log p_{T-t}$
which we need for the drift in equation 2. This can be circumvented by
approximating the drift using score-matching techniques (Hyvärinen, 2005; Ho
et al., 2020). These methods construct an estimate of the score function by
solving the score matching problem in the form of
$\mathsf{s}_{\theta}({\mathbf{x}}_{t},t)\approx\nabla_{{\mathbf{x}}_{t}}\log
p_{t}({\mathbf{x}}_{t})$. This score can also be used in the setting of DDPM
(Ho et al., 2020).

###  2.1 Conditional sampling for the linear inverse problem

In the preceding section we introduced diffusion models as a method to sample
from a target distribution $p_{\text{data}}$. We now suppose that we have
access to measurements, or observations
${\mathbf{y}}\in\mathbb{R}^{d_{{\textnormal{y}}}}$ of
${\mathbf{x}}_{0}\in\mathbb{R}^{d_{{\textnormal{x}}}}$:

${\mathbf{y}}={\mathbf{H}}{\mathbf{x}}_{0}+{\mathbf{u}},\quad{\mathbf{u}}\sim\mathcal{N}(0,\sigma_{{\textnormal{y}}}^{2}{\mathbf{I}}_{d_{{\textnormal{y}}}}).$
(3)

We would then be interested in sampling from the conditional distribution of
${\mathbf{x}}_{0}$ given ${\mathbf{y}}$, i.e.,
$p_{\text{data}}(~{}\cdot~{}|{\mathbf{y}})$. To that end, we have to modify
the reverse SDE. Optimally, we would like to sample from the reverse SDE to
the forward SDE started in $p_{\text{data}}(~{}\cdot~{}|{\mathbf{y}})$,
instead of the one started in $p_{\text{data}}$.

Optimally, we would want to replace the score $\nabla_{{\mathbf{z}}_{t}}\log
p_{T-t}({\mathbf{z}}_{t})$ in equation 2 with the posterior score
$\nabla_{{\mathbf{z}}_{t}}\log
p_{T-t|\mathbf{y}}({\mathbf{z}}_{t}|{\mathbf{y}})$. Written in terms of the
forward process, this coincides with

$\nabla_{{\mathbf{x}}_{t}}\log
p_{t|{\mathbf{y}}}({\mathbf{x}}_{t}|{\mathbf{y}})=\nabla_{{\mathbf{x}}_{t}}\log
p_{t}({\mathbf{x}}_{t})+\nabla_{{\mathbf{x}}_{t}}\log
p_{{\mathbf{y}}|t}({\mathbf{y}}|{\mathbf{x}}_{t}).$ (4)

The term $p_{{\mathbf{y}}|t}({\mathbf{y}}|{\mathbf{x}}_{t})$ is given by the
integral

$p_{{\mathbf{y}}|t}({\mathbf{y}}|{\mathbf{x}}_{t})=\int
p_{{\mathbf{y}}|0}({\mathbf{y}}|{\mathbf{x}}_{0})p_{0|t}({\mathbf{x}}_{0}|{\mathbf{x}}_{t})\mathrm{d}{\mathbf{x}}_{0},$
(5)

which involves a marginalization over ${\mathbf{x}}_{0}$. The term
$p_{0|t}({\mathbf{x}}_{0}|{\mathbf{x}}_{t})$ however is only defined
implicitly through running the diffusion model, making the above integral
difficult to evaluate. One way around this is to train a neural network to
directly approximate $\nabla\log
p_{t|{\mathbf{y}}}({\mathbf{x}}_{t}|{\mathbf{y}})$ (Batzolis et al., 2021).
Alternatively, if one already has access to an approximation of
$\nabla_{{\mathbf{x}}_{t}}\log p_{t}({\mathbf{x}}_{t})$, one can train an
auxiliary network to approximate the term $\nabla\log
p_{{\mathbf{y}}|t}({\mathbf{y}}|{\mathbf{x}}_{t})$ in equation 4, (Song et
al., 2020; Dhariwal & Nichol, 2021). However, these methods can be time and
training-data intensive, as one needs to retrain networks for each conditional
task as well have access to paired training data from the joint distribution
$({\mathbf{x}}_{0},{\mathbf{y}})$. Alternatively, one could indeed try to do a
Monte-Carlo approximation of the score corresponding to equation 5. However,
to run this estimator, one would need to evaluate the probability flow ODE
together with its derivative (Song et al., 2020, Section D.2) for each MC-
sample, which is computationally expensive. Furthermore, gradient estimators
of this kind are known to have high variance (Mohamed et al., 2020, Section
3).

##  3 Tweedie Moment Projected Diffusions

In this section, we first introduce Tweedie moment projections in Section 3.1
below. The central to our method is the approximation
$p_{0|t}({\mathbf{x}}_{0}|{\mathbf{x}}_{t})\approx\mathcal{N}\left({\mathbf{x}}_{0};{\mathbf{m}}_{0|t}({\mathbf{x}}_{t}),{\mathbf{C}}_{0|t}({\mathbf{x}}_{t})\right)$.
In that case, since the conditional distribution of ${\mathbf{y}}$ given
${\mathbf{x}}_{0}$ is also Gaussian, we can _compute the integral in equation
5 analytically_ — $p_{{\mathbf{y}}|t}({\mathbf{y}}|{\mathbf{x}}_{t})$ will
just be another Gaussian in that case, its mean and covariance being
determined through ${\mathbf{m}}_{0|t}$, ${\mathbf{C}}_{0|t}$, ${\mathbf{H}}$
and $\sigma_{y}$. In particular, we can then use this Gaussian to approximate
$\nabla\log p_{{\mathbf{y}}|t}({\mathbf{y}}|{\mathbf{x}}_{t})$, since the
score of a Gaussian is available in closed form.

###  3.1 Tweedie moment projections

Instead of just approximating the variance of
$p_{0|t}({\mathbf{x}}_{0}|{\mathbf{x}}_{t})$ heuristically, we approximate it
by projecting on to the closest Gaussian distribution using Tweedie’s formula
for the second moment. Our approximation at this stage consists of two main
steps: (1) Find the mean and covariance of
$p_{0|t}({\mathbf{x}}_{0}|{\mathbf{x}}_{t})$ using Tweedie’s formula, and (2)
approximate this density with a Gaussian using the mean and covariance of
$p_{0|t}({\mathbf{x}}_{0}|{\mathbf{x}}_{t})$ (moment projection). Due to this
approximation, we will refer to resulting methods as Tweedie Moment Projected
Diffusions (TMPD). We will first introduce Tweedie’s formula for the mean and
covariance and then describe the moment projection.

######  Proposition 1 (Tweedie’s formula).

Let ${{\mathbf{m}}}_{0|t}$ and ${\mathbf{C}}_{0|t}$ be the mean and the
covariance of $p_{0|t}({\mathbf{x}}_{0}|{\mathbf{x}}_{t})$, respectively. Then
given the marginal density $p_{t}({\mathbf{x}}_{t})$, the mean is given as

$\displaystyle{\mathbf{m}}_{0|t}=\mathbb{E}[{\mathbf{x}}_{0}|{\mathbf{x}}_{t}]=\frac{1}{\sqrt{\alpha_{t}}}(\mathbf{x}_{t}+v_{t}\nabla_{{\mathbf{x}}_{t}}\log
p_{t}({\mathbf{x}}_{t})).$

Then the covariance ${\mathbf{C}}_{0|t}$ is given by

$\displaystyle{\mathbf{C}}_{0|t}=\mathbb{E}\left[({\mathbf{x}}_{0}-{\mathbf{m}}_{0|t})({\mathbf{x}}_{0}-{\mathbf{m}}_{0|t})^{\top}\left.\right|{\mathbf{x}}_{t}\right]=\frac{v_{t}}{\alpha_{t}}({\mathbf{I}}_{d_{x}}+v_{t}\nabla^{2}\log
p_{t}({\mathbf{x}}_{t}))=\frac{v_{t}}{\sqrt{\alpha_{t}}}\nabla_{{\mathbf{x}}_{t}}{\mathbf{m}}_{0|t}.$
(6)

The proof is an adaptation of Meng et al. (2021, Theorem 1), see Appendix A.1.
While ${\mathbf{m}}_{0|t}$ and ${\mathbf{C}}_{0|t}$ gives us the moments of
the density $p_{0|t}({\mathbf{x}}_{0}|{\mathbf{x}}_{t})$, we do not have the
exact form of this density. At this stage, we employ moment projection, i.e.,
we choose the closest Gaussian in Kullback-Leibler (KL) divergence to a
distribution with same moments as formalised next.

######  Proposition 2 (Moment projection).

Let $p_{0|t}({\mathbf{x}}_{0}|{\mathbf{x}}_{t})$ be a distribution with mean
${\mathbf{m}}_{0|t}$ and covariance ${\mathbf{C}}_{0|t}$. Let
$\hat{p}_{0|t}({\mathbf{x}}_{0}|{\mathbf{x}}_{t})$ be the the closest Gaussian
in KL divergence to $p_{0|t}({\mathbf{x}}_{0}|{\mathbf{x}}_{t})$, i.e.,

$\displaystyle\hat{p}_{0|t}({\mathbf{x}}_{0}|{\mathbf{x}}_{t})=\arg\min_{q\in\mathcal{Q}}D_{\mathrm{KL}}(p_{0|t}({\mathbf{x}}_{0}|{\mathbf{x}}_{t})||q),$

where $\mathcal{Q}$ is the family of multivariate Gaussian distributions. Then

$\displaystyle\hat{p}_{0|t}({\mathbf{x}}_{0}|{\mathbf{x}}_{t})=\mathcal{N}({\mathbf{x}}_{0};{\mathbf{m}}_{0|t},{\mathbf{C}}_{0|t}).$

This is a well-known moment matching result, see, e.g., Bishop (2006, Section
10.7). Merging Propositions 1 and 2 leads to the following Tweedie moment
projection:

$\displaystyle p_{0|t}({\mathbf{x}}_{0}|{\mathbf{x}}_{t})$
$\displaystyle\approx\mathcal{N}\left({\mathbf{x}}_{0};{\mathbf{m}}_{0|t},\frac{v_{t}}{\sqrt{\alpha_{t}}}\nabla_{{\mathbf{x}}_{t}}{\mathbf{m}}_{0|t}\right),$
(7)

where ${\mathbf{m}}_{0|t}$ is given in Proposition 1. In the next section, we
demonstrate how to use this approximation to obtain approximate likelihoods.

###  3.2 Tweedie Moment Projected Likelihood Approximation

We next use the approximation in equation 7 to compute the following integral
analytically

$\displaystyle p_{{\mathbf{y}}|t}({\mathbf{y}}|{\mathbf{x}}_{t})$
$\displaystyle=\int
p_{{\mathbf{y}}|0}({\mathbf{y}}|{\mathbf{x}}_{0})p_{0|t}({\mathbf{x}}_{0}|{\mathbf{x}}_{t})\mathrm{d}{\mathbf{x}}_{t}=\mathcal{N}\left({\mathbf{y}};{\mathbf{H}}{\mathbf{m}}_{0|t},{\mathbf{H}}\frac{v_{t}}{\sqrt{\alpha_{t}}}\nabla_{{\mathbf{x}}_{t}}{\mathbf{m}}_{0|t}{\mathbf{H}}^{\top}+\sigma^{2}_{{\textnormal{y}}}{\mathbf{I}}_{d_{{\textnormal{y}}}}\right)$
(8)

which leads to the approximation

$f^{\mathbf{y}}({\mathbf{x}}_{t}):=\nabla_{{\mathbf{x}}_{t}}{\mathbf{m}}_{0|t}{\mathbf{H}}^{\top}({\mathbf{H}}\frac{v_{t}}{\sqrt{\alpha_{t}}}\nabla_{{\mathbf{x}}_{t}}{\mathbf{m}}_{0|t}{\mathbf{H}}^{\top}+\sigma^{2}_{{\textnormal{y}}}{\mathbf{I}}_{d_{{\textnormal{y}}}})^{-1}({\mathbf{y}}-{\mathbf{H}}{\mathbf{m}}_{0|t})\approx\nabla_{{\mathbf{x}}_{t}}\log
p_{{\mathbf{y}}|t}({\mathbf{y}}|{\mathbf{x}}_{t}).$ (9)

We note that this approximation is exact for the case $p_{\text{data}}$ is a
Gaussian. However, otherwise, there are several approximations behind equation
9. In particular, the approximation in equation 9 treats ${\mathbf{C}}_{0|n}$
constant w.r.t. ${\mathbf{x}}_{t}$ when computing the gradient (which is the
case if $p_{\rm{data}}$ is Gaussian). For non-Gaussian $p_{\rm{data}}$, this
results in computationally efficient sampler, as otherwise the resulting terms
can be expensive to compute.

###  3.3 Algorithms

Plugging the approximation in equation 9 into the reverse SDE in equation 2
together with the prior score as described in Section 2.1 results in a TMPD
for conditional sampling to solve inverse problems. The SDE we will
approximate numerically to sample from the conditional distribution is given
by

$\displaystyle\mathrm{d}{\mathbf{z}}_{t}=\frac{1}{2}\beta(T-t){\mathbf{z}}_{t}\mathrm{d}t+\beta(T-t)(\nabla_{{\mathbf{z}}_{t}}\log
p_{T-t}({\mathbf{z}}_{t})+f^{{\mathbf{y}}}_{T-t}({\mathbf{z}}_{t}))\mathrm{d}t+\sqrt{\beta(T-t)}\mathrm{d}\bar{{\mathbf{w}}}_{t}$
(10)

where ${\mathbf{z}}_{0}\sim p_{T}$ and
$f^{{\mathbf{y}}}_{t}({\mathbf{z}}_{t})\approx\nabla\log
p_{{\mathbf{y}}|T-t}({\mathbf{y}}|{\mathbf{z}}_{t})$ is our approximation to
the data likelihood, given by equation 9. We call this SDE the _TMPD SDE_.

We have two options to convert TMPD SDE into implementable methods: (1) score-
based samplers which we abbreviate also as TMPD as they are Euler-Maruyama
discretizations of the TMPD SDE akin to Song et al. (2020), (2) denoising
diffusion models (TMPD-D). The denoising diffusion approach is derived from
approximate reverse Markov Chains and is the approach of DDPM and DPS methods
(Ho et al., 2020; Chung et al., 2022a). We note that the Gaussian projection
can be used in this discrete setting, assuming that the conditional density is
available analytically as in Ho et al. (2020), and can be written as
$p_{n|0}({\mathbf{x}}_{n}|{\mathbf{x}}_{0})=\mathcal{N}({\mathbf{x}}_{n};\sqrt{\alpha_{n}}{\mathbf{x}}_{0},v_{n}{\mathbf{I}}_{d_{x}})$.
The idea is to update the unconditional mean
${\mathbf{m}}_{0|n}({\mathbf{x}}_{n})$ of the density
$p_{0|n}({\mathbf{x}}_{0}|{\mathbf{x}}_{n})$ with a Bayesian update:
$p({\mathbf{x}}_{0}|{\mathbf{x}}_{n},{\mathbf{y}})\propto
p({\mathbf{y}}|{\mathbf{x}}_{n})p_{0|n}({\mathbf{x}}_{0}|{\mathbf{x}}_{n})$.
Given a similar formulation as above, assuming we have a readily available
approximation
$p_{0|n}({\mathbf{x}}_{0}|{\mathbf{x}}_{n})\approx\mathcal{N}({\mathbf{x}}_{0};{\mathbf{m}}_{0|n},{\mathbf{C}}_{0|n})$
and a likelihood similar to equation 8 where $t$ can be replaced by $n$. Then
we can compute the moments of
$p({\mathbf{x}}_{0}|{\mathbf{x}}_{n},{\mathbf{y}})$ analytically, which we
denote ${\mathbf{m}}_{0|n}^{{\mathbf{y}}}$ and
${\mathbf{C}}_{0|n}^{\mathbf{y}}$. The Bayes update for Gaussians gives
(Bishop, 2006)

$\displaystyle{\mathbf{m}}_{0|n}^{{\mathbf{y}}}$
$\displaystyle={\mathbf{m}}_{0|n}+{\mathbf{C}}_{0|n}{\mathbf{H}}^{\top}({\mathbf{H}}{\mathbf{C}}_{0|n}{\mathbf{H}}^{\top}+\sigma_{y}^{2}{\mathbf{I}}_{d_{x}})^{-1}({\mathbf{y}}-{\mathbf{C}}_{0|n}{\mathbf{m}}_{0|n}),$
(11)

Incorporating equation 11 into the usual Ancestral sampling (Ho et al., 2020)
steps leads to Algorithm 1, termed TMPD-Denoising (TPMD-D). The update in
equation 11 can be used in any discrete sampler such as denoising diffusion
implicit models (DDIM) (Song et al., 2021a).

We show in our experiments promising results for TMPD that motivate the
exploration of less computationally expensive approximations to the full
Jacobian. We also empirically study a computationally cheaper method that is
applicable to inpainting and that does not require computing the full Jacobian
of the score network, see Appendix E.1.

##  4 Theoretical Guarantees

Because of the approximations, our method, as well as $\Pi$GDM Song et al.
(2023) and DPS (Chung et al., 2022a) do not sample the exact posterior for
general prior distributions. Therefore, one cannot hope for these methods to
sample the true posterior and a priori it is not even clear how the sampled
distribution relates to the true posterior. Without further justification,
such methods should only be interpreted as _guidance methods_ , where paths
are guided to regions in which make a given observation ${\mathbf{y}}$ more
likely, not as posterior sampling methods.

We justify our approximation by showing that the TMPD-SDE in equation 10 is
able to sample the _exact posterior_ in the Gaussian case. One can see that
this is in contrast to $\Pi$GDM and DPS in our numerical experiments or by
explicitly evaluating their approximations on simple one-dimensional examples.

######  Proposition 3.

Assume that $p_{\rm{data}}$ is Gaussian. Then, the posterior score expression
using equation 9 is exact, i.e., if there are no errors in the initial
condition and drift approximation
$s_{\theta}({\mathbf{x}}_{t},t)=\nabla_{{\mathbf{x}}_{t}}\log
p_{t}({\mathbf{x}}_{t})$, the TMPD will sample
$p_{\rm{data}}(\cdot|{\mathbf{y}})$ at its final time.

The proof is given in Appendix B.1. However, most distributions will not be
Gaussian. The following Theorem generalizes the above proposition to _non-
Gaussian distributions_ , as long as they have a density with respect to a
Gaussian. We study how close our sampled measure will be to the true posterior
distribution and give explicit bounds on the total variation distance in terms
of the regularity properties of the density:

######  Theorem 1.

Assume that the data distribution $p_{\textup{data}}$ can be written as

$p_{\rm{data}}({\mathbf{x}}_{0})=\exp(\Phi({\mathbf{x}}_{0}))\mathcal{N}({\mathbf{x}}_{0};\bm{\mu}_{0},\bm{\Sigma}_{0}),$
(12)

for some $\bm{\mu}_{0}$ and $\bm{\Sigma}_{0}$. We furthermore assume that for
some $M\geq 1$

$\displaystyle 1/M\leq\exp(\Phi({\mathbf{x}}))\leq
M\quad\quad\text{and}\quad\quad\|\nabla_{{\mathbf{x}}}\Phi({\mathbf{x}})\|\leq
L.$

for all ${\mathbf{x}}\in\mathbb{R}^{d_{x}}$. Then

$\|p_{\rm{data}}(\cdot|{\mathbf{y}})-q_{T}(\cdot)\|_{\textup{TV}}\leq
C(1+T^{1/2})\left((M^{5/2}-1)(L^{1/2}+1)+L^{1/2}\right),$

where $q_{t}$ denotes the law of the corresponding reverse-time process for
equation 10 at time $t$ and the constant $C$ that only depends on
${\mathbf{y}},{\mathbf{H}},\sigma_{y},\bm{\mu}_{0}$ and $\bm{\Sigma}_{0}$.

See Appendix B.2 for a proof.

######  Remark 1.

We note that this analysis can be extended to account for (i) discretization
error, (ii) the error in the initial distribution of the reverse SDE, (iii)
the score approximation error using the techniques from Chen et al. (2023). We
think this is a fruitful direction for our work and we leave it for future
work.

##  5 Related Work

In this section, we review two closely related methods which we use as
benchmarks and summarise the relationship between methods.

The first work is by Chung et al. (2022a), abbreviated DPS-D here222Since the
authors run DDPM-type sampler, we re-abbreviate DPS as DPS-D in this work (as
DPS approximation itself can be also run with Euler-Maruyama schemes). where
$p_{t}({\mathbf{x}}_{0}|{\mathbf{x}}_{t})$ is approximated by a Dirac delta
(point) distribution centred at ${\mathbf{m}}_{0|t}$ computed by Tweedie’s
formula, within the integral given in equation 5. The mean is chosen the same
way as our method. Authors choose a Dirac centred at this mean, in our
framework this corresponds to choosing a zero covariance, i.e.,

$\displaystyle{\mathbf{m}}_{0|t}^{\textup{DPS-D}}={\mathbf{m}}_{0|t}\quad\quad\text{and}\quad\quad{\mathbf{C}}_{0|t}^{\text{DPS-D}}=0.$
(13)

In the work Song et al. (2023), abbreviated $\Pi$GDM-D here333Since the
authors run DDIM-type sampler, we re-abbreviate $\Pi$GDM as $\Pi$GDM-D in this
work (as $\Pi$GDM approximation itself can also be run with Euler-Maruyama
schemes). the same estimator for the mean is chosen. However, the variance is
set to a multiple of the identity, corresponding to choices

$\displaystyle{\mathbf{m}}_{0|t}^{\text{$\Pi$G}}={\mathbf{m}}_{0|t}\quad\quad\text{and}\quad\quad{\mathbf{C}}_{0|t}^{\text{$\Pi$G}}=r_{t}^{2}{\mathbf{I}}_{d_{x}}$
(14)

The choice of $r_{t}$ is such that it matches the variance of the reverse SDE
if $p_{\rm{data}}$ would be a standard normal distribution. Since they employ
a different forward SDE to ours (a standard Brownian motion, also called
variance exploding SDE), $r_{t}$ is set to be equal to
$((v_{t})^{-1}+1)^{-1}$. In our case, with the OU/variance preserving SDE as a
forward process, $r_{t}$ would be equal to $v_{t}$.

Another relevant work is by Stevens et al. (2023) who consider maximum-a-
posteriori (MAP) approach to find the moments of $p_{0|t}$.

##  6 Experiments

In this section, we demonstrate our results as well as the peformance of other
approximations to the likelihood provided in Chung et al. (2022a); Song et al.
(2023). In particular, we perform comparisons for two of our methods TMPD (an
SGM using our approximation) and TMPD-D (a DDPM sampler using equation 11). We
compare these to DPS (an SGM sampler using the posterior approximation in
equation 13), DPS-D (Chung et al., 2022a) (a DDPM-type sampler using equation
13), $\Pi$GDM (Song et al., 2023) (an SGM sampler using the posterior
approximation in equation 14), and finally $\Pi$GDM-D (a DDIM-type sampler
using equation 14, but in our experiments, since we set the DDIM
hyperparameter $\eta=1.0$ this is equivalent to a DDPM-type sampler (Song et
al., 2021a)).

###  6.1 Gaussian Mixture Model

We now demonstrate a non-linear SDE example and follow the Gaussian Mixture
Model example of Cardoso et al. (2023) where the data distribution
$p_{0}({\mathbf{x}}_{0})$ is a mixture of 25 Gaussian distributions. The means
and variances of the components of the mixture are given in Appendix E.3. In
this case, for each choice of observation ${\mathbf{y}}$, observation map
${\mathbf{H}}$ and measurement noise standard deviation
$\sigma_{{\textnormal{y}}}$, the target posterior can be computed explicitly
(see Appendix E.3). To investigate the performance of posterior sampling
methods, for each pair of dimensions and observation noise
$(d_{{\textnormal{x}}},d_{{\textnormal{y}}},\sigma_{{\textnormal{y}}})\in\\{8,80,800\\}\times\\{1,2,4\\}\times\\{10^{-2},10^{-1},10^{0}\\}$
we randomly generate multiple measurement models
$({\mathbf{y}},{\mathbf{H}})\in\mathbb{R}^{d_{{\textnormal{y}}}}\times\mathbb{R}^{d_{{\textnormal{y}}}\times
d_{{\textnormal{x}}}},$ and equally weight each component of the Gaussian
mixture. Further details are given in Appendix E.3. We chose to control the
dimension to gain insight into the performance of posterior sampling methods
under varying dimensions. We chose to control and the noise level since the
different posterior sampling methods have accuracy that depends on the signal
to noise ratio. Through randomly varying the observation model, we gain an
insight into the performance of the posterior sampling methods with different
levels of posterior multimodality. This example is interesting because it
allows us to study the behaviour of our methods on non-linear problems in high
dimensions whilst having access to the target posterior with which to compare
(usually, obtaining a ‘ground-truth’ posterior is not feasible for non-linear
problems).

To compare the posterior distribution estimated by each algorithm with the
target posterior distribution, we use the sliced Wasserstein (SW) distance
defined in Appendix E.3. We use $10^{4}$ slices for the SW distance and
compare 1000 samples of TMPD-D, $\Pi$GDM-D and DPS-D in Tables 1 obtained
using 1000 denoising steps and 1000 samples of the true posterior
distribution.

Table 1 indicates the Central Limit Theorem (CLT) 95% confidence intervals
obtained by considering 20 randomly selected measurement models
(${\mathbf{H}}$) for each setting
($d_{{\textnormal{x}}},d_{{\textnormal{x}}},\sigma_{{\textnormal{y}}}$).
Figure 2 shows the first two dimensions of the estimated posterior
distributions corresponding to the configurations ($80,1$) from Table 1 for
one of the randomly generated measurement model
(${\mathbf{H}},\sigma_{{\textnormal{y}}}=0.1$). Illustration of other settings
are given in Appendix E.3. These illustrations give us insight into the
behaviour of the algorithms and their ability to accurately estimate the
posterior distribution. We see that TMPD-D estimates the target posterior well
compared to $\Pi$GDM-D and DPS-D. TMPD=D covers all of the modes, whereas
$\Pi$GDM-D and DPS-D do not.

We perform the same experiment using 1000 samples of TMPD, DTMPD, DPS and
$\Pi$GDM in Appendix E.3, obtained using 1000 Euler-Maruyama time-steps, and
results are shown in Appendix E.3.

Figure 2: We display the first two dimensions of the GMM inverse problem for one of the measurement models tested (${\mathbf{H}},\sigma_{{\textnormal{y}}}=0.1,(d_{{\textnormal{x}}},d_{{\textnormal{y}}})=(80,1)$). The blue dots represent samples from the target posterior, while the red dots correspond to samples generated by each of the algorithms used (the names of the algorithms are given at the bottom of each column). Table 1: Sliced Wasserstein for the GMM case. |  | $\sigma_{{\mathbf{y}}}=0.01$ | $\sigma_{{\mathbf{y}}}=0.1$ | $\sigma_{{\mathbf{y}}}=1.0$
---|---|---|---|---
$d_{x}$ | $d_{y}$ | TMPD-D |  $\Pi$GDM-D | DPS-D | TPMD-D |  $\Pi$GDM-D | DPS-D | TMPD-D |  $\Pi$GDM-D | DPS-D
$8$ | $1$ | 1.6 $\pm$ 0.5 | 2.6 $\pm$ 0.9 | 4.7 $\pm$ 1.5 | 1.4 $\pm$ 0.5 | 2.2 $\pm$ 0.9 | 4.7 $\pm$ 1.6 | 0.9 $\pm$ 0.3 | 1.5 $\pm$ 0.4 | 5.2 $\pm$ 1.3
$8$ | $2$ | 0.7 $\pm$ 0.3 | 2.1 $\pm$ 1.0 | 1.8 $\pm$ 1.5 | 0.9 $\pm$ 0.3 | 1.6 $\pm$ 0.6 | 1.5 $\pm$ 0.9 | 0.9 $\pm$ 0.2 | 2.3 $\pm$ 0.4 | 3.5 $\pm$ 1.2
$8$ | $4$ |  0.3 $\pm$ 0.3 | 3.8 $\pm$ 2.3 | 0.7 $\pm$ 0.6 |  0.3 $\pm$ 0.2 | 3.8 $\pm$ 2.2 | 0.8 $\pm$ 0.6 |  0.6 $\pm$ 0.2 | 1.8 $\pm$ 0.3 | 2.5 $\pm$ 0.9
$80$ | $1$ |  2.7 $\pm$ 0.7 | 3.2 $\pm$ 1.0 | 5.6 $\pm$ 1.8 |  2.3 $\pm$ 0.7 | 2.9 $\pm$ 0.8 | 5.1 $\pm$ 1.8 |  1.5 $\pm$ 0.7 | 1.6 $\pm$ 0.5 | 6.9 $\pm$ 1.8
$80$ | $2$ |  1.0 $\pm$ 0.5 | 2.8 $\pm$ 1.3 | 3.2 $\pm$ 1.9 |  1.2 $\pm$ 0.5 | 2.7 $\pm$ 1.2 | 3.1 $\pm$ 1.9 |  1.1 $\pm$ 0.2 | 1.4 $\pm$ 0.2 | 3.9 $\pm$ 1.2
$80$ | $4$ |  0.3 $\pm$ 0.1 | 0.6 $\pm$ 0.4 | 1.2 $\pm$ 1.1 |  0.4 $\pm$ 0.2 | 0.6 $\pm$ 0.4 | 1.0 $\pm$ 1.1 |  0.9 $\pm$ 0.3 |  0.9 $\pm$ 0.2 | 1.7 $\pm$ 0.6
$800$ | $1$ |  3.1 $\pm$ 0.7 | 3.5 $\pm$ 1.1 | 5.8 $\pm$ 1.6 |  2.9 $\pm$ 0.6 | 3.3 $\pm$ 0.9 | 5.7 $\pm$ 1.6 |  1.5 $\pm$ 0.4 | 2.0 $\pm$ 0.4 | 6.8 $\pm$ 1.0
$800$ | $2$ |  1.4 $\pm$ 0.4 | 3.1 $\pm$ 1.1 | 3.5 $\pm$ 1.7 |  1.3 $\pm$ 0.3 | 2.7 $\pm$ 0.9 | 3.1 $\pm$ 1.4 |  1.2 $\pm$ 0.3 | 2.0 $\pm$ 0.5 | 4.7 $\pm$ 1.3
$800$ | $4$ |  0.4 $\pm$ 0.2 |  0.4 $\pm$ 0.2 | 1.4 $\pm$ 1.0 |  0.4 $\pm$ 0.2 |  0.4 $\pm$ 0.2 | 1.3 $\pm$ 0.9 | 0.9 $\pm$ 0.2 |  0.6 $\pm$ 0.2 | 0.9 $\pm$ 0.4

###  6.2 Noisy observation inpainting

We consider inpainting problems on the FFHQ ($256\times 256\times 3$) (Karras
et al., 2019), CelebA-HQ ($256\times 256\times 3$) (Karras et al., 2017) and
CIFAR-10 ($32\times 32\times 3$) (Krizhevsky et al., 2009) datasets. We
compare TMPD to $\Pi$DGM and DPS. We also compare score-based diffusion models
with their denoising-diffusion counterparts (denoted with suffix, -D). We use
pretrained denoising networks that are available here for all the methods. The
methods TMPD $\Pi$GDM, DPS all have the same numerical solver of their
respective reverse-SDE, and TMPD-D $\Pi$GDM-D, DPS-D all use DDPM since DDIM
and DDPM are equivalent with our chosen hyperparameter $\eta=1.0$ (Song et
al., 2021a). Therefore the sampling methods being compared only differ in the
$\nabla_{{\mathbf{x}}_{t}}\log
p_{{\mathbf{y}}|t}({\mathbf{y}}|{\mathbf{x}}_{t})$ term in their reverse-SDE,
and so this experiment allows us to study the behaviour of our method on high-
dimensional image inpainting compared to the different approximations of the
smoothed likelihood. We use the ‘box’ mask on images from the CIFAR-10 in Fig.
3, and we use the ‘half’ mask on images from FFHQ and CelebA-HQ. Images are
normalized to the range $[0,1]$ and it is on this scale that we add Gaussian
measurement noise with standard deviation $\sigma_{y}\in\\{0.001,0.01\\}$.
Forward measurement operators are specified as follows: For half-type
inpainting, we mask out a $128\times 128$ (or $16\times 16$ for CIFAR-10)
right half region of the image; for box-type inpainting, we mask out a
$64\times 64$ (or $8\times 8$ for CIFAR-10) box region following Cardoso et
al. (2023). Whereas no hyperparameters are required for our method, more
experimental details including the hyper-parameters of DPS and $\Pi$GDM can be
found in Appendix E.4.

Figure 3: Results on solving inpainting problems with Gaussian noise. For more
details and examples, see Appendix E.4. Figure 4: Results on solving super-
resolution problems with Gaussian noise. For more details and examples, see
Appendix E.4.

###  6.3 Super-resolution

We consider noisy super-resolution problems on CIFAR-10 and CelebA-HQ. We set
$\sigma_{{\textnormal{y}}}\in\\{0.01,0.05,0.1\\}$ for all the datasets. The
results are shown in Fig. 4. For CIFAR-10, We use a downsampling ratio of 4
and a nearest-neighbour downsampling method. For CelebA-HQ dataset we use a
downsampling ratio of 16 and a bicubic downsampling method. More details
including the hyper-parameters of DPS and $\Pi$GDM can be found in Appendix
E.4.

##  7 Discussions, limitations and future work

In this paper, we introduced TMPD, a diffusion modelling approach to solve
inverse problems and sample from conditional distributions using unconditional
diffusion models. On various tasks, TMPD achieves competitive quality with
other methods that aim to solve the noisy, linear inverse problem while
avoiding the expensive, problem-specific training of conditional models. We
can use problem-agnostic diffusion models to solve certain problems that would
be cost-ineffective to address individually with conditional diffusion models,
leading to a much wider set of applications. The ability to handle measurement
noise also gives TMPD the potential to address certain real-world
applications, such as CT and MRI imaging with Gaussian noise (Sijbers & den
Dekker, 2004; Song et al., 2021b). Despite having on-par restoration results
with $\Pi$GDM (Song et al., 2023), TMPD is slower, as each iteration costs
more memory and compute due to the Jacobian over the score model. Therefore,
it would be helpful to explore more efficient sampling techniques, such as
ensemble methods or low rank approximations.

On the positive side, we found our approach more robust across various
problems due to the stabilising effect of better approximations. Moreover, our
approach is more principled compared to existing approaches, as shown in
Section 4. We provided a way to analyse similar methods and our moment
approximations can be analysed more rigorously to provide deeper theoretical
results for these kinds of methods. Our future work plans include expanding
the analysis we provided in this work.

## Acknowledgements

This work has been supported by The Alan Turing Institute through the Theory
and Methods Challenge Fortnights event “Accelerating generative models and
nonconvex optimisation”, which took place on 6-10 June 2022 and 5-9 Sep 2022
at The Alan Turing Institute headquarters. The authors would also like to
thank the Isaac Newton Institute for Mathematical Sciences, Cambridge, for
support and hospitality during the programme The Mathematical and Statistical
Foundation of Future Data-Driven Engineering where work on this paper was
undertaken. This work was supported by EPSRC grant no EP/R014604/1. JP and SR
have been funded by Deutsche Forschungsgemeinschaft (DFG) - Project-ID
318763901 - SFB1294. BB gratefully acknowledges the EPSRC for funding this
research through the EPSRC Centre for Doctoral Training in Future
Infrastructure and Built Environment: Resilience in a Changing World (EPSRC
grant reference number EP/S02302X/1); and nPlan for the invaluable
facilitation of work that was completed whilst on internship with nPlan.

## References

  * Batzolis et al. (2021) Georgios Batzolis, Jan Stanczuk, Carola-Bibiane Schönlieb, and Christian Etmann.  Conditional image generation with score-based diffusion models.  _arXiv preprint arXiv:2111.13606_ , 2021.
  * Bishop (2006) Christopher M Bishop.  _Pattern recognition and machine learning_ , volume 4.  Springer, 2006.
  * Bora et al. (2017) Ashish Bora, Ajil Jalal, Eric Price, and Alexandros G Dimakis.  Compressed sensing using generative models.  In _International conference on machine learning_ , pp. 537–546. PMLR, 2017.
  * Cardoso et al. (2023) Gabriel Cardoso, Yazid Janati El Idrissi, Sylvain Le Corff, and Eric Moulines.  Monte carlo guided diffusion for bayesian linear inverse problems, 2023\.
  * Chen et al. (2023) Sitan Chen, Sinho Chewi, Jerry Li, Yuanzhi Li, Adil Salim, and Anru R Zhang.  Sampling is as easy as learning the score: theory for diffusion models with minimal data assumptions.  In _International Conference on Learning Representations_ , 2023.
  * Chung et al. (2022a) Hyungjin Chung, Jeongsol Kim, Michael Thompson Mccann, Marc Louis Klasky, and Jong Chul Ye.  Diffusion posterior sampling for general noisy inverse problems.  In _The Eleventh International Conference on Learning Representations_ , 2022a.
  * Chung et al. (2022b) Hyungjin Chung, Byeongsu Sim, Dohoon Ryu, and Jong Chul Ye.  Improving diffusion models for inverse problems using manifold constraints.  _Advances in Neural Information Processing Systems_ , 35:25683–25696, 2022b.
  * Chung et al. (2023) Hyungjin Chung, Jeongsol Kim, Sehui Kim, and Jong Chul Ye.  Parallel diffusion models of operator and image for blind inverse problems.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pp. 6059–6069, 2023.
  * Daras et al. (2022) Giannis Daras, Yuval Dagan, Alex Dimakis, et al.  Score-guided intermediate level optimization: Fast langevin mixing for inverse problems.  In _Proceedings of the 39th International Conference on Machine Learning (ICML)_ , 2022.
  * Dhariwal & Nichol (2021) Prafulla Dhariwal and Alexander Nichol.  Diffusion models beat gans on image synthesis.  _Advances in neural information processing systems_ , 34:8780–8794, 2021.
  * Efron (2011) Bradley Efron.  Tweedie’s formula and selection bias.  _Journal of the American Statistical Association_ , 106(496):1602–1614, 2011.  doi: 10.1198/jasa.2011.tm11181.  URL https://doi.org/10.1198/jasa.2011.tm11181.  PMID: 22505788.
  * Feng & Bouman (2023) Berthy T Feng and Katherine L Bouman.  Efficient bayesian computational imaging with a surrogate score-based prior.  _arXiv preprint arXiv:2309.01949_ , 2023.
  * Feng et al. (2023) Berthy T Feng, Jamie Smith, Michael Rubinstein, Huiwen Chang, Katherine L Bouman, and William T Freeman.  Score-based diffusion models as principled priors for inverse imaging.  _arXiv preprint arXiv:2304.11751_ , 2023.
  * Givens & Shortt (1984) Clark R. Givens and Rae Michael Shortt.  A class of Wasserstein metrics for probability distributions.  _Michigan Mathematical Journal_ , 31(2):231 – 240, 1984.  doi: 10.1307/mmj/1029003026.  URL https://doi.org/10.1307/mmj/1029003026.
  * Ho et al. (2020) Jonathan Ho, Ajay Jain, and Pieter Abbeel.  Denoising diffusion probabilistic models, 2020.
  * Hyvärinen (2005) Aapo Hyvärinen.  Estimation of non-normalized statistical models by score matching.  _Journal of Machine Learning Research_ , 6(24):695–709, 2005.  URL http://jmlr.org/papers/v6/hyvarinen05a.html.
  * Jalal et al. (2021) Ajil Jalal, Marius Arvinte, Giannis Daras, Eric Price, Alexandros G Dimakis, and Jon Tamir.  Robust compressed sensing mri with deep generative priors.  _Advances in Neural Information Processing Systems_ , 34:14938–14954, 2021.
  * Kadkhodaie & Simoncelli (2021) Zahra Kadkhodaie and Eero Simoncelli.  Stochastic solutions for linear inverse problems using the prior implicit in a denoiser.  _Advances in Neural Information Processing Systems_ , 34:13242–13254, 2021.
  * Karras et al. (2017) Tero Karras, Timo Aila, Samuli Laine, and Jaakko Lehtinen.  Progressive growing of gans for improved quality, stability, and variation.  _arXiv preprint arXiv:1710.10196_ , 2017.
  * Karras et al. (2019) Tero Karras, Samuli Laine, and Timo Aila.  A style-based generator architecture for generative adversarial networks.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , pp. 4401–4410, 2019.
  * Kawar et al. (2021) Bahjat Kawar, Gregory Vaksman, and Michael Elad.  Snips: Solving noisy inverse problems stochastically.  _Advances in Neural Information Processing Systems_ , 34:21757–21769, 2021.
  * Kawar et al. (2022) Bahjat Kawar, Michael Elad, Stefano Ermon, and Jiaming Song.  Denoising diffusion restoration models.  _Advances in Neural Information Processing Systems_ , 35:23593–23606, 2022.
  * Krizhevsky et al. (2009) Alex Krizhevsky, Geoffrey Hinton, et al.  Learning multiple layers of features from tiny images.  2009\.
  * Laumont et al. (2022) Rémi Laumont, Valentin De Bortoli, Andrés Almansa, Julie Delon, Alain Durmus, and Marcelo Pereyra.  Bayesian imaging using plug & play priors: when Langevin meets Tweedie.  _SIAM Journal on Imaging Sciences_ , 15(2):701–737, 2022.
  * Mardani et al. (2023) Morteza Mardani, Jiaming Song, Jan Kautz, and Arash Vahdat.  A variational perspective on solving inverse problems with diffusion models.  _arXiv preprint arXiv:2305.04391_ , 2023.
  * Meng et al. (2021) Chenlin Meng, Yang Song, Wenzhe Li, and Stefano Ermon.  Estimating high order gradients of the data distribution by denoising.  _Advances in Neural Information Processing Systems_ , 34:25359–25369, 2021.
  * Mohamed et al. (2020) Shakir Mohamed, Mihaela Rosca, Michael Figurnov, and Andriy Mnih.  Monte carlo gradient estimation in machine learning.  _The Journal of Machine Learning Research_ , 21(1):5183–5244, 2020.
  * Pidstrigach et al. (2023) Jakiw Pidstrigach, Youssef Marzouk, Sebastian Reich, and Sven Wang.  Infinite-dimensional diffusion models for function spaces.  _arXiv preprint arXiv:2302.10130_ , 2023.
  * Robbins (1992) Herbert E. Robbins.  _An Empirical Bayes Approach to Statistics_ , pp. 388–394.  Springer New York, New York, NY, 1992.  ISBN 978-1-4612-0919-5.  doi: 10.1007/978-1-4612-0919-5˙26.  URL https://doi.org/10.1007/978-1-4612-0919-5_26.
  * Sijbers & den Dekker (2004) J. Sijbers and A.J. den Dekker.  Maximum likelihood estimation of signal amplitude and noise variance from mr data.  _Magnetic Resonance in Medicine_ , 51(3):586–594, 2004.  doi: https://doi.org/10.1002/mrm.10728.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/mrm.10728.
  * Song et al. (2021a) Jiaming Song, Chenlin Meng, and Stefano Ermon.  Denoising diffusion implicit models.  2021a.  URL https://openreview.net/forum?id=St1giarCHLP.
  * Song et al. (2023) Jiaming Song, Arash Vahdat, Morteza Mardani, and Jan Kautz.  Pseudoinverse-guided diffusion models for inverse problems.  In _International Conference on Learning Representations_ , 2023.  URL https://openreview.net/forum?id=9_gsMA8MRKQ.
  * Song et al. (2020) Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole.  Score-based generative modeling through stochastic differential equations.  _arXiv preprint arXiv:2011.13456_ , 2020.
  * Song et al. (2021b) Yang Song, Liyue Shen, Lei Xing, and Stefano Ermon.  Solving inverse problems in medical imaging with score-based generative models.  _arXiv preprint arXiv:2111.08005_ , 2021b.
  * Stevens et al. (2023) Tristan SW Stevens, Jean-Luc Robert, Faik C Yu, Jun Seob Shin, and Ruud JG van Sloun.  Removing structured noise with diffusion models.  _arXiv preprint arXiv:2302.05290_ , 2023.
  * Stuart (2010) Andrew M Stuart.  Inverse problems: a bayesian perspective.  _Acta numerica_ , 19:451–559, 2010.
  * Tarantola (2005) Albert Tarantola.  _Inverse problem theory and methods for model parameter estimation_.  SIAM, 2005.

##  Appendix A Proofs for Section 3

###  A.1 Proof of Proposition 1

In order to prove this result, we adapt Theorem 1 of Meng et al. (2021). We
write the proof for generic exponential family which can be adapted to our
case easily. Let us consider

$\displaystyle
p_{t|0}({\mathbf{x}}_{t}|{\mathbf{x}}_{0})=\mathcal{N}({\mathbf{x}}_{t};\sqrt{\alpha_{t}}{\mathbf{x}}_{0},v_{t}{\mathbf{I}}_{d_{x}}).$

and write $p_{t}({\mathbf{x}}_{t})=\int
p_{t|0}({\mathbf{x}}_{t}|{\mathbf{x}}_{0})p_{0}({\mathbf{x}}_{0})\mathrm{d}{\mathbf{x}}_{0}$.
We are interested in finding the moments of the posterior

$\displaystyle
p_{0|t}({\mathbf{x}}_{0}|{\mathbf{x}}_{t})=\frac{p_{t|0}({\mathbf{x}}_{t}|{\mathbf{x}}_{0})p_{0}({\mathbf{x}}_{0})}{p_{t}({\mathbf{x}}_{t})}.$

Let us redefine the right handside in this Bayes’ rule using the exponential
family parameterisation of the Gaussian
$p_{t|0}({\mathbf{x}}_{t}|{\mathbf{x}}_{0})$

$\displaystyle
p(\bm{\eta}_{0}|{\mathbf{x}}_{t})=\frac{p_{t|0}({\mathbf{x}}_{t}|\bm{\eta}_{0})q_{0}(\bm{\eta}_{0})}{p_{t}({\mathbf{x}}_{t})},$

where $\bm{\eta}_{0}={\mathbf{x}}_{0}\sqrt{\alpha}_{t}/v_{t}$ and

$\displaystyle
p_{t|0}({\mathbf{x}}_{t}|\bm{\eta}_{0})=e^{\bm{\eta}_{0}^{\top}{\mathbf{x}}_{t}-\psi(\bm{\eta}_{0})}q({\mathbf{x}}_{t}),$

where

$\displaystyle
q({\mathbf{x}}_{t})=((2\pi)^{d}v^{2d})^{-1/2}e^{-{\mathbf{x}}_{t}^{\top}{\mathbf{x}}_{t}/2v_{t}}$

Now let

$\displaystyle\lambda({\mathbf{x}}_{t})=\log\frac{p_{t}({\mathbf{x}}_{t})}{q({\mathbf{x}}_{t})},$

we can rewrite

$\displaystyle
p(\bm{\eta}_{0}|{\mathbf{x}}_{t})=e^{\bm{\eta}_{0}^{\top}{\mathbf{x}}_{t}-\psi({\mathbf{x}}_{t})-\lambda({\mathbf{x}}_{t})}q_{0}(\bm{\eta}_{0}).$

It is then easy to show that the moments of $\bm{\eta}_{0}$ are given by (Meng
et al., 2021)

$\displaystyle\mathbb{E}[\bm{\eta}_{0}|{\mathbf{x}}_{t}]$
$\displaystyle=\mathbf{J}_{\lambda}({\mathbf{x}}_{t}),$
$\displaystyle\mathbb{E}[\bm{\eta}_{0}^{\top}|{\mathbf{x}}_{t}]$
$\displaystyle=\mathbf{J}_{\lambda}({\mathbf{x}}_{t})^{\top},$
$\displaystyle\mathbb{E}[\bm{\eta}_{0}\bm{\eta}_{0}^{\top}|{\mathbf{x}}_{t}]$
$\displaystyle=\mathbf{S}_{\lambda}({\mathbf{x}}_{t})+\mathbf{J}_{\lambda}({\mathbf{x}}_{t})\mathbf{J}_{\lambda}({\mathbf{x}}_{t})^{\top},$

where $\mathbf{J}_{\lambda}$ is the Jacobian of $\lambda$ w.r.t.
${\mathbf{x}}_{t}$ and $\mathbf{S}_{\lambda}$ is the Hessian of $\lambda$
w.r.t. ${\mathbf{x}}_{t}$. Recall now that
$\bm{\eta}_{0}={\mathbf{x}}_{0}\sqrt{\alpha_{t}}/v_{t}$ and

$\displaystyle\lambda({\mathbf{x}}_{t})$ $\displaystyle=\log
p_{t}({\mathbf{x}}_{t})+\frac{{\mathbf{x}}_{t}^{\top}{\mathbf{x}}_{t}}{2v_{t}}+C$
$\displaystyle\mathbf{J}_{\lambda}({\mathbf{x}}_{t})$
$\displaystyle=\nabla_{{\mathbf{x}}_{t}}\log
p_{t}({\mathbf{x}}_{t})+\frac{{\mathbf{x}}_{t}}{v_{t}}$
$\displaystyle\mathbf{S}_{\lambda}({\mathbf{x}}_{t})$
$\displaystyle=\nabla^{2}_{{\mathbf{x}}_{t}}\log
p_{t}({\mathbf{x}}_{t})+\frac{1}{v_{t}}{\mathbf{I}}_{d_{x}}.$

This implies that

$\displaystyle\mathbb{E}[{\mathbf{x}}_{0}|{\mathbf{x}}_{t}]=\frac{1}{\sqrt{\alpha_{t}}}({\mathbf{x}}_{t}+v_{t}\nabla_{{\mathbf{x}}_{t}}\log
p_{t}({\mathbf{x}}_{t})),$

which proves the Tweedie’s formula for the mean. For the covariance, note that

$\displaystyle\text{Cov}(\bm{\eta_{0}}|{\mathbf{x}}_{t})$
$\displaystyle=\mathbb{E}[\bm{\eta}_{0}\bm{\eta}_{0}^{\top}|{\mathbf{x}}_{t}]-\mathbb{E}[\bm{\eta}_{0}|{\mathbf{x}}_{t}]\mathbb{E}[\bm{\eta}_{0}|{\mathbf{x}}_{t}]^{\top}$
$\displaystyle=\nabla^{2}\log
p_{t}({\mathbf{x}}_{t})+\frac{1}{v_{t}}{\mathbf{I}}_{d_{x}}.$

Since

$\displaystyle\text{Cov}(\bm{\eta_{0}}|{\mathbf{x}}_{t})=\frac{\alpha_{t}}{v_{t}^{2}}\text{Cov}({\mathbf{x}}_{0}|{\mathbf{x}}_{t}),$

we conclude

$\displaystyle\text{Cov}({\mathbf{x}}_{0}|{\mathbf{x}}_{t})=\frac{v_{t}}{\alpha_{t}}({\mathbf{I}}_{d_{x}}+v_{t}\nabla^{2}\log
p_{t}({\mathbf{x}}_{t})),$

which concludes the proof. $\square$

##  Appendix B Proofs for Section 4

###  B.1 Proof of Proposition 3

If $p_{\rm{data}}$ is Gaussian, then the full process ${\mathbf{x}}_{t}$ is a
Gaussian process. In particular, $({\mathbf{x}}_{0},{\mathbf{x}}_{t})$ are
jointly Gaussian:

$\displaystyle{\mathbf{x}}_{0}|{\mathbf{x}}_{t}\propto\mathcal{N}(\mathbb{E}[{\mathbf{x}}_{0}|{\mathbf{x}}_{t}],\mathbb{V}[{\mathbf{x}}_{0}|{\mathbf{x}}_{t}]),$

where $\mathbb{V}$ denotes the covariance. However, the right-hand side is
precisely the approximation we make, due to Proposition 1. Therefore, the
approximation we make in equation 9 is correct. Adding this to the learned
drift $\mathsf{s}_{\theta}({\mathbf{x}}_{t},t)=\nabla\log
p_{t}({\mathbf{x}}_{t})$, we get an expression for $\nabla\log
p_{t}({\mathbf{x}}_{t}|{\mathbf{y}})$ by equation 4.

###  B.2 Proof of Theorem 1

We first introduce three SDEs.

##### Conditioned SDEs

: The first pair of SDEs is given by an OU process, but started in the
conditional distribution:

$\mathrm{d}{\mathbf{x}}_{t}^{c}=-\frac{1}{2}{\mathbf{x}}_{t}^{c}\mathrm{d}t+\mathrm{d}{\mathbf{w}}_{t},\quad{\mathbf{x}}_{0}^{c}\sim
p_{0}^{c}:=p_{\rm{data}}({\mathbf{x}}_{0}|{\mathbf{y}}).$ (15)

We denote the marginals of ${\mathbf{x}}_{t}^{c}$ by $p_{t}^{c}$ to
differentiate them from the marginals $p_{t}$ of the OU process started in the
correct distribution, defined in equation 1. Note that
$p_{t}^{c}=p_{t}({\mathbf{x}}_{t}|{\mathbf{y}})$. The time reversal
${\mathbf{z}}_{t}={\mathbf{x}}_{T-t}$ then satisfies the SDE

$\displaystyle\mathrm{d}{\mathbf{z}}_{t}=$
$\displaystyle\frac{1}{2}{\mathbf{z}}_{t}\mathrm{d}t+\nabla_{{\mathbf{z}}_{t}}\log
p^{c}_{T-t}({\mathbf{z}}_{t})\mathrm{d}t+\mathrm{d}{\mathbf{w}}_{t}$ (16)
$\displaystyle=$
$\displaystyle\frac{1}{2}{\mathbf{z}}_{t}\mathrm{d}t+\nabla_{{\mathbf{z}}_{t}}\log
p_{T-t}({\mathbf{z}}_{t})\mathrm{d}t+\nabla_{{\mathbf{z}}_{t}}\log
p_{{\mathbf{y}}|T-t}({\mathbf{y}}|{\mathbf{z}}_{t})\mathrm{d}t+\mathrm{d}{\mathbf{w}}_{t}$
$\displaystyle\quad{\mathbf{z}}_{0}\sim$ $\displaystyle p^{c}_{T},$

where we used equation 4.

Solutions to the above reverse SDE will sample our target measure
$p_{0}^{c}=p_{\rm{data}}({\mathbf{x}}_{0}|{\mathbf{y}})$ at final time.
Therefore, we want to study how our algorithm approximates solutions of
${\mathbf{x}}^{c}$.

##### Intermediate Gaussian SDE:

Instead of bounding the distance of solutions to the above conditioned reverse
SDE to our algorithm, we will instead introduce an intermediate process, which
we will later use in a triangle inequality. The process will be a Gaussian
process. Therefore, we denote it with a superscript $G$.

The process ${\mathbf{x}}^{G}_{t}$ will be defined analogous to equation 1,
but assuming that it is started in
$\mathcal{N}({\mathbf{x}}_{0};\bm{\mu}_{0},\bm{\Sigma}_{0})$ instead of
$p_{\rm{data}}$. Since the forward SDE in equation 1 is linear, all of the
marginals of ${\mathbf{x}}_{t}^{G}$, which we denote by $p_{t}^{G}$, will also
be Gaussian.

Again, we define a conditioned version of ${\mathbf{x}}$, called
${\mathbf{x}}^{G,c}$ analogous to equation 15, i.e. ${\mathbf{x}}_{0}^{G,c}$
will be distributed as ${\mathbf{x}}_{0}^{G}$ conditioned on ${\mathbf{y}}=y$,
i.e. $p^{G}({\mathbf{x}}_{0}|{\mathbf{y}}=y)$. Since we have a linear
observation model, $p^{G}({\mathbf{x}}_{0}|{\mathbf{y}}=y)$ is still a
Gaussian, and therefore also ${\mathbf{x}}^{G,c}$ will be a Gaussian process.
Its reverse SDE ${\mathbf{z}}^{G,c}$ is defined through equation 16, just that
every appearance of $p$ will also have a superscript $G$,

$\displaystyle\mathrm{d}{\mathbf{z}}_{t}=$
$\displaystyle\frac{1}{2}{\mathbf{z}}_{t}\mathrm{d}t+\nabla_{{\mathbf{z}}_{t}}\log
p^{G}_{T-t}({\mathbf{z}}_{t})\mathrm{d}t+\nabla_{{\mathbf{z}}_{t}}\log
p^{G}_{T-t}({\mathbf{y}}|{\mathbf{z}}_{t})\mathrm{d}t+\mathrm{d}{\mathbf{w}}_{t}$
(17) $\displaystyle\quad{\mathbf{z}}_{0}\sim$ $\displaystyle p^{c}_{T},$

##### Algorithm SDE:

Finally, we define a third reverse SDE, which is the SDE that we are
discretizing when implementing our algorithm. It is given by

$\displaystyle\mathrm{d}{\mathbf{z}}_{t}=$
$\displaystyle\frac{1}{2}{\mathbf{z}}_{t}\mathrm{d}t+\nabla_{{\mathbf{z}}_{t}}\log
p_{T-t}({\mathbf{z}}_{t})\mathrm{d}t+f_{T-t}({\mathbf{z}}_{t})\mathrm{d}t+\mathrm{d}{\mathbf{w}}_{t}$
(18) $\displaystyle\quad{\mathbf{z}}_{0}\sim$ $\displaystyle p^{c}_{T},$

where

$f_{t}({\mathbf{z}}_{t})=\nabla_{{\mathbf{z}}_{t}}\mathbb{E}[{\mathbf{x}}_{0}|{\mathbf{x}}_{t}=z_{t}]{\mathbf{H}}^{T}\mathbb{V}[{\mathbf{x}}_{0}|{\mathbf{x}}_{t}=z_{t}]^{-1}({\mathbf{y}}-{\mathbf{H}}\mathbb{E}[{\mathbf{x}}_{0}|{\mathbf{x}}_{t}=z_{t}]),$
(19)

and $\mathbb{V}$ denotes the conditional covariance. Except for approximation
errors due to time discretization and in the initial conditions, this is the
SDE we are sampling from in our algorithm. We denote the marginal of
${\mathbf{z}}_{t}$ by $q_{t}$. Now we are ready to write our proof.

###### Proof of Theorem 1.

We have that $p_{\rm{data}}(\cdot|{\mathbf{y}})=p^{c}_{0}$

$\|p^{c}_{0}-q_{T}\|_{\text{TV}}\leq\|p^{c}_{0}-p^{G,c}_{0}\|_{\text{TV}}+\|p^{G,c}_{0}-q_{T}\|_{\text{TV}}$

##### Bounding the first term:

The first term is the total variation distance between
$p_{0}({\mathbf{x}}_{0}|{\mathbf{y}}=y)$ and
$p^{G}_{0}({\mathbf{x}}_{0}^{G}|{\mathbf{y}}=y)$. Now,

$\displaystyle p_{0}({\mathbf{x}}_{0}=x_{0}|{\mathbf{y}}=y)$
$\displaystyle=\frac{p_{{\mathbf{y}}|0}({\mathbf{y}}=y|{\mathbf{x}}_{0}=x_{0})p_{0}({\mathbf{x}}_{0}=x_{0})}{p_{{\mathbf{y}}}({\mathbf{y}}=y)}$
$\displaystyle=\frac{p^{G}_{{\mathbf{y}}|0}({\mathbf{y}}=y|{\mathbf{x}}_{0}=x_{0})p^{G}_{0}({\mathbf{x}}_{0}=x_{0})}{p^{G}_{{\mathbf{y}}}({\mathbf{y}}=y)}\frac{p_{0}({\mathbf{x}}_{0}=x_{0})}{p_{0}^{G}({\mathbf{x}}_{0}^{G}=x_{0}^{G})}\frac{p_{{\mathbf{y}}}({\mathbf{y}}=y)}{p^{G}_{{\mathbf{y}}}({\mathbf{y}}=y)}$
$\displaystyle=p_{{\mathbf{y}}|0}^{G}({\mathbf{x}}_{0}=x_{0}|{\mathbf{y}}=y)\exp(\Phi(x_{0}))\frac{p_{{\mathbf{y}}}({\mathbf{y}}=y)}{p_{{\mathbf{y}}}^{G}({\mathbf{y}}=y)},$

where we used that the conditional distribution of ${\mathbf{y}}$ given
${\mathbf{x}}_{0}$ does not depend on the distribution of ${\mathbf{x}}_{0}$.
For the last term we get

$\displaystyle p_{{\mathbf{y}}}({\mathbf{y}}=y)$ $\displaystyle=\int
p_{{\mathbf{y}}|0}({\mathbf{y}}=y|{\mathbf{x}}_{0}=x_{0})p_{0}({\mathbf{x}}_{0}=x_{0})\mathrm{d}x_{0}$
$\displaystyle=\int
p^{G}_{{\mathbf{y}}|0}({\mathbf{y}}=y|{\mathbf{x}}_{0}=x_{0})p^{G}_{0}({\mathbf{x}}_{0}=x_{0})\exp(\Phi(x_{0}))\mathrm{d}x_{0}$
$\displaystyle=Np^{G}_{{\mathbf{y}}}({\mathbf{y}}=y),$

with $N\in[1/M,M]$. Therefore,

$\displaystyle\|p^{c}_{0}-p^{G,c}_{0}\|_{\text{TV}}\leq\int|\frac{p^{c}_{0}}{p^{G,c}_{0}}(x_{0})-1|p^{G,c}_{0}(x_{0})\mathrm{d}x_{0}\leq
M^{2}-1,$

where we used that $M-1$ is always greater or equal to $1-\frac{1}{M}$, since
$M\geq 1$.

##### Bounding the second term:

For the second term we get that $p_{0}^{G,c}$ is the final marginal of
equation 17, while $q_{T}$ is the final marginal of equation 18. We now use
Pinsker’s inequality,

$\|p^{G,c}_{0}-q_{T}\|_{\text{TV}}\leq\sqrt{\text{KL}(p^{G,c}_{0}\|q_{T})}.$

We denote the path measures induced by equation 17 and equation 18 by
$\mathbb{P}^{G,c}$ and $\mathbb{Q}$ respectively. They have final time
marginals $p_{0}^{G,c}$ and $q_{T}$. Therefore, we can bound the KL-Divergence
of the marginals $p^{G,c}$ and $q_{T}$ by the KL-Divergence on the full path
space:

$\text{KL}(p^{G,c}_{0}\|q_{T})\leq\text{KL}(\mathbb{P}^{G,c}\|\mathbb{Q}).$

We will assume that we can apply Girsanovs theorem and show later that this is
justified. Using Girsanovs theorem, we can evaluate the Radon-Nikodym
derivative $\frac{\mathrm{d}\mathbb{Q}}{\mathrm{d}\mathbb{P}^{G,c}}$ on the
path space and therefore calculate the KL-divergence:

$\displaystyle\text{KL}(\mathbb{P}^{G,c}\|\mathbb{Q})=$
$\displaystyle\mathbb{E}_{\mathbb{P}^{G,c}}[\log\frac{\mathrm{d}\mathbb{P}^{G,c}}{\mathrm{d}\mathbb{Q}}]$
$\displaystyle=$
$\displaystyle\mathbb{E}_{\mathbb{P}^{G,c}}\left[-\int_{0}^{t}\nabla\log
p_{t}(x_{t})+f_{T-t}({\mathbf{x}}_{t})-\nabla\log
p_{T-t}^{G,c}({\mathbf{x}}_{t})\mathrm{d}w_{t}\right]$ $\displaystyle+$
$\displaystyle\mathbb{E}_{\mathbb{P}^{G,c}}\left[\int_{0}^{t}\|\nabla\log
p_{t}(x_{t})+f_{T-t}({\mathbf{x}}_{t})-\log
p_{T-t}^{G,c}({\mathbf{x}}_{t})\|^{2}\mathrm{d}t\right]$ $\displaystyle=$
$\displaystyle\int_{0}^{t}\mathbb{E}_{p_{t}^{G,c}}[\|\nabla\log
p_{t}(x_{t})+f_{T-t}({\mathbf{x}}_{t})-\log
p_{T-t}^{G,c}({\mathbf{x}}_{t})\|^{2}]\mathrm{d}t,$

where the term on the second line drops out because stochastic integrals have
expectation $0$. Now we have the drift of the Gaussian SDE given by

$\displaystyle\nabla\log p_{t}^{G,c}(x_{t})$ (20) $\displaystyle=$
$\displaystyle\nabla\log p_{t}^{G}(x_{t})+\nabla\log
p^{G}({\mathbf{y}}|{\mathbf{x}}_{t}^{G}=x_{t})$ $\displaystyle=$
$\displaystyle\nabla\log p_{t}^{G}(x_{t})$ $\displaystyle+$
$\displaystyle\nabla_{x_{t}}\mathbb{E}[{\mathbf{x}}^{G}_{0}|{\mathbf{x}}_{t}^{G}=x_{t}]{\mathbf{H}}^{\top}({\mathbf{H}}\mathbb{V}[{\mathbf{x}}^{G}_{0}|{\mathbf{x}}^{G}_{t}=x_{t}]{\mathbf{H}}^{T}+\sigma_{y}^{2}{\mathbf{I}}_{d_{y}})^{-1}({\mathbf{y}}-{\mathbf{H}}\mathbb{E}[{\mathbf{x}}^{G}_{0}|{\mathbf{x}}_{t}^{G}=x_{t}])$
$\displaystyle=:$ $\displaystyle\nabla\log
p_{t}^{G}(x_{t})+\tilde{f}_{t}(x_{t}),$

see equation 43, while the drift of the algorithm SDE is given by

$\displaystyle\nabla\log p_{t}(x_{t})+f_{t}(x_{t}),$ (21)

where $f_{t}$ is given by equation 19. We see that the difference in the SDE
drifts mainly consists of differences between conditional moments of $p_{t}$
and $p_{t}^{G}$, as well as the derivatives of the conditional expectations.
Therefore, the main difficulty of the proof is to bound these.

The density of the conditional distribution of $p_{0|t}$ is given by

$\displaystyle p_{0|t}(x_{0}|x_{t})=$
$\displaystyle\frac{p_{0,t}(x_{0},x_{t})}{p_{t}(x_{t})}=\frac{p^{G}_{0,t}(x_{0},x_{t})\exp(\Phi(x_{0}))}{p(x_{t})}=\frac{p^{G}_{0|t}(x_{0}|x_{t})p^{G}(x_{t})\exp(\Phi(x_{0}))}{p(x_{t})}$
(22) $\displaystyle=$ $\displaystyle
p^{G}_{0|t}\exp(\Phi(x_{0})-\Phi(x_{t})),$

where

$\displaystyle\exp(\Phi_{t}(x_{t}))=\frac{p_{t}(x_{t})}{p^{G}_{t}(x_{t})}=\mathbb{E}[\exp(\Phi({\mathbf{x}}_{0}))|{\mathbf{x}}_{t}=x_{t}].$
(23)

In the last equality, we used that

$\frac{\mathrm{d}\mathbb{P}}{\mathrm{d}\mathbb{P}^{G}}(x_{[0,t]})=\exp(\Phi(x_{0})),$

where we denoted by $\mathbb{P}$ the path measure induced by equation 15.
Therefore, also their marginals $p_{t}^{G}=\mathbb{P}^{G}_{t}$ and
$p_{t}=\mathbb{P}_{t}$ have relative densities, which are given by integrating
out the density to time $t$, as we did in equation 23.

By assumption $\exp(\Phi)$ is bounded from above and below by $M$ and $1/M$
respectively, and by equation 23 the same holds for $\exp(\Phi_{t})$.
Therefore, by equation 22, $p_{t|0}$ is absolutely continuous with respect to
$p_{t|0}^{G}$ with a density that is bounded from above and below by $M^{2}$
and $1/M^{2}$ respectively. We now obtain

$\displaystyle|\mathbb{E}[{\mathbf{x}}_{0}|{\mathbf{x}}_{t}=x_{t}]-\mathbb{E}[{\mathbf{x}}_{0}^{G}|{\mathbf{x}}_{t}^{G}=x_{t}]|$
$\displaystyle=$ $\displaystyle\left|\int
x_{0}[p(x_{0}|x_{t})-p^{G}(x_{0}|x_{t})]\mathrm{d}x_{0}\right|=\left|\int
x_{0}p^{G}(x_{0}|x_{t})(\exp(\Phi(x_{0})-\Phi(x_{t}))-1)\mathrm{d}x_{0}\right|$
$\displaystyle\leq$
$\displaystyle(M^{2}-1)\mathbb{E}[|{\mathbf{x}}_{0}^{G}||{\mathbf{x}}_{t}^{G}=x_{t}^{G}]$

The same holds for every entry of the covariance matrix. We denote

$\displaystyle V_{ij}$
$\displaystyle:=\mathbb{V}[{\mathbf{x}}_{0}|{\mathbf{x}}_{t}]_{ij},$
$\displaystyle V^{G}_{ij}$
$\displaystyle:=\mathbb{V}[{\mathbf{x}}_{0}^{G}|{\mathbf{x}}_{t}^{G}]_{ij}=N_{ij},$
$\displaystyle|V_{ij}-V^{G}_{ij}|$ $\displaystyle\leq(M^{2}-1)N_{ij}$
$\displaystyle N_{ij}$
$\displaystyle:=\mathbb{E}[|{\mathbf{x}}^{G}_{0,i}{\mathbf{x}}^{G}_{0,j}|~{}|{\mathbf{x}}_{t}^{G}]$

Now we get that

$\displaystyle\|V-V^{G}\|_{F}$
$\displaystyle=\left(\sum_{ij}(V_{ij}-V^{G}_{ij})^{2}\right)^{1/2}\leq(M^{2}-1)\|V^{G}\|_{F}\max_{ij}(N_{ij}).$

We can bound $N_{ij}$ by

$N_{ij}\leq\mathbb{E}[({\mathbf{x}}_{0,i}^{G})^{2}+({\mathbf{x}}_{0,j}^{G})^{2}|{\mathbf{x}}_{t}^{G}]\leq\text{Tr}(\mathbb{V}[{\mathbf{x}}_{0}^{G}|{\mathbf{x}}_{t}^{G}=x_{t}]).$
(24)

The latter term does not actually depend on $x_{t}$, but only on $t$. By using
the formulas in equation 39 we see thta it can be bounded independently of $t$
(only depending on $m_{0}$ and $\Sigma_{0}$). Therefore, we can put it into a
constant. We get that

$\|V-V^{G}\|_{F}\lesssim(M^{2}-1)\|V^{G}\|_{F},$

where

$a\lesssim b\quad\Leftrightarrow\quad a\leq Cb,$

with a constant $C$ only depending on $m_{0}$, $\Sigma_{0}$, $H$, $\sigma^{y}$
and the observation $y$. Now let $v$ be a vector with $\|v\|=1$. It holds that

$\displaystyle v^{T}V_{ij}v=$
$\displaystyle\int(v^{T}(x_{0}-\mathbb{E}[{\mathbf{x}}_{0}|{\mathbf{x}}_{t}=x_{t}]))^{2}p^{G}_{0|t}(x_{0}|x_{t})\exp(\Phi(x_{0})-\Phi_{t}(x_{t}))\mathrm{d}x_{0}$
(25) $\displaystyle=$
$\displaystyle\frac{1}{M}v^{T}V_{ij}^{G}v+N\int(v^{T}(\mathbb{E}[{\mathbf{x}}_{0}^{G}|{\mathbf{x}}_{t}^{G}=x_{t}]-\mathbb{E}[{\mathbf{x}}_{0}|{\mathbf{x}}_{t}=x_{t}]))^{2}p_{0|t}(x_{0}|x_{t})\mathrm{d}x_{0}$
$\displaystyle\geq$ $\displaystyle v^{T}V_{ij}^{G}v\frac{1}{M}$

with $N\geq 1/M$. Since $V$ and $V^{G}$ are positive semidefinite, symmetric
matrices, this implies that all eigenvalues of $V$ are bounded by the lowest
eigenvalue of $V^{G}$ times $1/M$. We define

$\displaystyle T_{1}$ $\displaystyle=HV^{G}H^{T}+\sigma_{y}^{2}I_{d_{y}}$ (26)
$\displaystyle T_{2}$ $\displaystyle=HVH^{T}+\sigma_{y}^{2}I_{d_{y}}.$ (27)

We then need to bound

$\|T_{1}^{-1}-T_{2}^{-1}\|\leq\|T_{1}-T_{2}\|(\|T_{1}^{-1}\|+\|T_{2}^{-1}\|).$

We start with

$\displaystyle\|T_{1}-T_{2}\|_{\text{op}}\leq\|H\|^{2}\|V-V^{G}\|_{\text{op}}\lesssim\|V-V^{G}\|_{F}\lesssim(M^{2}-1)\|V^{G}\|_{F}.$

We also used the equivalence of the Frobenius and operator norm here. Due to
our eigenvalue bound (equation 25) on $V$, we also get an analogous bound on
$V+\sigma_{y}I_{d_{y}}$. Therefore,

$\|T_{1}^{-1}-T_{2}^{-1}\|_{\text{op}}\lesssim(M^{2}-1)\|V^{G}\|_{F}\|V^{G}\|_{F}^{-1}(1+M)\leq(M^{2}-1)(M+1)\leq
M^{3}-1$

where we used that the operator norm of the inverse, is the inverse of the
operator norm and the equivalence of the operator norm to the Frobenius norm.

Finally, we need to bound

$\displaystyle\nabla_{x_{t}}\int x_{0}p(x_{0}|x_{t})\mathrm{d}x_{0}$ (28)
$\displaystyle=$ $\displaystyle\int
x_{0}\nabla_{x_{t}}p^{G}(x_{0}|x_{t})\exp(\Phi(x_{0})-\Phi(x_{t}))\mathrm{d}x_{0}$
$\displaystyle=$ $\displaystyle\int
x_{0}\exp(\Phi(x_{0})-\Phi(x_{t}))\nabla_{x_{t}}p^{G}(x_{0}|x_{t})\mathrm{d}x_{0}+\int
x_{0}\nabla_{x_{t}}\Phi(x_{t})\exp(\Phi(x_{0})-\Phi(x_{t}))p^{G}(x_{0}|x_{t})\mathrm{d}x_{0}$
$\displaystyle=$ $\displaystyle\int
x_{0}\exp(\Phi(x_{0})-\Phi(x_{t}))\nabla_{x_{t}}p^{G}(x_{0}|x_{t})\mathrm{d}x_{0}+\nabla_{x_{t}}\Phi(x_{t})\mathbb{E}[{\mathbf{x}}_{0}|{\mathbf{x}}_{t}=x_{t}].$

Also

$\displaystyle\left\|\int
x_{0}\exp(\Phi(x_{0})-\Phi(x_{t}))\nabla_{x_{t}}p^{G}(x_{0}|x_{t})\mathrm{d}x_{0}-\nabla_{x_{t}}\mathbb{E}[{\mathbf{x}}_{0}^{G}|{\mathbf{x}}_{t}^{G}=x_{t}]\right\|$
(29) $\displaystyle=$
$\displaystyle\left\|\int\nabla_{x_{t}}p^{G}(x_{0}|x_{t})x_{0}^{T}(\exp(\Phi(x_{0})-\Phi(x_{t}))-1)\mathrm{d}x_{0}\right\|$
$\displaystyle\leq$ $\displaystyle(M^{2}-1)\int p^{G}(x_{0}|x_{t})\|\nabla\log
p^{G}(x_{0}|x_{t})x_{0}^{T}\|\mathrm{d}x_{0}$ $\displaystyle\leq$
$\displaystyle(M^{2}-1)(\int p^{G}(x_{0}|x_{t})\|\nabla\log
p^{G}(x_{0}|x_{t})x_{0}^{T}\|\mathrm{d}x_{0})^{1/2}$ $\displaystyle\leq$
$\displaystyle(M^{2}-1)(\int
p^{G}(x_{0}|x_{t})x_{0}^{T}\Sigma_{0|t}^{-1}(x_{0}-m_{0|t}^{x_{t}})\mathrm{d}x_{0})^{1/2}$
$\displaystyle\lesssim$ $\displaystyle(M^{2}-1)(\int
p^{G}(x_{0}|x_{t})(x_{0}-m_{0|t}^{x_{t}})^{T}\Sigma_{0|t}^{-1}(x_{0}-m_{0|t}^{x_{t}})\mathrm{d}x_{0})^{1/2}\lesssim\|x_{t}\|(M^{2}-1)\sqrt{d_{x}}$
$\displaystyle\lesssim$ $\displaystyle\|x_{t}\|(M^{2}-1),$

here we used that we can upper bound the operator norm of a positive
semidefinite matrix by its trace from the third to the fourth line. We denoted
by $m_{0|t}^{x_{t}}$ and $\Sigma_{0|t}$ the mean and covariance of
$p^{G}({\mathbf{x}}_{0}^{G}|{\mathbf{x}}_{t}^{G}=x_{t})$. From the second to
last to the last line we used that the $m_{0|t}^{x_{t}}$ depends on $x_{t}$
linearly, and the magnitude of the linear dependence can be bounded uniformly
in $t$ (see equation 39). The integral in the last line is the variance of a
standard normal random variable, which evaluates to $d_{x}$. Furthermore,

$\nabla\Phi_{t}(x_{t})=\mathbb{E}[\nabla\Phi_{0}({\mathbf{x}}_{0})|{\mathbf{x}}_{t}=x_{t}]\leq
L,$

see for example the Proof of Theorem 1 in Pidstrigach et al. (2023).
Therefore,

$\displaystyle\|\nabla_{x_{t}}\mathbb{E}[{\mathbf{x}}_{0}^{G}|{\mathbf{x}}_{t}^{G}=x_{t}]-\nabla_{x_{t}}\mathbb{E}[{\mathbf{x}}_{0}|{\mathbf{x}}_{t}=x_{t}]\|$
$\displaystyle\lesssim$
$\displaystyle\|x_{t}\|(M^{2}-1)+L\|\mathbb{E}[{\mathbf{x}}_{0}|{\mathbf{x}}_{t}=x_{t}]\|$
$\displaystyle\leq$
$\displaystyle\|x_{t}\|(M^{2}-1)+L(\|\mathbb{E}[{\mathbf{x}}_{0}|{\mathbf{x}}_{t}=x_{t}]-\mathbb{E}[{\mathbf{x}}_{0}^{G}|{\mathbf{x}}_{t}^{G}=x_{t}]\|+\|\mathbb{E}[{\mathbf{x}}_{0}^{G}|{\mathbf{x}}_{t}^{G}=x_{t}]\|)$
$\displaystyle\lesssim$
$\displaystyle\|x_{t}\|(M^{2}-1)+L(M^{2}-1+\|x_{t}\|)=\|x_{t}\|(M^{2}-1+L)+(M^{2}-1)L=:B_{x_{t}}$

Furthermore, from equation 28 we see that

$\displaystyle\|\nabla_{x_{t}}\mathbb{E}[{\mathbf{x}}_{0}|{\mathbf{x}}_{t}=x_{t}]\|\leq$
$\displaystyle\|x_{t}\|(M^{2}-1)+L(M^{2}-1+\|x_{t}\|)=B_{x_{t}}$

too. We now subtract the full drifts in equation 20 and equation 21 and arrive
at

$\displaystyle\|\nabla\tilde{f}(x_{t})-f_{t}(x_{t})\|$ $\displaystyle\leq$
$\displaystyle
B_{x_{t}}\times\|T_{1}^{-1}(y-H\mathbb{E}[{\mathbf{x}}^{G}_{0}|{\mathbf{x}}_{t}^{G}=x_{t}])\|$
$\displaystyle+B_{x_{t}}\|T_{1}^{-1}-T_{2}^{-1}\|_{\text{op}}\|y-H\mathbb{E}[{\mathbf{x}}^{G}_{0}|{\mathbf{x}}_{t}^{G}=x_{t}]\|$
$\displaystyle+B_{x_{t}}\|T_{2}^{-1}\|_{\text{op}}\|H\mathbb{E}[{\mathbf{x}}^{G}_{0}|{\mathbf{x}}_{t}^{G}=z_{t}]-H\mathbb{E}[{\mathbf{x}}_{0}|{\mathbf{x}}_{t}=z_{t}]\|$
$\displaystyle\lesssim$ $\displaystyle
B_{x_{t}}+B_{x_{t}}(M^{3}-1)+B_{x_{t}}M(M^{2}-1)$ $\displaystyle\leq$
$\displaystyle B_{x_{t}}(M^{3}-1)=(M^{5}-1)(\|x_{t}\|+L)+(M^{3}-1)L\|x_{t}\|.$

Furthermore, we get that

$\displaystyle\|\nabla\log p_{t}^{G}(x_{t})-\nabla\log
p_{t}(x_{t})\|\leq\|\nabla\log\frac{p_{t}^{G}(x_{t})}{p_{t}(x_{t})}\|\leq\|\nabla_{x_{t}}\Phi_{t}(x_{t})\|\leq
L.$

Plugging in equation 20, we get that

$\|\nabla\log
p_{t}(x_{t})+f_{t}(x_{t})-p_{t}^{G,c}(x_{t})\|\leq(M^{5}-1)(\|x_{t}\|+L)+(M^{3}-1)L\|x_{t}\|+L.$

Using this expression, we see that the drift change is linear in $\|x_{t}\|$.
Since we changed the drift with respect to the reference measure
$\mathbb{P}^{G,c}$, under which $x_{t}$ is Gaussian, we get that

$\mathbb{E}_{\mathbb{P}^{G,c}}\left[\exp\left(\frac{1}{2}\int_{0}^{T}\|\nabla\log
p_{t}(x_{t})+f_{t}(x_{t})-p_{t}^{G,c}(x_{t})\|^{2}\mathrm{d}t\right)\right]<\infty.$

Therefore, by Novikovs condition, our application of Girsanovs theorem was
justified. Furthermore, we see that we the likelihood KL-divergence can be
bounded y

$\displaystyle\mathbb{E}_{\mathbb{P}^{G,c}}\left[\frac{1}{2}\int_{0}^{T}\|\nabla\log
p_{t}(x_{t})+f_{t}(x_{t})-p_{t}^{G,c}(x_{t})\|^{2}\mathrm{d}t\right]$
$\displaystyle\lesssim T((M^{5}-1)(L+1)+(M^{3}-1)L+L)$ $\displaystyle\leq
T((M^{5}-1)(L+1)+L),$

where we used that the mean and covariance of the Gaussian process $x_{t}$
(under $\mathbb{P}^{G,c}$) can be bounded by a constant only depending on $y$,
$\Sigma_{t}$ and $m_{t}$. Taking the square root proves our theorem. ∎

##  Appendix C Variance Exploding SDE (Time rescaled Brownian motion)

The second main type of SDE used in denoising-diffusion is the Variance
Exploding (Song et al., 2020) or time-rescaled Brownian motion, is

$d{\mathbf{x}}_{t}=\sqrt{\frac{dv(t)}{dt}}d{\mathbf{w}}_{t},\quad{\mathbf{x}}_{0}\sim
p_{0}=p_{\text{data}}$ (30)

where $({\mathbf{w}}_{t})_{t\in[0,T]}$ is a Brownian motion and
$v:\mathbb{R}_{\geq 0}\rightarrow\mathbb{R}_{\geq 0}$ is an increasing
function where $v(0)=0$. The transition kernel of the time rescaled Brownian
motion is

$p_{t|0}({\mathbf{x}}_{t}|{\mathbf{x}}_{0})=\mathcal{N}({\mathbf{x}}_{0},v_{t}{\mathbf{I}}_{d_{{\textnormal{x}}}}),$

denoting $v(t)=v_{t}$. The signal ${\mathbf{x}}_{0}$ is sampled from a second
SDE, the reverse process, given in forward time as

$d{\mathbf{x}}_{t}=-\nabla_{{\mathbf{x}}_{t}}\log
p_{t}({\mathbf{x}}_{t})dt-\sqrt{\frac{dv_{t}}{dt}}d\bar{{\mathbf{w}}}_{t},\quad{\mathbf{x}}_{T}\sim
p_{T}.$

##  Appendix D Tweedie’s formula

In this section, we give an alternative derivation of the Tweedie’s formula
for VPSDEs and VESDEs.

Tweedie’s formula (Robbins, 1992; Efron, 2011) gives the minimum mean squared
error (MMSE) estimator of ${\mathbf{x}}_{0}|{\mathbf{x}}_{t}$ when
${\mathbf{x}}_{t}|{\mathbf{x}}_{0}$ is Gaussian and the score
$\nabla_{{\mathbf{x}}_{t}}\log p_{t}({\mathbf{x}}_{t})$ is available.
Tweedie’s formula applied to the time rescaled SDEs is given below.

###  D.1 Variance Preserving SDE (time rescaled Ornstein-Uhlenbeck process)

Using the Variance Preserving SDE or time-rescaled Ornstein-Uhlenbeck process
in equation 2, given the random variable

$\sqrt{\frac{\alpha_{t}}{v_{t}}}{\mathbf{x}}_{0}:={\mathbf{v}}_{0}\sim
p_{{\mathbf{v}}_{0}},$

where, we now define $p_{{\mathbf{v}}_{0}}({\mathbf{v}}_{0})$ as the
probability density of the random variable ${\mathbf{v}}_{0}$ evaluated at the
realization ${\mathbf{v}}_{0}$, then the random variable

$\frac{{\mathbf{x}}_{t}}{\sqrt{v_{t}}}:={\mathbf{v}}_{t}\sim\mathcal{N}({\mathbf{v}}_{0},{\mathbf{I}}_{d_{{\textnormal{x}}}}),$

has a marginal probability density evaluated at the realization
${\mathbf{v}}_{t}$ is the convolution of the random variable
${\mathbf{v}}_{0}$ with the Gaussian kernel

$p_{{\mathbf{v}}_{t}}({\mathbf{v}}_{t})=\int\phi({\mathbf{v}}_{t}-{\mathbf{v}}_{0})p_{{\mathbf{v}}_{0}}({\mathbf{v}}_{0})d{\mathbf{v}}_{0}.$

Then, Tweedie’s formula (Robbins, 1992; Efron, 2011) gives

$\mathbb{E}_{{\mathbf{v}}_{0}\sim
p_{{\mathbf{v}}_{0}|{\mathbf{v}}_{t}}}[{\mathbf{v}}_{0}]={\mathbf{v}}_{t}+\nabla_{{\mathbf{v}}_{t}}\log
p_{{\mathbf{v}}_{t}}({\mathbf{v}}_{t})$

From continuous change of random variables,
$p_{{\mathbf{v}}_{t}}({\mathbf{v}}_{t})=p_{{\mathbf{x}}_{t}}({\mathbf{v}}_{t}\sqrt{v_{t}})\left|\frac{d{\mathbf{x}}_{t}}{d{\mathbf{v}}_{t}}\right|=p_{{\mathbf{x}}_{t}}({\mathbf{x}}_{t})\sqrt{v_{t}},$
giving $\log p_{{\mathbf{v}}_{t}}({\mathbf{v}}_{t})=\log
p_{{\mathbf{x}}_{t}}({\mathbf{x}}_{t})+\log\sqrt{v_{t}},$ which yields

$\nabla_{{\mathbf{v}}_{t}}\log
p_{{\mathbf{v}}_{t}}({\mathbf{v}}_{t})=\sqrt{v_{t}}\nabla_{{\mathbf{x}}_{t}}\log
p_{{\mathbf{x}}_{t}}({\mathbf{x}}_{t}).$

Now, substituting into Tweedie’s formula,

$\mathbb{E}_{{\mathbf{v}}_{0}\sim
p_{{\mathbf{v}}_{0}|{\mathbf{v}}_{t}}}[{\mathbf{v}}_{0}]=\frac{{\mathbf{x}}_{t}}{\sqrt{v_{t}}}+\sqrt{v_{t}}\nabla_{{\mathbf{x}}_{t}}\log
p_{{\mathbf{x}}_{t}}({\mathbf{x}}_{t}),$

giving

$\mathbb{E}_{{\mathbf{x}}_{0}\sim
p_{{\mathbf{x}}_{0}|{\mathbf{x}}_{t}}}[{\mathbf{x}}_{0}]=\frac{1}{\sqrt{\alpha_{t}}}{\mathbf{x}}_{t}+\frac{1}{\sqrt{\alpha_{t}}}v_{t}\nabla_{{\mathbf{x}}_{t}}\log
p_{{\mathbf{x}}_{t}}({\mathbf{x}}_{t}).$

The covariance is given as (Robbins, 1992; Efron, 2011),
$\mathbb{V}_{{\mathbf{v}}_{0}\sim
p_{{\mathbf{v}}_{0}|{\mathbf{v}}_{t}}}[{\mathbf{v}}_{0}|{\mathbf{v}}_{t}]=I+\nabla_{{\mathbf{v}}_{t}}\nabla_{{\mathbf{v}}_{t}}\log
p_{{\mathbf{v}}_{t}}({\mathbf{v}}_{t}),$ where

$\displaystyle\nabla_{{\mathbf{v}}_{t}}\nabla_{{\mathbf{v}}_{t}}\log
p_{{\mathbf{v}}_{t}}({\mathbf{v}}_{t})$
$\displaystyle=\nabla_{{\mathbf{v}}_{t}}(\sqrt{v_{t}}\nabla_{{\mathbf{x}}_{t}}\log
p_{{\mathbf{x}}_{t}}({\mathbf{x}}_{t}))$ (31)
$\displaystyle=v_{t}\nabla^{2}_{{\mathbf{x}}_{t}}\log
p_{{\mathbf{x}}_{t}}({\mathbf{x}}_{t})$ (32)

and so

$\mathbb{V}_{{\mathbf{v}}_{0}\sim
p_{{\mathbf{v}}_{0}|{\mathbf{v}}_{t}}}[{\mathbf{v}}_{0}]={\mathbf{I}}_{d_{{\textnormal{x}}}}+v_{t}\nabla^{2}_{{\mathbf{x}}_{t}}\log
p_{{\mathbf{x}}_{t}}({\mathbf{x}}_{t})$

giving

$\mathbb{V}_{{\mathbf{x}}_{0}\sim
p_{{\mathbf{x}}_{0}|{\mathbf{x}}_{t}}}[{\mathbf{x}}_{0}]=\frac{v_{t}}{\alpha_{t}}({\mathbf{I}}_{d_{{\textnormal{x}}}}+v_{t}\nabla_{{\mathbf{x}}_{t}}\log
p_{{\mathbf{x}}_{t}}({\mathbf{x}}_{t})).$

In practice, we make the approximation $\nabla_{{\mathbf{x}}_{t}}\log
p_{{\mathbf{x}}_{t}}({\mathbf{x}}_{t})\approx\mathsf{s}_{\theta}({\mathbf{x}}_{t},t)$,
giving

$\mathbb{E}_{{\mathbf{x}}_{0}\sim
p_{{\mathbf{x}}_{0}|{\mathbf{x}}_{t}}}[{\mathbf{x}}_{0}]\approx{\mathbf{m}}_{0|t}:=\frac{1}{\sqrt{\alpha_{t}}}{\mathbf{x}}_{t}+\frac{1}{\sqrt{\alpha_{t}}}v_{t}\mathsf{s}_{\theta}({\mathbf{x}}_{t},t)$

and

$\mathbb{V}_{{\mathbf{x}}_{0}\sim
p_{{\mathbf{x}}_{0}|{\mathbf{x}}_{t}}}[{\mathbf{x}}_{0}]\approx\frac{v_{t}}{\alpha_{t}}({\mathbf{I}}_{d_{{\textnormal{x}}}}+v_{t}\nabla_{{\mathbf{x}}_{t}}\mathsf{s}_{\theta}({\mathbf{x}}_{t},t)),$

which yields a Gaussian approximation that matches the first and second
moments,

$\displaystyle
p_{{\mathbf{x}}_{0}|{\mathbf{x}}_{t}}({\mathbf{x}}_{0}|{\mathbf{x}}_{t})$
$\displaystyle\approx\mathcal{N}\left({\mathbf{m}}_{0|t},\frac{v_{t}}{\alpha_{t}}({\mathbf{I}}_{d_{{\textnormal{x}}}}+v_{t}\nabla_{{\mathbf{x}}_{t}}\mathsf{s}_{\theta}({\mathbf{x}}_{t},t))\right)$
(34)
$\displaystyle=\mathcal{N}\left({\mathbf{m}}_{0|t},\frac{v_{t}}{\sqrt{\alpha_{t}}}\nabla_{{\mathbf{x}}_{t}}{\mathbf{m}}_{0|t}\right).$
(35)

The matrix
$\frac{v_{t}}{\alpha_{t}}({\mathbf{I}}_{d_{{\textnormal{x}}}}+v_{t}\nabla_{{\mathbf{x}}_{t}}\mathsf{s}_{\theta}({\mathbf{x}}_{t},t))$
needs to be inverted to calculate the log likelihood, and therefore must be
both symmetric and positive definite for all time, which puts the requirement
that $\mathsf{s}_{\theta}({\mathbf{x}}_{t},t)$ can be written as the negative
gradient of a potential, however it is only approximated as such
$\mathsf{s}_{\theta}({\mathbf{x}}_{t},t)\approx\nabla\log
p_{{\mathbf{x}}_{t}}({\mathbf{x}}_{t})$.

###  D.2 Variance Exploding SDE (time rescaled Brownian motion)

Using the Variance Exploding SDE (Song et al., 2020) or time rescaled Brownian
motion, given the random variable

$\sqrt{\frac{1}{v_{t}}}{\mathbf{x}}_{0}:={\mathbf{v}}_{0}\sim
p_{{\mathbf{z}}_{0}},$ (36)

a similar calculation as in D.1 gives

$\mathbb{E}_{{\mathbf{x}}_{0}\sim
p_{{\mathbf{x}}_{0}|{\mathbf{x}}_{t}}}[{\mathbf{x}}_{0}]={\mathbf{x}}_{t}+v_{t}\nabla_{{\mathbf{x}}_{t}}\log
p_{{\mathbf{x}}_{t}}({\mathbf{x}}_{t}),$

and

$\mathbb{V}_{{\mathbf{x}}_{0}\sim
p_{{\mathbf{x}}_{0}|{\mathbf{x}}_{t}}}[{\mathbf{x}}_{0}]=v_{t}({\mathbf{I}}_{d_{{\textnormal{x}}}}+v_{t}\nabla_{{\mathbf{x}}_{t}}\log
p_{{\mathbf{x}}_{t}}({\mathbf{x}}_{t})).$

Again, in practice, we make the approximation $\nabla_{{\mathbf{x}}_{t}}\log
p_{{\mathbf{x}}_{t}}({\mathbf{x}}_{t})=\mathsf{s}_{\theta}({\mathbf{x}}_{t},t)$,
giving

$\mathbb{E}_{{\mathbf{x}}_{0}\sim
p_{{\mathbf{x}}_{0}|{\mathbf{x}}_{t}}}[{\mathbf{x}}_{0}]\approx{\mathbf{m}}_{0|t}:={\mathbf{x}}_{t}+v_{t}\mathsf{s}_{\theta}({\mathbf{x}}_{t},t)$

and

$\mathbb{V}_{{\mathbf{x}}_{0}\sim
p_{{\mathbf{x}}_{0}|{\mathbf{x}}_{t}}}[{\mathbf{x}}_{0}]\approx
v_{t}({\mathbf{I}}_{d_{{\textnormal{x}}}}+v_{t}\nabla_{{\mathbf{x}}_{t}}\mathsf{s}_{\theta}({\mathbf{x}}_{t},t)),$

which yields a Gaussian approximation that matches the first and second
moments,

$\displaystyle
p_{{\mathbf{x}}_{0}|{\mathbf{x}}_{t}}({\mathbf{x}}_{0}|{\mathbf{x}}_{t})$
$\displaystyle\approx\mathcal{N}\left({\mathbf{m}}_{0|t},v_{t}({\mathbf{I}}_{d_{{\textnormal{x}}}}+v_{t}\nabla_{{\mathbf{x}}_{t}}\mathsf{s}_{\theta}({\mathbf{x}}_{t},t))\right)$
(37)
$\displaystyle=\mathcal{N}\left({\mathbf{m}}_{0|t},v_{t}\nabla_{{\mathbf{x}}_{t}}{\mathbf{m}}_{0|t}\right).$
(38)

##  Appendix E Algorithmic details and numerics

The code for all of the experiments and instructions to run them will be made
available after publication.

Algorithm 1 TMPD-D (Ancestral sampling, VP)

${\mathbf{x}}_{N}\sim\mathcal{N}(0,{\mathbf{I}}_{d_{{\mathbf{x}}}})$

for $n={N-1,\ldots,0}$ do

${\mathbf{m}}_{0|t}=\frac{1}{\sqrt{\alpha_{n}}}({\mathbf{x}}_{n}+v_{n}\mathsf{s}_{\theta}(\hat{{\mathbf{x}}}_{n},t_{n}))$

${\mathbf{m}}_{0|t}^{{\mathbf{y}}}={\mathbf{m}}_{0|t}+\frac{{v_{n}}}{\sqrt{\alpha_{n}}}\nabla_{{\mathbf{x}}_{n}}{\mathbf{m}}_{0|t}{\mathbf{H}}^{\top}({\mathbf{H}}\frac{{v_{n}}}{\sqrt{\alpha_{n}}}\nabla_{{\mathbf{x}}_{n}}{\mathbf{m}}_{0|t}{\mathbf{H}}^{\top}+\sigma_{y}^{2}{\mathbf{I}}_{d_{{\mathbf{y}}}})^{-1}({\mathbf{y}}-{\mathbf{H}}{\mathbf{m}}_{0|t})$

${\mathbf{z}}_{n}\sim\mathcal{N}(0,{\mathbf{I}}_{d_{{\mathbf{x}}}})$

$\sigma_{n}=\sqrt{(1-\alpha_{n-1})\beta_{n}/(1-\alpha_{n})}$

${\mathbf{x}}_{n-1}\leftarrow\frac{\sqrt{1-\beta_{n}}(1-\alpha_{n-1})}{1-\alpha_{n}}{\mathbf{x}}_{n}+\frac{\sqrt{\alpha_{n-1}}\beta_{n}}{1-\alpha_{n}}{\mathbf{x}}^{\prime}_{0}+\sigma_{n}{\mathbf{z}}_{n}$

end for

return $\hat{{\mathbf{x}}}_{0}$

###  E.1 Computationally cheaper approximation of Gaussian Projection

We discuss the computational time and memory cost of TMPD compared to $\Pi$GDM
methods here. Let $d_{{\textnormal{y}}}$ is the dimensions of the observation,
$d_{{\textnormal{x}}}$ are the dimension of the image, $T_{s}$ is the time
complexity of evaluating the score network. Computing a jacobian-vector
product has time complexity $T_{s}$. Then, the time complexity of TMPD (not
including fast matrix-vector products) is
$\mathcal{O}(N(d_{{\textnormal{y}}}^{3}+T_{s}d_{{\textnormal{y}}}+T_{s}+T_{s}))$.
$\Pi$GDM comes at a smaller computational cost due to needing only
$\mathcal{O}(1)$ vector-jacobian-products, instead of
$\mathcal{O}(d_{{\mathbf{y}}})$, resulting in a time complexity (bar the
matrix-vector products) of $\mathcal{O}(N(T_{s}+T_{s}))$. DPS comes at the
same complexity as $\Pi$GDM.

Whereas TMPD requires calculating the Jacobian which has memory complexity of
atleast $\mathcal{O}(d_{{\textnormal{x}}}d_{{\textnormal{y}}})$. This is too
large for high resolution image problems where the dimension of the
observation is large. In comparison, the memory complexity of $\Pi$GDM depends
on the observation operator, but for the class of problems that are explored
in Song et al. (2023) is $\mathcal{O}(d_{{\textnormal{x}}})$. DPS comes at the
same memory complexity as $\Pi$GDM.

To make the computational cost of TMPD smaller, we can make an approximation
of the Gaussian Projection that requires fewer jacobian-vector products and
does not require linear solves. One approximation that we found useful for
sampling from high dimensional diffusion models, e.g., high resolution images,
is denoted here as diagonal Gaussian Projection (DTMPD). Instead of the full
second moment, DTMPD uses the diagonal of the second moment. Intuitively, this
approximation will perform well empirically since it is a similar
approximation to $\Pi$GDM that assumes dimensional independence of the
distribution $p({\mathbf{x}}_{0}|{\mathbf{x}}_{t})$,

$\text{DTMPD}\quad\nabla_{{\mathbf{x}}_{t}}\log
p_{{\mathbf{y}}|t}({\mathbf{y}}|{\mathbf{x}}_{t})\approx\nabla_{{\mathbf{x}}_{t}}{\mathbf{m}}_{0|t}{\mathbf{H}}^{\top}({\mathbf{H}}\frac{v_{t}}{\sqrt{\alpha_{t}}}\text{diag}(\nabla_{{\mathbf{x}}_{t}}{\mathbf{m}}_{0|t}){\mathbf{H}}^{\top}+\sigma^{2}_{{\textnormal{y}}}{\mathbf{I}}_{d_{{\textnormal{y}}}})^{-1}({\mathbf{y}}-{\mathbf{H}}{\mathbf{m}}_{0|t}).$

Whilst this approximation does not require a linear solve, taking out the
$\mathcal{O}(Nd_{{\mathbf{y}}}^{3})$ time complexity term, we would still like
to take out the $\mathcal{O}(NT_{s}d_{{\textnormal{y}}})$ time complexity and
$\mathcal{O}(d_{{\textnormal{x}}}d_{{\textnormal{y}}})$ memory complexity term
from calculating and storing the Jacobian, respectively, since this is too
large for solving high resolution image applications. A further approximation
approximates the diagonal of the Jacobian by the row sum of the Jacobian which
only requires one vector jacobian product and brings the memory and time
complexity of DTMPD down to that of $\Pi$GDM. We use this approximation in the
image experiments.

###  E.2 Gaussian

When the data distribution $p_{0}({\mathbf{x}}_{0})$ is a (multivariate)
Gaussian, then the reverse SDE is a linear SDE and we can calculate all of the
terms needed to sample from the target posterior using diffusion explicitly.
Moreover, we can sample from the target posterior using a direct or implicit
method such as Cholesky decomposition. In this simple example, we compare
samples from the direct method to various conditional diffusion methods (TMPD,
$\Pi$GDM, DPM), by plotting a sample estimate of the $L^{2}$ Wasserstein
distance between the sample and the target Gaussian measures (Givens & Shortt,
1984) by using the analytical mean and covariance of the target distribution
and empirical estimate of the mean and covariance of the sample distribution.
To generate $p_{0}({\mathbf{x}}_{0})$, we use an equally spaced grid of
vectors ${\mathbf{u}}_{i}\in[-5.0,5.0]^{2}$ for $i\in\\{1,2,...,32\\}^{2}$,
pick a Matern 5/2 kernel for the covariance function
$k({\mathbf{u}}_{i},{\mathbf{u}}_{j})=\left(1+\sqrt{5}|{\mathbf{u}}_{i}-{\mathbf{u}}_{j}|+\frac{5}{3}|{\mathbf{u}}_{i}-{\mathbf{u}}_{j}|^{2}\right)\exp\left(-\sqrt{3}|{\mathbf{u}}_{i}-{\mathbf{u}}_{j}|\right)$
which defines the prior
$p_{0}({\mathbf{x}}_{0})=\mathcal{N}({\mathbf{m}}_{0},{\mathbf{C}}_{0})$
covariance ${\mathbf{C}}_{0}\in\mathbb{R}^{1024\times 1024}$ where
${{\mathbf{C}}_{0}}_{ij}=k({\mathbf{u}}_{i},{\mathbf{u}}_{j})$ and we define
the mean as a zero vector ${\mathbf{m}}_{0}=\mathbf{0}\in\mathbb{R}^{1024}$.
To compute analytically the distribution of
$p_{0|{\mathbf{y}}}({\mathbf{x}}_{0}|{\mathbf{y}})$, we sample once
${\mathbf{y}}={\mathbf{H}}{\mathbf{x}}_{0}+{\mathbf{z}},{\mathbf{z}}\sim\mathcal{N}(\mathbf{0},\sigma_{{\textnormal{y}}}^{2}{\mathbf{I}}_{d_{{\textnormal{y}}}}),$
where ${\mathbf{H}}$ is an inpainting type matrix that we have chosen to
subsample $d_{y}<d_{x}$ elements of ${\mathbf{x}}_{0}$. We use the standard
Gaussian formula to calculate the mean and covariance of
${\mathbf{x}}_{0}|{\mathbf{y}}$ which in this case are a complete description
of $p_{0|{\mathbf{y}}}({\mathbf{x}}_{0}|{\mathbf{y}})$. The $L^{2}$
Wasserstein estimate is plotted over an increasing sample size $N\in[9,1500]$
and for $\sigma_{{\textnormal{y}}}=0.1$ in Figure 1.

We provide an illustration of the mean and uncertainty captured by the
diffusion samples using $1500$ samples of each diffusion model to produce a
Monte-Carlo estimate of the mean and diagonal variance vector, and compare
these to the exact mean and diagonal variance.

Below, we provide a calculation comparing the exact diffusion posterior for
the linear diffusion posterior sde to the approximations used in $\Pi$GDM
(Song et al., 2023) and TMPD. Let us assume that the target distribution be
known
$p_{0}({\mathbf{x}}_{0})=\mathcal{N}({\mathbf{m}}_{0},{\mathbf{C}}_{0})$.
Then, via Bayes’ rule,

$\displaystyle\log p_{0|t}({\mathbf{x}}_{0}|{\mathbf{x}}_{t})=$
$\displaystyle\log p_{t|0}({\mathbf{x}}_{t}|{\mathbf{x}}_{0})+\log
p_{0}({\mathbf{x}}_{0})+{\text{constant}}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}({\mathbf{x}}_{t}-\sqrt{\alpha_{t}}{\mathbf{x}}_{0})^{\top}(v_{t}{\mathbf{I}}_{d_{{\textnormal{x}}}})^{-1}({\mathbf{x}}_{t}-\sqrt{\alpha_{t}}{\mathbf{x}}_{0})$
$\displaystyle-\frac{1}{2}({\mathbf{x}}_{0}-{\mathbf{m}}_{0})^{\top}{\mathbf{C}}_{0}^{-1}({\mathbf{x}}_{0}-{\mathbf{m}}_{0})+\text{constant}$
$\displaystyle=$
$\displaystyle-\frac{1}{2}(\sqrt{\alpha_{t}}{\mathbf{x}}_{0})^{\top}(v_{t}{\mathbf{I}}_{d_{{\textnormal{x}}}})^{-1}(\sqrt{\alpha_{t}}{\mathbf{x}}_{0})+{\mathbf{x}}_{t}^{\top}(v_{t}{\mathbf{I}}_{d_{{\textnormal{x}}}})^{-1}(\sqrt{\alpha_{t}}{\mathbf{x}}_{0})$
$\displaystyle-\frac{1}{2}{\mathbf{x}}_{0}^{\top}{\mathbf{C}}_{0}^{-1}{\mathbf{x}}_{0}+{\mathbf{m}}_{0}^{\top}{\mathbf{C}}_{0}^{-1}{\mathbf{x}}_{0}+\text{constant}$
$\displaystyle=$
$\displaystyle-\frac{1}{2}({\mathbf{x}}_{0}-{\mathbf{m}}_{t})^{\top}{\bm{\Sigma}}_{t}^{-1}({\mathbf{x}}_{0}-{\mathbf{m}}_{t}),$

so

$p_{0|t}({\mathbf{x}}_{0}|{\mathbf{x}}_{t})=\mathcal{N}({\mathbf{m}}_{t},{\bm{\Sigma}}_{t})$
(39)

where

${\bm{\Sigma}}_{t}=\left((\frac{v_{t}}{\alpha_{t}}{\mathbf{I}}_{d_{{\textnormal{x}}}})^{-1}+{\mathbf{C}}_{0}^{-1}\right)^{-1}$
(40)

and

${\mathbf{m}}_{t}={\bm{\Sigma}}_{t}\left(\frac{\sqrt{\alpha_{t}}}{v_{t}}{\mathbf{x}}_{t}+{\mathbf{C}}_{0}^{-1}{\mathbf{m}}_{0}\right).$

We also have that,

$p_{t}({\mathbf{x}}_{t})=\mathcal{N}(\sqrt{\alpha_{t}}{\mathbf{m}}_{0},{\mathbf{C}}_{t})$

where
${\mathbf{C}}_{t}=\alpha_{t}{\mathbf{C}}_{0}+v_{t}{\mathbf{I}}_{d_{{\textnormal{x}}}},$
and thus

$\log
p_{t}({\mathbf{x}}_{t})=-\frac{1}{2}({\mathbf{x}}_{t}-{\mathbf{m}}_{0}\sqrt{\alpha_{t}})^{\top}{\mathbf{C}}_{t}^{-1}({\mathbf{x}}_{t}-{\mathbf{m}}_{0}\sqrt{\alpha_{t}}).$

From

$p_{{\mathbf{y}}|t}({\mathbf{y}}|{\mathbf{x}}_{t})=\int
p_{{\mathbf{y}}|0}({\mathbf{y}}|{\mathbf{x}}_{0})p_{0|t}({\mathbf{x}}_{0}|{\mathbf{x}}_{t})d{\mathbf{x}}_{0}=\mathcal{N}({\mathbf{H}}{\mathbf{m}}_{t},{\mathbf{H}}{\bm{\Sigma}}_{t}{\mathbf{H}}^{\top}+\sigma_{{\textnormal{y}}}^{2}{\mathbf{I}}_{d_{{\textnormal{y}}}}),$

we have that

$\log
p_{t}({\mathbf{y}}|{\mathbf{x}}_{t})=-\frac{1}{2}({\mathbf{y}}-{\mathbf{H}}{\mathbf{m}}_{t})^{\top}({\mathbf{H}}{\bm{\Sigma}}_{t}{\mathbf{H}}^{T}+\sigma_{{\textnormal{y}}}^{2}{\mathbf{I}}_{d_{{\textnormal{y}}}})^{-1}({\mathbf{y}}-{\mathbf{H}}{\mathbf{m}}_{t}).$

The posterior score can be calculated (see Cardoso et al. (2023)) from
applying Bayes’ theorem,

$\log p_{t|{\mathbf{y}}}({\mathbf{x}}_{t}|{\mathbf{y}})=\log
p_{t}({\mathbf{x}}_{t})+\log
p_{{\mathbf{y}}|t}({\mathbf{y}}|{\mathbf{x}}_{t})+\text{Constant},$

and taking gradients with respect to the state ${\mathbf{x}}_{t}$ gives

$\displaystyle\nabla_{{\mathbf{x}}_{t}}\log
p_{t|y}({\mathbf{x}}_{t}|{\mathbf{y}})=$
$\displaystyle\nabla_{{\mathbf{x}}_{t}}\log
p_{t}({\mathbf{x}}_{t})+\nabla_{{\mathbf{x}}_{t}}\log
p_{{\mathbf{y}}|t}({\mathbf{y}}|{\mathbf{x}}_{t})$ (41) $\displaystyle=$
$\displaystyle-(\alpha_{t}{\mathbf{C}}_{0}+v_{t}{\mathbf{I}}_{d_{{\textnormal{x}}}})^{-1}({\mathbf{x}}_{t}-{\mathbf{m}}_{0}\sqrt{\alpha_{t}})$
(42)
$\displaystyle+\frac{\sqrt{\alpha_{t}}}{v_{t}}{\bm{\Sigma}}_{t}^{\top}{\mathbf{H}}^{\top}({\mathbf{H}}{\bm{\Sigma}}_{t}{\mathbf{H}}^{\top}+\sigma^{2}_{{\textnormal{y}}}{\mathbf{I}}_{d_{{\textnormal{y}}}})^{-1}({\mathbf{y}}-{\mathbf{H}}{\bm{\Sigma}}_{t}(\frac{\sqrt{\alpha_{t}}}{v_{t}}{\mathbf{x}}_{t}+{\mathbf{C}}_{0}^{-1}{\mathbf{m}}_{0})).$
(43)

Setting ${\mathbf{m}}_{0}=\mathbf{0}$ in (43) gives,

$\displaystyle\nabla_{{\mathbf{x}}_{t}}\log
p_{t|{\mathbf{y}}}({\mathbf{x}}_{t}|{\mathbf{y}})=$
$\displaystyle-(\alpha_{t}{\mathbf{C}}_{0}+v_{t}{\mathbf{I}}_{d_{{\textnormal{x}}}})^{-1}{\mathbf{x}}_{t}$
(44)
$\displaystyle+\frac{\sqrt{\alpha_{t}}}{v_{t}}{\bm{\Sigma}}_{t}^{\top}{\mathbf{H}}^{\top}({\mathbf{H}}{\bm{\Sigma}}_{t}{\mathbf{H}}^{\top}+\sigma_{{\textnormal{y}}}^{2}{\mathbf{I}}_{d_{{\textnormal{y}}}})^{-1}({\mathbf{y}}-\frac{\sqrt{\alpha_{t}}}{v_{t}}{\mathbf{H}}{\bm{\Sigma}}_{t}{\mathbf{x}}_{t})$
(45)

####  E.2.1 Comparison of the exact posterior to the approximation made in
Pseudo-Inverse-Guidance

We aim to compare (45), to the approximation used in $\Pi$GDM (Song et al.,
2023);

$p_{0|t}({\mathbf{x}}_{0}|{\mathbf{x}}_{t})\approx\mathcal{N}({\mathbf{m}}_{0|t},r_{t}^{2}{\mathbf{I}}_{d_{{\textnormal{x}}}}),$

where ${\mathbf{m}}_{0|t}={\mathbf{x}}_{t}+v_{t}\nabla_{{\mathbf{x}}_{t}}\log
p_{t}({\mathbf{x}}_{t})$, which is the MMSE estimate of
${\mathbf{x}}_{0}|{\mathbf{x}}_{t}$ and $r_{t}$ is chosen empirically. This
gives
$p_{{\mathbf{y}}|t}({\mathbf{y}}|{\mathbf{x}}_{t})\approx\mathcal{N}({\mathbf{H}}{\mathbf{m}}_{0|t},r_{t}^{2}{\mathbf{H}}{\mathbf{H}}^{T}+\sigma_{{\textnormal{y}}}^{2}{\mathbf{I}}_{d_{{\textnormal{y}}}})$,
which in turn gives

$\displaystyle\nabla_{{\mathbf{x}}_{t}}\log
p_{t|y}({\mathbf{x}}_{t}|{\mathbf{y}})\approx$
$\displaystyle\nabla_{{\mathbf{x}}_{t}}\log p_{t}({\mathbf{x}}_{t})$ (46)
$\displaystyle+\nabla_{{\mathbf{x}}_{t}}{\mathbf{m}}_{0|t}{\mathbf{H}}^{\top}(r_{t}^{2}{\mathbf{H}}{\mathbf{H}}^{\top}+\sigma_{{\textnormal{y}}}^{2}{\mathbf{I}}_{d_{{\textnormal{y}}}})^{-1}({\mathbf{y}}-{\mathbf{H}}{\mathbf{m}}_{0|t})$
(47)

which is computationally tractable in the standard diffusion model setting
where the score is non-linear and approximated as
$\nabla_{{\mathbf{x}}_{t}}\log
p_{t}({\mathbf{x}}_{t})\approx\mathsf{s}_{\theta}({\mathbf{x}}_{t},t)$.
Substituting in the known score and again setting
${\mathbf{m}}_{0}=\mathbf{0}$ score to compare to the linear case (45) gives,

$\displaystyle\nabla_{{\mathbf{x}}_{t}}\log
p_{t|{\mathbf{y}}}({\mathbf{x}}_{t}|{\mathbf{y}})\approx$
$\displaystyle\nabla_{{\mathbf{x}}_{t}}\log p_{t}({\mathbf{x}}_{t})$ (48)
$\displaystyle+\nabla_{{\mathbf{x}}_{t}}{\mathbf{m}}_{0|t}{\mathbf{H}}^{\top}(r_{t}^{2}{\mathbf{H}}{\mathbf{H}}^{\top}+\sigma_{{\textnormal{y}}}^{2}I_{d_{{\textnormal{y}}}})^{-1}(y-{\mathbf{H}}{\mathbf{m}}_{0|t})$
(49) $\displaystyle=$
$\displaystyle-(\alpha_{t}{\mathbf{C}}_{0}+v_{t}{\mathbf{I}}_{d_{{\textnormal{x}}}})^{-1}{\mathbf{x}}_{t}$
(50)
$\displaystyle+\frac{1}{\sqrt{\alpha_{t}}}({\mathbf{I}}_{d_{{\textnormal{x}}}}-v_{t}(\alpha_{t}{\mathbf{C}}_{0}+v_{t}{\mathbf{I}}_{d_{{\textnormal{x}}}})^{-1}){\mathbf{H}}^{\top}(r_{t}^{2}{\mathbf{H}}{\mathbf{H}}^{T}+\sigma_{{\textnormal{y}}}^{2}{\mathbf{I}}_{d_{{\textnormal{y}}}})^{-1}$
(51)
$\displaystyle\quad\quad({\mathbf{y}}-\frac{1}{\sqrt{\alpha_{t}}}{\mathbf{H}}({\mathbf{I}}_{d_{{\textnormal{x}}}}-v_{t}(\alpha_{t}{\mathbf{C}}_{0}+v_{t}{\mathbf{I}}_{d_{{\textnormal{x}}}})^{-1}){\mathbf{x}}_{t})$
(52)

comparing therms, $\Pi$GDM is making the approximation

$\nabla_{{\mathbf{x}}_{t}}{\mathbf{m}}_{0|t}=\frac{1}{\sqrt{\alpha_{t}}}{\mathbf{I}}_{d_{{\textnormal{x}}}}-\frac{v_{t}}{\sqrt{\alpha_{t}}}(\alpha_{t}{\mathbf{C}}_{0}+v_{t}{\mathbf{I}}_{d_{{\textnormal{x}}}})^{-1}\approx\frac{\sqrt{\alpha_{t}}}{v_{t}}{\bm{\Sigma}}_{t}.$

Note that since,

$\displaystyle\frac{\sqrt{\alpha_{t}}}{v_{t}}{\bm{\Sigma}}_{t}$
$\displaystyle=\frac{\sqrt{\alpha_{t}}}{v_{t}}({\mathbf{C}}_{0}^{-1}+(\frac{v_{t}}{\alpha_{t}}I_{d_{{\textnormal{x}}}})^{-1})^{-1}$
(53)
$\displaystyle=\frac{1}{v_{t}\sqrt{\alpha_{t}}}((\alpha_{t}{\mathbf{C}}_{0})^{-1}+(v_{t}{\mathbf{I}}_{d_{{\textnormal{x}}}})^{-1})^{-1}$
(54)
$\displaystyle=\frac{1}{v_{t}\sqrt{\alpha_{t}}}(v_{t}{\mathbf{I}}-v_{t}^{2}(\alpha_{t}{\mathbf{C}}_{0}+v_{t}{\mathbf{I}}_{d_{{\textnormal{x}}}})^{-1})\quad\text{Woodbury
identity}$ (55)
$\displaystyle=\frac{1}{\sqrt{\alpha_{t}}}({\mathbf{I}}-v_{t}(\alpha_{t}{\mathbf{C}}_{0}+v_{t}{\mathbf{I}}_{d_{{\textnormal{x}}}})^{-1})$
(56) $\displaystyle=\nabla_{{\mathbf{x}}_{t}}{\mathbf{m}}_{0|t},$ (57)

this is exact. The other approximation $\Pi$GDM is making is
$r_{t}^{2}{\mathbf{I}}\approx{\bm{\Sigma}}_{t}$, and note that this
approximation is accurate with $r^{2}_{t}=v_{t}$ as $t\rightarrow 0$ but
inaccurate as $t\rightarrow 1$. But the exact linear SDE can be recovered by
instead using TMPD,

$\displaystyle\nabla_{{\mathbf{x}}_{t}}\log
p_{t|y}({\mathbf{x}}_{t}|{\mathbf{y}})=$
$\displaystyle\nabla_{{\mathbf{x}}_{t}}\log p_{t}({\mathbf{x}}_{t})$ (58)
$\displaystyle+\nabla_{{\mathbf{x}}_{t}}{\mathbf{m}}_{0|t}{\mathbf{H}}^{\top}({\mathbf{H}}\frac{v_{t}}{\sqrt{\alpha_{t}}}\nabla_{{\mathbf{x}}_{t}}{\mathbf{m}}_{0|t}{\mathbf{H}}^{\top}+\sigma_{{\textnormal{y}}}^{2}I_{d_{{\textnormal{y}}}})^{-1}({\mathbf{y}}-{\mathbf{H}}{\mathbf{m}}_{0|t})$
(59) $\displaystyle=$ $\displaystyle\nabla_{{\mathbf{x}}_{t}}\log
p_{t}({\mathbf{x}}_{t})$ (60)
$\displaystyle+\frac{\sqrt{\alpha_{t}}}{v_{t}}{\bm{\Sigma}}_{t}{\mathbf{H}}^{\top}({\mathbf{H}}{\bm{\Sigma}}_{t}{\mathbf{H}}^{\top}+\sigma_{{\textnormal{y}}}^{2}I_{d_{{\textnormal{y}}}})^{-1}({\mathbf{y}}-\frac{\sqrt{\alpha_{t}}}{v_{t}}{\mathbf{H}}{\bm{\Sigma}}_{t}{\mathbf{x}}_{t})$
(61)

###  E.3 GMM

For a given dimension $d_{{\textnormal{x}}}$, we consider $p_{0}$ a mixture of
25 Gaussian random variables. The components have mean
$\mu_{i,j}:=(8i,8j,...,8i,8j)\in\mathbb{R}^{d_{{\textnormal{x}}}}\text{for}(i,j)\in\\{-2,-1,0,1,2\\}^{2}$
and unit variance. We have set the associated unnormalized weights
$\omega_{i,j}=1.0$. We have set $\sigma_{\delta}^{2}=10^{-4}$.

Note that $p_{t}({\mathbf{x}}_{t})=\int
p_{t|0}({\mathbf{x}}_{t}|{\mathbf{x}}_{0})p_{0}({\mathbf{x}}_{0})d{\mathbf{x}}_{0}$.
As $p_{0}({\mathbf{x}}_{0})$ is a mixture of Gaussians,
$p_{t}({\mathbf{x}}_{t})$ is also a mixture of Gaussians with means
$\sqrt{\alpha_{t}}\mu_{i,j}$ and unitary variances. Therefore, using automatic
differentiation libraries, we can calculate $\nabla_{{\mathbf{x}}_{t}}\log
p_{t}({\mathbf{x}}_{t})$. We chose $\beta_{\text{max}}=500.0$ and
$\beta_{\text{min}}=0.1$. We use 1000 timesteps for the time-discretization.
For the pair of dimensions and chosen observation noise standard deviation
$(d_{{\textnormal{x}}},d_{{\textnormal{y}}},\sigma_{{\textnormal{y}}})$ the
measurement model $(y,{\mathbf{H}})$ is drawn as follows:

  * •

${\mathbf{H}}$: We first draw
$\tilde{{\mathbf{H}}}\sim\mathcal{N}(\mathbf{0}_{d_{{\textnormal{y}}}\times
d_{{\textnormal{x}}}},{\mathbf{I}}_{d_{{\textnormal{y}}}\times
d_{{\textnormal{x}}}})$ and compute the SVD decomposition of
$\tilde{{\mathbf{H}}}={\bm{U}}{\bm{S}}{\bm{V}}^{\top}$. Then, we sample for
$(i,j)\in\\{-2,-1,0,1,2\\}^{2}$, $s_{i,j}$ according to a uniform in $[0,1]$.
Finally, we set
${\mathbf{H}}={\bm{U}}\text{diag}(\\{s_{i,j}\\}_{(i,j)\in\\{-2,-1,0,1,2\\}^{2}}){\bm{V}}^{\top}$.

  * •

${\mathbf{y}}$: We then draw ${\mathbf{x}}_{*}\sim p_{0}$ and set
${\mathbf{y}}:={\mathbf{H}}x_{*}+{\mathbf{z}}$ where
${\mathbf{z}}\sim\mathcal{N}(\mathbf{0},\sigma_{{\textnormal{y}}}^{2}{\mathbf{I}}_{d_{{\textnormal{y}}}})$.

Once we have drawn both ${\mathbf{x}}_{*}\sim p_{0}$ and
$({\mathbf{y}},{\mathbf{H}},\sigma_{{\textnormal{y}}})$, the posterior can be
exactly calculated using Bayes formula and gives a mixture of Gaussians with
mixture components $c_{i,j}$ and associated weights $\tilde{\omega}_{i,j}$,

$\displaystyle c_{i,j}$
$\displaystyle:=\mathcal{N}({\bm{\Sigma}}({\mathbf{H}}^{\top}{\mathbf{y}}/\sigma_{{\textnormal{y}}}^{2}+\mu_{i,j}),{\bm{\Sigma}}),$
(62) $\displaystyle\tilde{\omega}_{i}$
$\displaystyle:=\omega_{i}\mathcal{N}({\mathbf{y}};{\mathbf{H}}\mu_{i,j},\sigma_{\delta}^{2}{\mathbf{I}}_{d_{{\textnormal{x}}}}+{\mathbf{H}}{\mathbf{H}}^{\top}),$
(63)

where
${\bm{\Sigma}}=({\mathbf{I}}_{d_{{\textnormal{x}}}}+\frac{1}{\sigma_{\delta}^{2}}{\mathbf{H}}^{\top}{\mathbf{H}})^{-1}$.

####  E.3.1 Euler-Maruyama solver

To compare the posterior distribution estimated by each algorithm with the
target posterior distribution, we use $10^{4}$ slices for the SW distance and
compare 1000 samples of the continuous SDEs defined by the TMPD, DTMPD, Song
et al. (2023) and Chung et al. (2022a) approximations obtained using 1000
Euler-Maruyama time-steps with 1000 samples of the true posterior
distribution.

Table 2 indicate the Central Limit Theorem (CLT) 95% confidence intervals
obtained by considering 20 randomly selected measurement models
(${\mathbf{H}},\sigma_{{\textnormal{y}}}$) for each setting
($d_{{\textnormal{x}}},d_{{\textnormal{x}}},\sigma_{{\textnormal{y}}}$).
Figure 2 shows the first two dimensions of the estimated posterior
distributions corresponding to the configurations ($80,1$) and ($800,1$) from
Table 1 for one of the randomly generated measurement model
(${\mathbf{H}},\sigma_{{\textnormal{y}}}$). Illustration of other settings are
given in Appendix E.3. These illustrations give us insight into the behaviour
of the algorithms and their ability to accurately estimate the posterior
distribution. We see that no one method estimates the posterior distribution
well.

Table 2: Sliced Wasserstein for the GMM case for the score-based SDEs discretized with Euler-Maruyama. |  | $\sigma_{{\mathbf{y}}}=0.01$ | $\sigma_{{\mathbf{y}}}=0.1$ | $\sigma_{{\mathbf{y}}}=1.0$
---|---|---|---|---
$d_{x}$ | $d_{y}$ | TMPD | DTMPD |  $\Pi$GDM | DPS | TMPD | DTMPD |  $\Pi$GDM | DPS | TMPD | DTMPD |  $\Pi$GDM | DPS
$8$ | $1$ | 1.5 $\pm$ 0.5 | 1.5 $\pm$ 0.5 | 1.5 $\pm$ 0.4 | 5.7 $\pm$ 2.2 | 1.4 $\pm$ 0.5 | 1.4 $\pm$ 0.5 | 1.2 $\pm$ 0.4 | 5.6 $\pm$ 2.1 | 0.9 $\pm$ 0.3 | 0.9 $\pm$ 0.3 | 0.9 $\pm$ 0.3 | 0.9 $\pm$ 0.3
$8$ | $2$ | 0.7 $\pm$ 0.3 | 3.2 $\pm$ 1.4 | 0.4 $\pm$ 0.3 | 6.2 $\pm$ 0.8 | 0.9 $\pm$ 0.3 | 2.7 $\pm$ 1.1 | 0.5 $\pm$ 0.3 | 6.2 $\pm$ 2.4 | 0.9 $\pm$ 0.2 | 1.8 $\pm$ 0.8 | 1.0 $\pm$ 0.3 | 1.2 $\pm$ 0.4
$8$ | $4$ | 0.3 $\pm$ 0.3 | 0.6 $\pm$ 0.4 | 0.1 $\pm$ 0.1 | - | 0.3 $\pm$ 0.2 | 0.7 $\pm$ 0.4 | 0.1 $\pm$ 0.0 | 8.4 $\pm$ 3.1 | 0.6 $\pm$ 0.2 | 0.9 $\pm$ 0.5 | 0.2 $\pm$ 0.1 | 0.3 $\pm$ 0.2
$80$ | $1$ | 2.7 $\pm$ 0.7 | 2.7 $\pm$ 0.7 | 2.9 $\pm$ 1.4 | 9.1 $\pm$ 1.3 | 2.3 $\pm$ 0.7 | 2.3 $\pm$ 0.7 | 2.1 $\pm$ 1.1 | 4.7 $\pm$ 1.8 | 1.5 $\pm$ 0.7 | 1.5 $\pm$ 0.7 | 1.8 $\pm$ 0.8 | 1.9 $\pm$ 0.9
$80$ | $2$ | 1.0 $\pm$ 0.5 | 3.3 $\pm$ 1.0 | 0.8 $\pm$ 0.7 | 2.2 $\pm$ 0.9 | 1.2 $\pm$ 0.5 | 3.3 $\pm$ 1.0 | 0.8 $\pm$ 0.7 | 6.0 $\pm$ 2.1 | 1.1 $\pm$ 0.2 | 2.2 $\pm$ 1.0 | 1.3 $\pm$ 0.5 | 1.5 $\pm$ 0.5
$80$ | $4$ | 0.3 $\pm$ 0.1 | 0.9 $\pm$ 0.5 | 0.1 $\pm$ 0.0 | - | 0.4 $\pm$ 0.2 | 1.0 $\pm$ 0.5 | 0.1 $\pm$ 0.1 | 4.4 $\pm$ 1.6 | 0.9 $\pm$ 0.2 | 1.0 $\pm$ 0.4 | 0.4 $\pm$ 0.2 | 0.5 $\pm$ 0.3
$800$ | $1$ | 3.1 $\pm$ 0.7 | 3.1 $\pm$ 0.7 | 3.2 $\pm$ 1.0 | 6.8 $\pm$ 1.2 | 2.9 $\pm$ 0.6 | 2.9 $\pm$ 0.6 | 2.8 $\pm$ 0.7 | 6.4 $\pm$ 1.5 | 1.5 $\pm$ 0.4 | 1.5 $\pm$ 0.4 | 1.3 $\pm$ 0.3 | 1.3 $\pm$ 0.3
$800$ | $2$ | 1.3 $\pm$ 0.4 | 3.6 $\pm$ 1.2 | 0.8 $\pm$ 0.5 | 7.4 $\pm$ 0.9 | 1.3 $\pm$ 0.3 | 3.2 $\pm$ 1.1 | 0.8 $\pm$ 0.4 | 6.4 $\pm$ 1.9 | 1.2 $\pm$ 0.3 | 1.9 $\pm$ 0.5 | 1.1 $\pm$ 0.3 | 1.1 $\pm$ 0.3
$800$ | $4$ | 0.3 $\pm$ 0.2 | 0.9 $\pm$ 0.6 | 0.6 $\pm$ 0.5 | - | 0.4 $\pm$ 0.2 | 0.9 $\pm$ 0.6 | 0.1 $\pm$ 0.0 | 5.8 $\pm$ 1.4 | 0.9 $\pm$ 0.2 | 1.1 $\pm$ 0.5 | 0.4 $\pm$ 0.2 | 0.4 $\pm$ 0.2

We now proceed to illustrate the first 2 components for one of the measurement
models for all the different combinations of
$(d_{{\textnormal{x}}},d_{{\textnormal{y}}})$ combinations used in Table 2.

###  E.4 Inpainting and super-resolution

Our experiments compare score-based methods (discretized with Euler-Maruyama)
and diffusion methods (DDPM or DDIM). Our experiments also compare the
different sampling methods discussed. In each comparison, we use the same
score network for each method and a comparable numerical method. All methods
are discretized using 1000 denoising steps. For the Markov chain methods we
use DDIM with $\eta=1.0$ and for the SDE methods we use an Euler-Maruyama
discretization. We show our full results in Fig. 5-12.

We found the hyperparameter suggested for the DDIM method for the VP SDE in
Song et al. (2023) $r_{t}=v_{t}$ to give unstable solutions for the Algorithm
given in Song et al. (2023). To make the method stable, we instead plug the
$\Pi$GDM posterior score approximation into a DDIM sampler in a similar way to
Algorithm 1, which, for the VPSDE, brings the algorithm $\Pi$GDM-D closer to
our method; we are then able to choose $r^{2}_{t}=v_{t}$ for both VP DDIM and
VP SDE methods. We are able to use the hyperparameter
$r^{2}_{t}=v_{t}/(1+v_{t})$ as suggested by Song et al. (2023) for both VE
DDIM and VE SDE methods. For DDPM we use the step-size constant suggested in
Chung et al. (2022a) for inpainting,
$\zeta_{i}=\zeta^{\prime}/\|{\mathbf{y}}-{\mathbf{H}}{\mathbf{m}}_{0|t}\|$,
where we use a value in the suggested range that worked well
$\zeta^{\prime}=0.4$, and don’t tune it further towards better results. In
some cases plugging in the DPS posterior score estimate to Euler-Maruyama
schemes lead to instable results in our experiments, usually when the
observation noise variance $\sigma_{y}^{2}$ was small. Consequently, we could
not always make a comparison to DPS for the Euler-Maruyama discretization of
the reverse-SDE.

Figure 5: Inpainting samples from the VE SDE on CIFAR-10. The observation
model was box mask with Gaussian ($\sigma_{y}=0.001$) noise. Figure 6:
Inpainting samples from the VP SDE on CIFAR-10. The observation model was box
mask with Gaussian ($\sigma_{y}=0.01$) noise. Figure 7: Inpainting samples
from the VE SDE on FFHQ. The observation model was box mask with Gaussian
($\sigma_{y}=0.001$) noise. Figure 8: Inpainting samples from the VE SDE on
FFHQ. The observation model was half mask with Gaussian ($\sigma_{y}=0.001$)
noise. Figure 9: Inpainting samples from the VE SDE on CelebA-HQ. The
observation model was a half mask with Gaussian ($\sigma_{y}=0.001$) noise.
Figure 10: Super-resolution samples from the VE SDE on CelebA-HQ. The
observation model was a bicubic downsampler with downsampling ratio of 16 (on
each axis), and Gaussian ($\sigma_{y}=0.001$) noise. Figure 11: Super-
resolution samples from the VP SDE on CIFAR-10. The observation model was a
nearest-neighbour downsampler with downsampling ratio of 2 (on each axis), and
Gaussian ($\sigma_{y}=0.1$) noise. Figure 12: Super-resolution samples from
the VP SDE on CIFAR-10. The observation model was a nearest-neighbour
downsampler with downsampling ratio of 2 (on each axis), and Gaussian
($\sigma_{y}=0.05$) noise.