# A Theoretical Justification for Image Inpainting using
Denoising Diffusion Probabilistic Models

Litu Rout<EMAIL_ADDRESS>The University of Texas at Austin  Advait
Parulekar<EMAIL_ADDRESS>The University of Texas at Austin  Constantine
Caramanis<EMAIL_ADDRESS>The University of Texas at Austin  Sanjay
Shakkottai<EMAIL_ADDRESS>The University of Texas at Austin

###### Abstract

We provide a theoretical justification for sample recovery using diffusion
based image inpainting in a linear model setting. While most inpainting
algorithms require retraining with each new mask, we prove that diffusion
based inpainting generalizes well to unseen masks without retraining. We
analyze a recently proposed popular diffusion based inpainting algorithm
called RePaint [LDRYTV22], and show that it has a bias due to misalignment
that hampers sample recovery even in a two-state diffusion process. Motivated
by our analysis, we propose a modified RePaint algorithm we call RePaint+ that
provably recovers the underlying true sample and enjoys a linear rate of
convergence. It achieves this by rectifying the misalignment error present in
drift and dispersion of the reverse process. To the best of our knowledge,
this is the first linear convergence result for a diffusion based image
inpainting algorithm.

##  1 Introduction

We study the mathematical principles that underlie the empirical success of
image inpainting using Denoising Diffusion Probabilistic Models (DDPMs)
[SWMG15, HJA20, SSKKEP20], which are the backbone of large-scale generative
models including DALL-E [RPGGVRCS21, RDNCC22], Imagen [SCSLWDGAML+22], and
Stable Diffusion [RBLEO22]. The goal of image inpainting is to reconstruct
missing parts that are semantically consistent with the known portions of an
image. A major challenge in this task is to generate parts that are consistent
with the available parts of the image; generative modeling is the key tool for
this inpainting goal.

Figure 1: Example demonstrating the bias in RePaint [LDRYTV22]. Starting from
the Gaussian prior $\overleftarrow{{\bm{x}}_{1}}$ (green circles), reverse SDE
as proposed by RePaint generates $\overleftarrow{{\bm{x}}_{0}}$ (brown circle)
that matches with the true data $\overrightarrow{{\bm{x}}_{0}}$ (blue square)
in known coordinates, but differs in the inpainted region. In the figure, the
blue squares along the $\sim 56^{\circ}$ line represents the true samples,
whereas the brown circles along the $\sim 10^{\circ}$ line represents the
samples recovered by RePaint. Note that the recovered samples match the true
samples along the x-coordinate (known data), but have a bias along the
y-coordinate (missing data).

Variational-Auto Encoders (VAEs) [KW13] and Generative Adversarial Networks
(GANs) [GPMXWOCB14] have been the basis for many successful inpainting
techniques in recent years. However, with the advent of Score based Generative
Models (SGMs) [SE19] or DDPMs [SWMG15, HJA20], current focus has shifted
towards an alternate paradigm of image inpainting. In this paradigm, most
methods fall into one of two categories. In the first category, one learns a
diffusion process specific to a downstream task, such as inpainting or super-
resolution [WDTSDM22, SCCLHSFN22]. In the second, a general purpose diffusion
based generative model is learned and the diffusion process is guided in the
inference phase catering to the downstream task [JADPDT21, CKJGY21, SSXE21,
DDDD22, KEES22]. In this paper, we analyze an inpainting approach [LDRYTV22]
that falls in the second category. We show that the generative prior of DDPM
is sufficient to fully reconstruct missing parts of an image.

DDPM represents a class of generative models that learn to diffuse a clean
image into tractable noise and then follow a reverse Markov process to produce
a clean image by progressively denoising pure noise [SWMG15, HJA20, SSKKEP20]
(see 2). Owing to their high expressive power, these models show appealing
results in conditional/unconditional image generation [HJA20, SSKKEP20, DN21,
KAAL22], text-to-image synthesis [RBLEO22], time series modeling [TSSE21],
audio synthesis [KPHZC21], image-to-image translation [SCCLHSFN22],
controllable text generation [LTGLH22], and image restoration [KEES22]. For
restoration tasks such as image inpainting, DDPMs are interesting because they
easily adapt to unknown tasks, such as new masks without having to go through
the entire retraining process. This gives DDPM an edge over prior restoration
techniques based upon GANs [WYWGLDQC18, YLYSLH18] and VAEs [RVV19, PLXL21].

However, one challenging aspect limiting their pervasive usage stems from the
choice of hyper-parameters in the reverse Markov process. Although the forward
diffusion process can be explicitly computed beforehand, the reverse process
is computationally very expensive as it requires sampling at every
intermediate state. Furthermore, popular diffusion based inpainting
algorithms, such as RePaint [LDRYTV22] require additional resampling at each
of these states to harmonize inpainted parts with the rest of the image. These
resampling steps increase the computational burden at the cost of semantically
meaningful reconstruction. Also, Figure 1 shows that RePaint generates biased
samples that hamper perfect recovery. Therefore, it is becoming increasingly
important to address these bottlenecks to facilitate their successful real-
world deployment. A persistent challenge in this respect stems from the fact
that despite the remarkable progress of diffusion based inpainting, our
theoretical understanding remains in its early stages.

To address these issues, one emerging line of theoretical research aims to
provide convergence guarantees for Score-based Generative Models (SGMs) that
constitute DDPMs [CCLLSZ22, LLT22, LLT22a, CLL22]. A crucial assumption of
these approaches is that the target distribution
$q(\overrightarrow{{\bm{x}}_{0}})$ assumes a density with respect to Lebesgue
measure. Another line of work [Bor22, Pid22] studies convergence of SGMs under
the famous manifold hypothesis [TSL00, FMN16]. Starting from a prior
$p(\overleftarrow{{\bm{x}}_{T}})$, these lines of work analyze convergence of
the reverse process $p_{{\bm{\theta}}}(\overleftarrow{{\bm{x}}_{0}})$ to the
target distribution $q(\overrightarrow{{\bm{x}}_{0}})$. Different from these
lines of research, our paper focuses on the convergence of resampling step
used to harmonize the inpainted image (3).

The assumption that the target distribution admits a density according
Lebesgue measure indicates that $q(\overrightarrow{{\bm{x}}_{0}})>0$ for all
$\overrightarrow{{\bm{x}}_{0}}\in\mathbb{R}^{d}$. In other words, every
$\overrightarrow{{\bm{x}}_{0}}$ is a possible sample drawn from
$q(\overrightarrow{{\bm{x}}_{0}})$, which contradicts the fact that most
natural images reside on a low dimensional manifold with compact support
[TSL00, FMN16]. This assumption would require a target distribution over
digits to allocate nonzero probability mass to very unlikely samples, such as
human faces, animals, and bedroom scenes. Although the prior analysis offers
some insights, it is of little practical significance in this manifold setting
[KSSKM22]. As a step towards circumventing this issue, [Pid22, Bor22] prove
that SGMs learn distributions supported on a low dimension substructure,
satisfying the manifold hypothesis.

For inpainting, however, we need an additional structure on the data
distribution similar to prior works in related disciplines [BJPD17, DGE18,
JLDC20, JADPDT21]. This is because our goal is not just to show minimum
discrepancy between $q(\overrightarrow{{\bm{x}}_{0}})$ and
$p_{\bm{\theta}}(\overleftarrow{{\bm{x}}_{0}})$ with respect to some
divergence ${\mathcal{D}}$, i.e.,
${\mathcal{D}}\left(p_{\bm{\theta}}(\overleftarrow{{\bm{x}}_{0}}),q(\overrightarrow{{\bm{x}}_{0}})\right)\leq\epsilon$,
but also to prove sample wise convergence, i.e.,
$\left\|\overleftarrow{{\bm{x}}_{0}}-\overrightarrow{{\bm{x}}_{0}}\right\|\leq\epsilon$.
It is worth mentioning that the problem of sample recovery is well studied in
optimization literature [Tib96, FNW07, BHSN10, FN03, BJPD17]. The focus of our
paper, instead, is to use this well understood setting to provide a
theoretical justification for sample recovery using a DDPM-based non-
optimization technique [LDRYTV22]. In this technique, we leverage the
expressive power of diffusion based generative models [HJA20] without
attempting to directly solve a constrained optimization problem. Following
conventional wisdom [BJPD17, JLDC20, JADPDT21], we implant a specific model on
the data generating distribution that satisfies the manifold hypothesis, and
show that the learning algorithm recovers this model by detecting the
underlying substructure: see Figure 2.

Figure 2: Comparison between RePaint ($\overleftarrow{\bar{{\bm{x}}}_{0}}$)
and RePaint+ ($\overleftarrow{\hat{{\bm{x}}}_{0}}$). RePaint+ discovers
$\epsilon$-accurate solutions whereas RePaint suffers from the misalignment
bias. Green circles indicate Gaussian prior $\overleftarrow{{\bm{x}}_{1}}$ of
the reverse Markov process. In the figure, the blue squares along the $\sim
56^{\circ}$ line represent the true samples. As before, the brown circles
along the $\sim 10^{\circ}$ line represents the samples recovered by RePaint
(this has a bias) and the orange stars along the $\sim 56^{\circ}$ line
represent the samples recovered by Repaint+. Note that the recovered samples
by RePaint+ has no observable bias (recovered samples are overlapping on the
true samples).

###  1.1 Contributions

Our main contribution is to provide a theoretical justification for diffusion
based image inpainting in a linear model setting. By analyzing diffusion over
two states, our analysis explains previously not understood phenomena.
Importantly, we derive algorithmic insights which, as we demonstrate later,
deliver improvement beyond the two-state processes.

First, unlike prior inpainting methods, we prove that diffusion based image
inpainting easily adapts to each new mask without retraining (Theorem 4 in 3).

Next, inspired by the empirical success of a recently proposed diffusion based
image inpainting method [LDRYTV22], we analyze its theoretical properties in a
system of two-state diffusion processes. We observe that there exists a bias
due to misalignment that hampers perfect recovery: see Figure 1. Our analysis
motivates us to rectify this misalignment that helps eliminate the bias. We
refer to this method as RePaint+. We provide a simplified version of RePaint+
in Algorithm 2 and a general version in Algorithm 3. We conduct toy
experiments in Appendix B.

In the linear model setting, we derive a closed-form solution of the
generative model using the transition kernels of DDPM (Theorem 3). Using this
solution, we prove that RePaint+ enjoys a linear rate of convergence in the
resampling phase (Theorem 4). This allows us to appropriately choose the
number of resampling rounds in practice. We generalize this notion of
convergence to a setting where the generative model is approximate (Theorem
6). Further, we justify the benefits of resampling over slowing down the
diffusion process in the context of inpainting (3.2.4).

Notation: We denote by bold upper-case letter ${\bm{A}}$ a matrix, bold lower-
case letter ${\bm{x}}$ a vector, and normal lower-case letter $x$ a scalar.
${\bm{I}}_{k}$ denotes a $k\times k$ identity matrix. The set of integers
$\\{1,\dots,N\\}$ are captured in $[N]$. Element-wise product is represented
by $\odot$. The operator ${\bm{D}}\left({\bm{x}}\right)$ diagonalizes a vector
${\bm{x}}$. $\left\|{\bm{A}}\right\|$ denotes the spectral norm and
$\left\|{\bm{A}}\right\|_{F}$, the Frobenius norm of a matrix ${\bm{A}}$. For
a vector, $\left\|{\bm{x}}\right\|$ denotes its Euclidean norm.
${\mathcal{U}}\left(\cdots\right)$ denotes uniform distribution.

##  2 Background on Diffusion Models

Like other generative models, such as GANs [GPMXWOCB14], VAEs [KW13], and Flow
[DSB17], Score-based generative models (SGM) [SE19] learn to sample from an
unknown distribution given a set of samples drawn from this distribution. In
this section, we first briefly introduce key ingredients of SGM. Then, we
detail the DDPM interpretation of SGM in 2.2, which holds the foundation of
RePaint+.

###  2.1 Score-based Generative Model

The central part of SGM consists of two Stochastic Differential Equations
(SDEs). The forward SDE (1), an Ornstein-Uhlenbeck (OU) process in its
simplest form, transforms the data distribution
$q(\overrightarrow{{\bm{x}}_{0}})$ to a reference distribution, which is
${\mathcal{N}}\left(\mathbf{0},{\bm{I}}\right)$ in most cases. Here, a sample
$\overrightarrow{{\bm{x}}_{t}}$ at time $t$ follows:

$\displaystyle
d\overrightarrow{{\bm{x}}_{t}}=-\frac{1}{2}\overrightarrow{{\bm{x}}_{t}}dt+dW_{t};~{}\overrightarrow{{\bm{x}}_{0}}\sim
q(\overrightarrow{{\bm{x}}_{0}}).$ (1)

With time rescaling parameters $\alpha_{t}$ and $\beta_{t}$, where
$\alpha_{t}=\int_{0}^{t}\beta_{s}ds$, and normalized time
$t=\\{0,\cdots,1\\}$, the transition kernel of (1) becomes:

$\displaystyle
q\left(\overrightarrow{{\bm{x}}_{t}}\mid\overrightarrow{{\bm{x}}_{0}}\right)={\mathcal{N}}\left(\overrightarrow{{\bm{x}}_{t}};\exp{(-\frac{\alpha_{t}}{2})\overrightarrow{{\bm{x}}_{0}},\left(1-\exp{(-\alpha_{t})}\right){\bm{I}}}\right),$

which leads to the well known form of forward SDE:

$\displaystyle
d\overrightarrow{{\bm{x}}_{t}}=-\frac{1}{2}\beta_{t}\overrightarrow{{\bm{x}}_{t}}+\sqrt{\beta_{t}}d\mathbf{B}_{t};~{}\overrightarrow{{\bm{x}}_{0}}\sim
q(\overrightarrow{{\bm{x}}_{0}}).$

Here, $\left\\{d\mathbf{B}_{t}\right\\}_{t\geq 0}$ represents the standard
Brownian motion in $\mathbb{R}^{d}$. The reverse SDE has a form similar to the
forward SDE with time reversal:

$\displaystyle
d\overleftarrow{{\bm{x}}_{t}}=\frac{1}{2}\beta_{1-t}\overleftarrow{{\bm{x}}_{t}}+\beta_{1-t}\nabla\log
q(\overrightarrow{{\bm{x}}_{t}}|\overrightarrow{{\bm{x}}_{0}})+\sqrt{\beta_{1-t}}d\mathbf{B}_{t},$

where $\overleftarrow{{\bm{x}}_{T}}\sim
p\left(\overleftarrow{{\bm{x}}_{T}}\right)\coloneqq{\mathcal{N}}\left(\mathbf{0},{\bm{I}}\right)$.
Usually, the Euler-Maryuama discretization scheme is employed while
implementing these SDEs in practice. Since we do not have access to
$q(\overrightarrow{{\bm{x}}_{0}})$, a neural network
$s_{\bm{\theta}}\left(\overrightarrow{{\bm{x}}_{t}},t\right)$ is used to
approximate the score function $\nabla\log
q(\overrightarrow{{\bm{x}}_{t}}|\overrightarrow{{\bm{x}}_{0}})$. Training is
performed by minimizing a score-matching objective [HD05, Vin11, SE19]:

$\displaystyle\min_{\theta}\mathop{\mathbb{E}}_{\begin{subarray}{c}\overrightarrow{{\bm{x}}_{0}}\sim
q\left(\overrightarrow{{\bm{x}}_{0}}\right)\\\
\overrightarrow{{\bm{x}}_{t}}\sim
q\left(\overrightarrow{{\bm{x}}_{t}}|\overrightarrow{{\bm{x}}_{0}}\right)\end{subarray}}\left[\left\|\nabla\log
q(\overrightarrow{{\bm{x}}_{t}}|\overrightarrow{{\bm{x}}_{0}})-s_{\bm{\theta}}\left(\overrightarrow{{\bm{x}}_{t}},t\right)\right\|^{2}\right].$

Next, we discuss DDPM and its connection with SGM.

###  2.2 Denoising Diffusion Probabilistic Model

Diffusion model [SWMG15, HJA20] is an emerging class of generative models that
share strikingly similar properties with score-based generative models [SE19].
It consists of two stochastic processes. First, the forward process (diffusion
process) gradually adds Gaussian noise to an image according to a variance
schedule. This is a Markov chain with stationary distribution typically set to
a tractable distribution that is easy to sample from, e.g., a standard
Gaussian ${\mathcal{N}}\left(\mathbf{0},{\bm{I}}\right)$. Second, the reverse
process (denoising process) learns to gradually denoise a sample drawn from
the tractable distribution. The denoising process continues until it produces
a high-quality image from the original data distribution.

####  2.2.1 Forward Process

Let $\overrightarrow{{\bm{x}}_{0}}\sim q(\overrightarrow{{\bm{x}}_{0}})$
denote a sample in $\mathbb{R}^{d}$ drawn from the data distribution
$q(\overrightarrow{{\bm{x}}_{0}})$. For discrete time steps $1\leq t\leq T$,
$q({\bm{x}}_{t})$ represents the distribution at time $t$. In the forward
process, the diffusion takes place according to a fixed Gaussian transition
kernel, i.e.,

$\displaystyle
q\left(\overrightarrow{{\bm{x}}_{1:T}}|\overrightarrow{{\bm{x}}_{0}}\right)\coloneqq\prod_{t=1}^{T}q(\overrightarrow{{\bm{x}}_{t}}|\overrightarrow{{\bm{x}}_{t-1}});~{}q(\overrightarrow{{\bm{x}}_{t}}|\overrightarrow{{\bm{x}}_{t-1}})\coloneqq\mathcal{N}\left(\overrightarrow{{\bm{x}}_{t}};\sqrt{1-\beta_{t}}\overrightarrow{{\bm{x}}_{t-1}},\beta_{t}{\bm{I}}_{d}\right),$
(2)

where $\beta_{1},\dots,\beta_{T}$ denote a deterministic variance schedule.
Thus, a sample from
$q(\overrightarrow{{\bm{x}}_{t}}|\overrightarrow{{\bm{x}}_{t-1}})$ is given by
$\overrightarrow{{\bm{x}}_{t}}=\sqrt{1-\beta_{t}}\overrightarrow{{\bm{x}}_{t-1}}+\sqrt{\beta_{t}}\mathbf{\mathbf{\epsilon}}$,
where
$\mathbf{\epsilon}\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}_{d}\right)$. The
first term is called drift and the second, dispersion. An important property
of the Gaussian diffusion process is that it has a simple form when sampling
from any intermediate time steps. For $\alpha_{t}\coloneqq 1-\beta_{t}$ and
$\bar{\alpha}_{t}\coloneqq\prod_{s=1}^{t}\alpha_{s}$, the conditional
probability at time $t$ becomes:

$\displaystyle
q\left(\overrightarrow{{\bm{x}}_{t}}|\overrightarrow{{\bm{x}}_{t-1}}\right)={\mathcal{N}}\left(\overrightarrow{{\bm{x}}_{t}};\sqrt{\bar{\alpha}_{t}}\overrightarrow{{\bm{x}}_{0}},\left(1-\bar{\alpha}_{t}\right){\bm{I}}_{d}\right).$

####  2.2.2 Reverse Process

The reverse process is a Markov chain initialized at the stationary
distribution of the forward process, i.e.,
$p(\overleftarrow{{\bm{x}}_{T}})={\mathcal{N}}\left(\mathbf{0},{\bm{I}}_{d}\right)$.
Unlike the forward process, the reverse process
$p_{\theta}\left(\overleftarrow{{\bm{x}}_{0:T}}\right)$ is generated from a
learned Gaussian kernel, i.e.,

$\displaystyle
p_{\theta}\left(\overleftarrow{{\bm{x}}_{0:T}}\right)=p(\overleftarrow{{\bm{x}}_{T}})\prod_{t=1}^{T}p_{\theta}\left(\overleftarrow{{\bm{x}}_{t-1}}|\overleftarrow{{\bm{x}}_{t}}\right);~{}p_{\theta}\left(\overleftarrow{{\bm{x}}_{t-1}}|\overleftarrow{{\bm{x}}_{t}}\right)\coloneqq{\mathcal{N}}\left(\overleftarrow{{\bm{x}}_{t-1}};\mu_{\theta}\left(\overleftarrow{{\bm{x}}_{t}},t\right),\Sigma_{\theta}\left(\overleftarrow{{\bm{x}}_{t}},t\right)\right),$
(3)

where $\mu_{\theta}\left(\overleftarrow{{\bm{x}}_{t}},t\right)$ and
$\Sigma_{\theta}\left(\overleftarrow{{\bm{x}}_{t}},t\right)$ are neural
networks parameterized by $\theta$ to predict the mean and the variance of
$p_{\theta}\left(\overleftarrow{{\bm{x}}_{t-1}}|\overleftarrow{{\bm{x}}_{t}}\right)$,
respectively.

####  2.2.3 Training and Inference

DDPM aims to maximize the likelihood of a sample generated by the reverse
process. As per equation (3), the probability assigned to such a sample is
obtained by marginalizing over the remaining random variables, i.e.,
$\overleftarrow{{\bm{x}}_{1}},\overleftarrow{{\bm{x}}_{2}},\cdots,\overleftarrow{{\bm{x}}_{T}}$
as denoted by the following expression:

$\displaystyle p(\overleftarrow{{\bm{x}}_{0}})=\int
p(\overleftarrow{{\bm{x}}_{0:T}})d\overleftarrow{{\bm{x}}_{1:T}},$

which is not easy to compute. It becomes tractable by considering relative
probability between the forward and the reverse processes [Jar97, SWMG15],
leading to the variational lower bound,

$\displaystyle\mathbb{E}\left[-\log
p_{\theta}\left(\overleftarrow{{\bm{x}}_{0}}\right)\right]\leq\mathbb{E}_{q}\left[-\log\frac{p_{\theta}\left(\overleftarrow{{\bm{x}}_{0:T}}\right)}{q\left(\overrightarrow{{\bm{x}}_{1:T}}|\overrightarrow{{\bm{x}}_{0}}\right)}\right]=\mathbb{E}_{q}\left[-\log
p\left(\overleftarrow{{\bm{x}}_{T}}\right)-\sum_{t\geq
1}\log\frac{p_{\theta}\left(\overleftarrow{{\bm{x}}_{t-1}}|\overleftarrow{{\bm{x}}_{t}}\right)}{q\left(\overrightarrow{{\bm{x}}_{t}}|\overrightarrow{{\bm{x}}_{t-1}}\right)}\right].$
(4)

This objective is further simplified to

$\displaystyle\mathbb{E}_{q}\big{[}{\mathcal{D}}_{KL}\left(q\left(\overrightarrow{{\bm{x}}_{T}}|\overrightarrow{{\bm{x}}_{0}}\right)\parallel
p\left(\overleftarrow{{\bm{x}}_{T}}\right)\right)+\sum_{t>1}{\mathcal{D}}_{KL}\left(q\left(\overrightarrow{{\bm{x}}_{t-1}}|\overrightarrow{{\bm{x}}_{t}},\overrightarrow{{\bm{x}}_{0}}\right)\parallel
p_{\theta}\left(\overleftarrow{{\bm{x}}_{t-1}}|\overleftarrow{{\bm{x}}_{t}}\right)\right)-\log
p_{\theta}\left(\overleftarrow{{\bm{x}}_{0}}|\overleftarrow{{\bm{x}}_{1}}\right)\Big{]},$
(5)

where
$q\left(\overrightarrow{{\bm{x}}_{t-1}}|\overrightarrow{{\bm{x}}_{t}},\overrightarrow{{\bm{x}}_{0}}\right)=\mathcal{N}\left(\overrightarrow{{\bm{x}}_{t-1}};\tilde{\mu}_{t}\left(\overrightarrow{{\bm{x}}_{t}},\overrightarrow{{\bm{x}}_{0}}\right),\tilde{\beta}_{t}{\bm{I}}_{d}\right)$,
$\tilde{\mu}_{t}\left(\overrightarrow{{\bm{x}}_{t}},\overrightarrow{{\bm{x}}_{0}}\right)\coloneqq\frac{\sqrt{\bar{\alpha}_{t-1}}\beta_{t}}{1-\bar{a}_{t}}\overrightarrow{{\bm{x}}_{0}}+\frac{\sqrt{\alpha}_{t}\left(1-\bar{\alpha}_{t-1}\right)}{1-\bar{\alpha}_{t}}\overrightarrow{{\bm{x}}_{t}}$
and
$\tilde{\beta}_{t}{\bm{I}}_{d}\coloneqq\frac{1-\bar{\alpha}_{t-1}}{1-\bar{\alpha}_{t}}\beta_{t}$.
Equation (5) contains three crucial terms. The first term measures the
divergence between the stationary distribution of the forward process
conditioned on a clean sample $\overrightarrow{{\bm{x}}_{0}}\sim
q\left(\overrightarrow{{\bm{x}}_{0}}\right)$ and the initial distribution of
the reverse process. The divergence is negligible in practice, thanks to the
exponential convergence of OU processes. This is ignored as there are no
trainable parameters, ${\bm{\theta}}$.

The second term measures the divergence between the forward and the reverse
process at intermediate time steps, $1<t\leq T$. Since the conditional
probabilities are Gaussians, this can be explicitly computed as:

$\displaystyle\mathbb{E}_{q}\left[\frac{1}{2\beta_{t}^{2}}\left\|\tilde{\mu}_{t}\left(\overrightarrow{{\bm{x}}_{t}},\overrightarrow{{\bm{x}}_{0}}\right)-\mu_{\theta}\left(\overleftarrow{{\bm{x}}_{t}},t\right)\right\|^{2}\right]+Constant,$
(6)

where
$p_{\theta}\left(\overleftarrow{{\bm{x}}_{t-1}}|\overleftarrow{{\bm{x}}_{t}}\right)={\mathcal{N}}\left(\overleftarrow{{\bm{x}}_{t-1}};\mu_{\theta}\left(\overleftarrow{{\bm{x}}_{t}},t\right),\beta_{t}{\bm{I}}\right)$.
In the training phase, $\overleftarrow{{\bm{x}}_{t}}$ is replaced with
$\overrightarrow{{\bm{x}}_{t}}$ since the first term in (5) is negligible.
Further, reparameterization of the posterior mean as
$\mu_{\theta}\left(\overrightarrow{{\bm{x}}_{t}},t\right)=\frac{1}{\sqrt{\alpha_{t}}}\left(\overrightarrow{{\bm{x}}_{t}}-\frac{\beta_{t}}{\sqrt{1-\bar{\alpha}_{t}}}~{}\mathbf{\epsilon}_{\theta}\left(\overrightarrow{{\bm{x}}_{t}},t\right)\right)$
yields better results and a simplified loss [HJA20]:

$\displaystyle\mathbb{E}_{{\bm{x}}_{0},\mathbf{\epsilon}}\left[\frac{1}{2\alpha_{t}\left(1-\bar{\alpha}_{t}\right)}\left\|\mathbf{\epsilon}-\mathbf{\epsilon}_{\theta}\left(\sqrt{\bar{\alpha}_{t}}~{}\overrightarrow{{\bm{x}}_{0}}+\sqrt{1-\bar{\alpha}_{t}}~{}\mathbf{\epsilon},t\right)\right\|^{2}\right].$
(7)

During inference, for
$\mathbf{\epsilon}\sim{\mathcal{N}}\left(0,{\bm{I}}\right)$, a sample is
generated by
$\overleftarrow{{\bm{x}}_{t-1}}=\frac{1}{\sqrt{\alpha_{t}}}\left(\overleftarrow{{\bm{x}}_{t}}-\frac{\beta_{t}}{\sqrt{1-\bar{\alpha}_{t}}}~{}\mathbf{\epsilon}_{\theta}\left(\overleftarrow{{\bm{x}}_{t}},t\right)\right)+\sqrt{\beta_{t}}~{}\mathbf{\epsilon}$.

The third term is the standard maximum likelihood estimator. With this
reparameterization, the problem of learning probability is converted to a
practically implementable minimum mean squared error (MMSE) problem. The goal
of training DDPMs is to obtain an optimal $\theta^{*}$ that solves the MMSE
problem with posterior mean (6) or noise (7).

During inference, the reverse process is initialized at the stationary
distribution of the forward process, usually
${\mathcal{N}}\left(\mathbf{0},{\bm{I}}\right)$. Then, the reverse Gaussian
transition kernel is followed using
$\mu_{{\bm{\theta}}^{*}}\left(\overleftarrow{{\bm{x}}_{t}},t\right)$. In the
final step,
$\mu_{{\bm{\theta}}^{*}}\left(\overleftarrow{{\bm{x}}_{1}},1\right)$ is
displayed without noise, which is a sample from the original data
distribution. Ignoring some constant factors, the forward and reverse SDEs of
SGM and the forward and reverse processes of DDPM are idential if we use
Gaussian transition kernel. Throughout this paper, we interchangeably use
these terms.

Input: Initialized weights ${\bm{\theta}}$, Stepsize $\eta$

repeat

Draw $\overrightarrow{{\bm{x}}_{0}}\sim q(\overrightarrow{{\bm{x}}_{0}})$,
$\epsilon\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}_{d}\right)$,

$t\sim Uniform\left(\\{1,\dots,T\\}\right)$

1 Loss (6):
${\mathcal{L}}\left({\bm{\theta}}\right)=\left\|\tilde{\mu}_{t}\left(\overrightarrow{{\bm{x}}_{t}},\overrightarrow{{\bm{x}}_{0}}\right)-\mu_{\theta}\left(\overleftarrow{{\bm{x}}_{t}},t\right)\right\|^{2}$
or

(7):
${\mathcal{L}}\left({\bm{\theta}}\right)=\left\|\mathbf{\epsilon}-\mathbf{\epsilon}_{\theta}\left(\sqrt{\bar{\alpha}_{t}}~{}\overrightarrow{{\bm{x}}_{0}}+\sqrt{1-\bar{\alpha}_{t}}~{}\mathbf{\epsilon},t\right)\right\|^{2}$

Run gradient descent
${\bm{\theta}}_{t+1}={\bm{\theta}}_{t}-\eta\nabla{\mathcal{L}}\left({\bm{\theta}}\right)$

until Convergence

Output: Trained weights ${\bm{\theta}}$

Algorithm 1 DDPM Training

Input: Given $\overrightarrow{{\bm{x}}_{0}}\sim
q\left(\overrightarrow{{\bm{x}}_{0}}\right)$, Inpainting mask ${\bm{m}}$, DDPM
trained weights ${\bm{\theta}}$, Fixed variance $\beta>0$

Draw
$\overleftarrow{{\bm{x}}_{1}}\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}_{d}\right)$

for $r=1,\dots,R$ do

Drift from reverse Markov process (3):
$\mu_{{\bm{\theta}}}\left(\overleftarrow{{\bm{x}}_{1}}\right)$

Alignment of drift:
$\overleftarrow{{\bm{x}}_{0}}=\mu_{{\bm{\theta}}}\left(\overleftarrow{{\bm{x}}_{1}}\right)\times\frac{1}{\sqrt{1-\beta}}$

Inpainting:
$\overleftarrow{{\bm{x}}_{0}}={\bm{m}}\odot\overleftarrow{{\bm{x}}_{0}}+\left(1-{\bm{m}}\right)\odot\overrightarrow{{\bm{x}}_{0}}$

Update $\overleftarrow{{\bm{x}}_{1}}\leftarrow\overleftarrow{{\bm{x}}_{0}}$

end for

Drift from reverse Markov process (3):
$\mu_{{\bm{\theta}}}\left(\overleftarrow{{\bm{x}}_{0}}\right)$

Alignment of drift:
$\overleftarrow{\hat{{\bm{x}}}_{0}}=\mu_{{\bm{\theta}}}\left(\overleftarrow{{\bm{x}}_{0}}\right)\times\frac{1}{\sqrt{1-\beta}}$

Output: $\overleftarrow{\hat{{\bm{x}}}_{0}}$

Algorithm 2 RePaint+ Inference Special Case

##  3 Theoretical Results

In this section, we present our main generative and inpainting results, first
for a two-state model (3.2.1, 3.2.2, 3.2.3), and subsequently generalize to a
multi-state diffusion model (3.2.4).

###  3.1 Problem Setup

The goal is to minimize the divergence
${\mathcal{D}}\left(p_{\bm{\theta}}(\overleftarrow{{\bm{x}}_{0}}),q(\overrightarrow{{\bm{x}}_{0}})\right)$,
where ${\mathcal{D}}(.,.)$ may be a total variation (TV) distance [CCLLSZ22],
KL-divergence [HJA20], or Wasserstein distance [Bor22]. For a generative
model, the objective is met as long as the distribution of
$\overleftarrow{{\bm{x}}_{0}}$ matches with
$q(\overrightarrow{{\bm{x}}_{0}})$. For diffusion based inpainting, however,
it is necessary to generate an $\epsilon$-close image, i.e.,
$\left\|\overleftarrow{{\bm{x}}_{0}}-\overrightarrow{{\bm{x}}_{0}}\right\|\leq\epsilon$
and $\overleftarrow{{\bm{x}}_{0}}$ lies in the support of
$q(\overrightarrow{{\bm{x}}_{0}})$. Suppose the data distribution
$q\left(\overrightarrow{{\bm{x}}_{0}}\right)$ is supported on a
$k$-dimensional subspace of $\mathbb{R}^{d}$. For
${\bm{A}}\in\mathbb{R}^{d\times k}$, we have the original samples
$\overrightarrow{{\bm{x}}_{0}}={\bm{A}}{\bm{z}}_{0}$, where
${\bm{z}}_{0}\in\mathbb{R}^{k}$ is distributed according to
${\mathcal{N}}\left(\mathbf{0},{\bm{I}}_{k}\right)$ and ${\bm{A}}$ is full
rank, i.e., $rank({\bm{A}})=k\leq d$. We ask if a model trained for a
generative modeling task helps recover the missing parts in an image
inpainting task. As noted earlier in 1, recovery under this model using
optimization techniques is well understood. Our goal is to use this well
established framework as a vehicle to analyze and provide new insights into
sample recovery using DDPMs.

Input: Given $\overrightarrow{{\bm{x}}_{0}}\sim
q\left(\overrightarrow{{\bm{x}}_{0}}\right)$, DDPM trained weights
${\bm{\theta}}$, Initial variance $\beta_{0}=0$, Variance schedule
$\\{\beta_{t}\\}_{t\geq 0}$, Alignment coefficients for drift
($\\{\omega_{t}\\}_{t\geq 0}$) and dispersion ($\\{\xi_{t}\\}_{t\geq 0}$),
Initial
$\overleftarrow{{\bm{x}}_{T}}\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}_{d}\right)$

for $t=T,\dots,1$ do

for $r=1,\dots,R$ do

if $t>1$ then

$\overrightarrow{\mathbf{\epsilon}}\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}_{d}\right)$
and
$\overleftarrow{\mathbf{\epsilon}}\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}_{d}\right)$

else

$\overrightarrow{\epsilon}=0$ and $\overleftarrow{\epsilon}=0$

end if

Known part:
$\overrightarrow{{\bm{x}}_{t-1}}=\sqrt{\bar{\alpha}_{t-1}}\overrightarrow{{\bm{x}}_{0}}+\sqrt{1-\bar{\alpha}_{t-1}}\overrightarrow{\mathbf{\epsilon}}$

Drift
$\mu_{{\bm{\theta}}}\left(\overleftarrow{{\bm{x}}_{t}},t\right)=\frac{1}{\sqrt{\alpha_{t}}}\left(\overleftarrow{{\bm{x}}_{t}}-\frac{\beta_{t}}{\sqrt{1-\bar{\alpha}_{t}}}~{}\mathbf{\epsilon}_{\theta}\left(\overleftarrow{{\bm{x}}_{t}},t\right)\right)$
and dispersion $\sqrt{\beta_{t}}\overleftarrow{\mathbf{\epsilon}}$ from
reverse process (3)

Alignment of drift:
$\mu_{{\bm{\theta}}}\left(\overleftarrow{{\bm{x}}_{t}},t\right)\times\omega_{t}$

Alignment of dispersion:
$\sqrt{\beta_{t}}\overleftarrow{\mathbf{\epsilon}}\times\xi_{t}$

Compute
$\overleftarrow{{\bm{x}}_{t-1}}=\mu_{{\bm{\theta}}}\left(\overleftarrow{{\bm{x}}_{t}},t\right)\times\omega_{t}+\sqrt{\beta_{t}}\overleftarrow{\mathbf{\epsilon}}\times\xi_{t}$

Inpainting:
$\overleftarrow{{\bm{x}}_{t-1}}={\bm{m}}\odot\overleftarrow{{\bm{x}}_{t-1}}+\left(1-{\bm{m}}\right)\odot\overrightarrow{{\bm{x}}_{t-1}}$

if $r<R$ and $t>1$ then

Draw $\mathbf{\epsilon}\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}_{d}\right)$

Pushforward:
$\overleftarrow{{\bm{x}}_{t}}=\sqrt{1-\beta_{t}}\overleftarrow{{\bm{x}}_{t-1}}+\sqrt{\beta_{t}}\mathbf{\epsilon}$

else

Update $\overleftarrow{{\bm{x}}_{1}}\leftarrow\overleftarrow{{\bm{x}}_{0}}$

end if

end for

end for

Drift
$\mu_{{\bm{\theta}}}\left(\overleftarrow{{\bm{x}}_{1}},1\right)=\frac{1}{\sqrt{\alpha_{1}}}\left(\overleftarrow{{\bm{x}}_{1}}-\frac{\beta_{1}}{\sqrt{1-\bar{\alpha}_{1}}}~{}\mathbf{\epsilon}_{\bm{\theta}}\left(\overleftarrow{{\bm{x}}_{1}},1\right)\right)$
from reverse Markov process (3)

Alignment of drift:
$\mu_{{\bm{\theta}}}\left(\overleftarrow{{\bm{x}}_{1}},1\right)\times\omega_{1}$

Compute
$\overleftarrow{\hat{{\bm{x}}}_{0}}=\mu_{{\bm{\theta}}}\left(\overleftarrow{{\bm{x}}_{1}},1\right)\times\omega_{1}$

Output: $\overleftarrow{\hat{{\bm{x}}}_{0}}$

Algorithm 3 RePaint+ Inference General Case

###  3.2 Major Insights from Analysis

Our main result is that given the weights of a DDPM trained for generative
modeling tasks, we can repurpose its objective to image inpainting without
retraining for the inpainting task and show perfect recovery with linear rate
of convergence. Under Assumption 1, we first show that the weights learned by
DDPM training Algorithm 1 successfully capture the underlying data generating
distribution, see Theorem 3. Then, we prove in Theorem 4 under Assummption 1
and 2 that the optimal solution of DDPM training Algorithm 1 generalizes well
to image inpainting tasks with unknown masks, see Algorithm 2 and Algorithm 3.
Similar to the variance schedule $\\{\beta_{t}\\}_{t\geq 0}$, the alignment
coefficients for drift $\\{\omega_{t}\\}_{t\geq 0}$ and dispersion
$\\{\xi_{t}\\}_{t\geq 0}$ can be computed beforehand. For instance, in
Algorithm 2, we choose $\beta_{t}=\beta$ and
$\omega_{t}=1/\sqrt{1-\beta}$111The proposed modifications over RePaint
[LDRYTV22] are highlighted in Algorithm 2 and Algorithm 3..

######  Assumption 1.

The column vectors of data generating model ${\bm{A}}$ are orthonormal, i.e.,
${\bm{A}}^{T}{\bm{A}}={\bm{I}}_{k}$.

######  Assumption 2.

Given an inpainting mask ${\bm{m}}\in\\{0,1\\}^{d}$, the following holds true:
${\bm{A}}{\bm{A}}^{T}{\bm{D}}({\bm{m}})\prec{\bm{I}}_{d}$.

Assumption 1 is a mild assumption used to simplify the expressions. For
Assumption 2, recall that ${\bm{D}}({\bm{m}})$ is a $d\times d$ diagonal
matrix (a mask) with elements in diagonal entries set to $1$ whenever data is
missing (i.e., data is masked) and $0$ otherwise (data is available). The
symbol $\prec$ indicates that the spectral norm of
${\bm{A}}{\bm{A}}^{T}{\bm{D}}({\bm{m}})$, denoted by $\lambda_{\max}$, is
strictly less than 1.

Assumption 2 has a physical interpretation that agrees with the intuition.
Since masking an image reduces its free energy, it is important that the mask
${\bm{D}}({\bm{m}})$ is well behaved so that sufficient energy is left in the
masked image for a faithful reconstruction. If ${\bm{D}}({\bm{m}})={\bm{I}}$,
then all the energy of the original signal is lost. In this scenario, it is
impossible to recover the original signal. When
${\bm{D}}({\bm{m}})=\mathbf{0}$, it is a trivial case because the algorithm
has access to the original signal itself. Our results on recoverability hold
for all the remaining cases when ${\bm{m}}\sim\\{0,1\\}^{d}$ such that
$\lambda_{\max}<1$.

An interesting avenue for further research is to precisely characterize how
much information we need in the form of ${\bm{D}}({\bm{m}})$ for perfect
recovery. This problem has been extensively studied in compressed sensing
literature [BJPD17, WRL19], which is not the focus of this paper. We defer
such questions to future work.

####  3.2.1 Generative Modeling using DDPM

Here, we present Theorem 3 for computing the analytical solution of DDPM for a
two-state model. In general, computing the analytical solution is a futile
exercise as we do not know $q(\overrightarrow{{\bm{x}}_{0}})$ a priori.
However, the additional structure on the data generating distribution, as
typically assumed in downstream tasks such as image inpainting, allows us to
derive an explicit form of the solution [BJPD17, WRL19]. We consider the
objective function (6) that estimates the mean of the posterior. One may wish
to estimate the noise instead with the reparameterized objective (7). Both the
choices produce similar experimental results [HJA20].

######  Theorem 3 (Generative Modeling).

Suppose Assumption 1 holds. Let us denote
${\bm{\theta}}^{*}=\arg_{\min}{\mathcal{L}}\left({\bm{\theta}}\right)$, where
${\mathcal{L}}\left({\bm{\theta}}\right)$ is defined as:

$\mathbb{E}_{\overrightarrow{{\bm{x}}_{0}},\overrightarrow{\mathbf{\epsilon}}}\left[\left\|\tilde{\mu}_{1}\left(\overrightarrow{{\bm{x}}_{1}}(\overrightarrow{{\bm{x}}_{0}},\overrightarrow{\mathbf{\epsilon}}),\overrightarrow{{\bm{x}}_{0}}\right)-\mu_{\theta}\left(\overrightarrow{{\bm{x}}_{1}}\left(\overrightarrow{{\bm{x}}_{0}},\overrightarrow{\mathbf{\epsilon}}\right)\right)\right\|^{2}\right].$

For a fixed variance $\beta>0$, if we consider a function approximator
$\mu_{\bm{\theta}}\left(\overrightarrow{{\bm{x}}_{1}}\left(\overrightarrow{{\bm{x}}_{0}},\overrightarrow{\mathbf{\epsilon}}\right)\right)={\bm{\theta}}\overrightarrow{{\bm{x}}_{1}}\left(\overrightarrow{{\bm{x}}_{0}},\overrightarrow{\mathbf{\epsilon}}\right)$,
then the closed-form solution
${\bm{\theta}}^{*}=\sqrt{1-\beta}{\bm{A}}{\bm{A}}^{T}$, which upon
renormalization by $\left(1/\sqrt{1-\beta}\right)$ recovers the true subspace
of $q\left(\overrightarrow{{\bm{x}}_{0}}\right)$.

###### Proof.

The proof is included in Appendix A.1. ∎

For any
$\overleftarrow{{\bm{x}}_{1}}\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}_{d}\right)$,
the reverse process with the optimal solution ${\bm{\theta}}^{*}$ generates
$\overleftarrow{{\bm{x}}_{0}}=\sqrt{1-\beta}{\bm{A}}{\bm{A}}^{T}\overleftarrow{{\bm{x}}_{1}}$.
After renormalization, it gives
$\overleftarrow{{\bm{x}}_{0}}={\bm{A}}\left({\bm{A}}^{T}\overleftarrow{{\bm{x}}_{1}}\right)$
which is unbiased and has identity covariance, i.e.,
$\mathbb{E}\left[{\bm{A}}^{T}\overleftarrow{{\bm{x}}_{0}}\right]=\mathbf{0}$
and
$\mathbb{E}\left[\left({\bm{A}}^{T}\overleftarrow{{\bm{x}}_{1}}\right)\left({\bm{A}}^{T}\overleftarrow{{\bm{x}}_{1}}\right)^{T}\right]=\mathbb{E}\left[{\bm{A}}^{T}\left(\overleftarrow{{\bm{x}}_{1}}\overleftarrow{{\bm{x}}_{1}}^{T}\right){\bm{A}}\right]={\bm{I}}_{k}$.
This verifies that Algorithm 1 recovers the subspace of the underlying data
distribution, which is supported on a linear manifold with Gaussian marginals.
In what follows, we provide a detailed analysis of the inpainting Algorithm 2.

####  3.2.2 Image Inpainting using RePaint+

Continuing with the two-state model, there are essentially three key
ingredients to diffusion based image inpainting. First, we initialize the
reverse process at
$\overleftarrow{{\bm{x}}_{1}}\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}_{d}\right)$.
Using the reverse Gaussian transition kernel obtained in Theorem 3, we
generate
$\overleftarrow{{\bm{x}}_{0}}=\left(1/\sqrt{1-\beta}\right)\mu_{{\bm{\theta}}^{*}}\left(\overleftarrow{{\bm{x}}_{1}}\right)$.
As per Theorem 3, $\overleftarrow{{\bm{x}}_{0}}$ is a sample that lies on the
manifold. However, we note that $\overleftarrow{{\bm{x}}_{0}}$ is not the
sample we are looking for as it could potentially lie far away from
$\overrightarrow{{\bm{x}}_{0}}$.

Second, we replace certain parts of $\overleftarrow{{\bm{x}}_{0}}$ as per the
given mask ${\bm{m}}$ with known information from
$\overrightarrow{{\bm{x}}_{0}}$ using the following rule:
$\overrightarrow{{\bm{x}}_{0}^{1}}={\bm{m}}\odot\overleftarrow{{\bm{x}}_{0}}+(\bm{1}-{\bm{m}})\odot\overrightarrow{{\bm{x}}_{0}}$.
This produces an arbitrary sample
$\overrightarrow{{\bm{x}}_{0}^{1}}\in\mathbb{R}^{d}$. Since the reverse SDE
pushes an arbitrary sample in $\mathbb{R}^{d}$ onto the manifold, we feed
$\overrightarrow{{\bm{x}}_{0}^{1}}$ as an input to the reverse SDE in the
first resampling round. Alternatively, one may pass
$\overrightarrow{{\bm{x}}_{0}^{1}}$ through the forward SDE without adding any
noise. In this case, the generated sample $\overrightarrow{{\bm{x}}_{1}^{1}}$
serves as the input to the reverse SDE. This provides just an additional
scaling that is easy to handle in our analysis.

Finally, the reverse SDE generates $\overleftarrow{{\bm{x}}_{0}^{1}}$, which
upon realignment is expected to be $\mathbf{\epsilon}$-close to the original
image, i.e.,
$\left\|\overleftarrow{{\bm{x}}_{0}^{1}}-\overrightarrow{{\bm{x}}_{0}}\right\|\leq\epsilon$.
This completes one resampling step.

Usually, one resampling step is not sufficient to obtain an
$\mathbf{\epsilon}$-close solution. We ask how many resampling steps Algorithm
2 needs to generate a satisfactory inpainted image. Furthermore, what is the
rate at which it converges to the $\mathbf{\epsilon}$-neighborhood of the
original image. We answer these questions favorably in Theorem 4.

######  Theorem 4 (Image Inpainting).

Let Assumption 1 and Assumption 2 hold. Suppose
$\lambda_{\max}\coloneqq\left\|{\bm{A}}{\bm{A}}^{T}{\bm{D}}({\bm{m}})\right\|$
and ${\bm{\theta}}^{*}=\arg_{\min}{\mathcal{L}}\left({\bm{\theta}}\right)$,
where ${\mathcal{L}}\left({\bm{\theta}}\right)$ is defined as:

$\mathbb{E}_{\overrightarrow{{\bm{x}}_{0}},\overrightarrow{\mathbf{\epsilon}}}\left[\left\|\tilde{\mu}_{1}\left(\overrightarrow{{\bm{x}}_{1}}(\overrightarrow{{\bm{x}}_{0}},\overrightarrow{\mathbf{\epsilon}}),\overrightarrow{{\bm{x}}_{0}}\right)-\mu_{\theta}\left(\overrightarrow{{\bm{x}}_{1}}\left(\overrightarrow{{\bm{x}}_{0}},\overrightarrow{\mathbf{\epsilon}}\right)\right)\right\|^{2}\right].$

For any mask ${\bm{m}}\sim\\{0,1\\}^{d}$, a fixed variance $\beta>0$, a
partially known image $\overrightarrow{{\bm{x}}_{0}}\sim
q\left(\overrightarrow{{\bm{x}}_{0}}\right)$, and reverse SDE initialized at
$\overleftarrow{{\bm{x}}_{1}}\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}_{d}\right)$,
the maximum number of resampling steps ($r$) needed by Algorithm 2 to recover
an $\epsilon$-accurate inpainted image, i.e.,

$\left\|\overleftarrow{{\bm{x}}_{0}^{r}}-\overrightarrow{{\bm{x}}_{0}}\right\|\leq\lambda_{\max}^{r}\biggl{(}\frac{\left\|{\bm{\theta}}^{*}\right\|\left\|\overleftarrow{{\bm{x}}_{1}}-\overrightarrow{{\bm{x}}_{0}}\right\|}{\sqrt{1-\beta}}\biggr{)}\leq\epsilon$

is upper bounded by

$\displaystyle{\mathcal{O}}\left(\frac{\log\left(\frac{\left\|\theta^{*}\right\|\left\|\overleftarrow{x_{1}}-\overrightarrow{{\bm{x}}_{0}}\right\|}{\epsilon\sqrt{1-\beta}}\right)}{\log\left(\frac{1}{\lambda_{\max}}\right)}\right).$

###### Proof.

The proof is included in Appendix A.2. ∎

We draw several key insights from Theorem 4.

Universal mask prinicple: An important observation is that the rate does not
depend upon the mask ${\bm{m}}$ as long as it is a valid inpainting mask.
Thus, Algorithm 2 recovers the original sample irrespective of whether it has
been trained on such masked images or not. This is important because it allows
us to repurpose the objective of a diffusion based generative model to address
image inpainting.

Linear rate of convergence:  A major implication of Theorem 4 is that
Algorithm 2 enjoys a linear rate of convergence. For this reason, we only need
a small increase in the number of resampling steps to get a significant
improvement in terms of $\epsilon$-accuracy.

Information bottleneck: One interesting controlling parameter for the rate of
convergence is $\lambda_{\max}$. From Theorem 4, it is evident that a large
value of $\lambda_{\max}$ requires more resampling steps. Since
$\lambda_{\max}$ correlates with missing information, it is understandable
that Algorithm 2 needs more iterations for a reasonable harmonization.

Low norm solution: Among other controlling parameters, Algorithm 2 seems to
favor low norm solutions from generative modeling. This indicates that having
an inductive bias in DDPM training Algorithm 1 to prefer a low norm solution
may assist in diffusion based inpainting.

Distance from initialization:  Finally, the distance from the original sample,
i.e., $\left\|\overleftarrow{x_{1}}-\overrightarrow{{\bm{x}}_{0}}\right\|$
increases the iteration complexity only logarithmically.

####  3.2.3 Image Inpainting with Noisy Generator

An immediate consequence of Theorem 4 follows when Algorithm 1 returns an
approximate solution
$\hat{{\bm{\theta}}}\coloneqq{\bm{\theta}}^{*}+\mathbf{\delta}\sqrt{1-\beta}$
in contrast to the exact analytical solution ${\bm{\theta}}^{*}$. This is a
reasonable setting because the closed-form solution is not exactly computable
when the data generating distribution
$q\left(\overrightarrow{{\bm{x}}_{0}}\right)$ is not known a priori, which is
the case for most interesting practical applications. We present Theorem 6
under Assumption 5 taking into account the approximate error in generative
modeling.

######  Assumption 5.

The perturbation $\mathbf{\delta}$ of the approximate solution is such that
$\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}({\bm{m}})\prec{\bm{I}}_{d}$.

######  Theorem 6.

Suppose Assumption 5 holds and
$\hat{\lambda}_{\max}\coloneqq\left\|\left({\bm{A}}{\bm{A}}^{T}+\delta\right){\bm{D}}({\bm{m}})\right\|$.
Then, running Algorithm 2 for $r$ resampling steps with a $\delta$-approximate
model $\hat{{\bm{\theta}}}$ yields:

$\displaystyle\left\|\overleftarrow{{\bm{x}}_{0}^{r}}-\overrightarrow{{\bm{x}}_{0}}\right\|\leq\hat{\lambda}_{\max}^{r}\biggl{(}\frac{\left\|{\bm{\theta}}^{*}\right\|\left\|\overleftarrow{{\bm{x}}_{1}}-\overrightarrow{{\bm{x}}_{0}}\right\|}{\sqrt{1-\beta}}+\frac{\left\|\hat{{\bm{\theta}}}-{\bm{\theta}}^{*}\right\|\left\|\overrightarrow{{\bm{x}}_{0}}\right\|}{\left(1-\hat{\lambda}_{\max}\right)\sqrt{1-\beta}}\biggr{)}+\frac{\left\|\hat{{\bm{\theta}}}-{\bm{\theta}}^{*}\right\|\left\|\overrightarrow{{\bm{x}}_{0}}\right\|}{\left(1-\hat{\lambda}_{\max}\right)\sqrt{1-\beta}}.$

###### Proof.

The proof is included in Appendix A.3. ∎

Theorem 6 states that since $\hat{\lambda}_{\max}^{r}<1$, we recover a
$\zeta$-approximate solution in the limit, where
$\zeta=\frac{\left\|\hat{{\bm{\theta}}}-{\bm{\theta}}^{*}\right\|\left\|\overrightarrow{{\bm{x}}_{0}}\right\|}{\left(1-\hat{\lambda}_{\max}^{r}\right)\sqrt{1-\beta}}$.
The reconstruction error is proportional to the approximation error of DDPM in
generative modeling. In the absence of noise, we have
$\left\|\hat{{\bm{\theta}}}-{\bm{\theta}}^{*}\right\|\rightarrow 0$, which
leads to recovery of the true underlying sample
$\overrightarrow{{\bm{x}}_{0}}$.

######  Corollary 7.

Suppose $q\left(\overrightarrow{{\bm{x}}_{0}}\right)$ has a compact support
with $\left\|\overrightarrow{{\bm{x}}_{0}}\right\|\leq\kappa$. If
$\delta={\mathcal{O}}\left(\epsilon\left(1-\hat{\lambda}_{\max}\right)/\kappa\right)$,
then

$\displaystyle\mathop{\sup}_{\overrightarrow{{\bm{x}}_{0}}\sim
q(\overrightarrow{{\bm{x}}_{0}})}\left\|\overleftarrow{{\bm{x}}_{0}^{r}}-\overrightarrow{{\bm{x}}_{0}}\right\|\leq\epsilon.$

###### Proof.

The proof is included in Appendix A.4. ∎

In practice, most of the data generating distributions have a compact support.
For these distributions, it is reasonable to have
$\left\|\overrightarrow{{\bm{x}}_{0}}\right\|\leq\kappa$. Then, we recover an
$\epsilon$-accurate solution as long as the error of DDPM generator is of the
order ${\mathcal{O}}\left(\epsilon\right)$.

####  3.2.4 Resampling vs Slowing Down Diffusion

We now generalize to a multi-state model. In 3.2.2 and 3.2.3, we discussed
that resampling plays a vital role in harmonizing the missing parts of an
image. We ask whether a similar result is achievable by slowing down the
diffusion process over $(T+1)$ states without resampling, instead of $2$
states with many resampling steps. To this end, we first derive a closed-form
solution for generative modeling with $(T+1)$ diffusion states in Theorem 8.
Then, we provide a simplified closed-form solution in Corollary 10 assuming
that the noise in the forward SDE and the reverse SDE are Independent and
Identically Distributed (IID). Using this solution, we show in Theorem 12 that
inpainting without resampling at intermediate states yields inferior results
compared to less states with many resampling steps.

######  Theorem 8.

Suppose Assumption 1 holds. For a deterministic variance schedule
$\\{\beta_{t}\\}_{t\geq 0}$, let $\alpha_{t}=1-\beta_{t}$,
$\bar{\alpha}_{t}=\prod_{s=0}^{t}\alpha_{s}$,
$\gamma\coloneqq\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\frac{1}{\left(1-\bar{\alpha}_{t}\right)}\right)\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\sqrt{\alpha}_{t}\left(1-\bar{\alpha}_{t-1}\right)\right)$
and
$\nu\coloneqq\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\frac{1}{\left(1-\bar{\alpha}_{t}\right)}\right)\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\bar{\alpha}_{t-1}\sqrt{\alpha_{t}}\right)-\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\frac{\bar{\alpha}_{t}}{\left(1-\bar{\alpha}_{t}\right)}\right)\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\sqrt{\alpha}_{t}\right)$,
where $t\sim Uniform\\{1,\cdots,T\\}$. Let us denote
${\bm{\theta}}^{*}=\arg_{\min}{\mathcal{L}}\left({\bm{\theta}}\right)$, where
${\mathcal{L}}\left({\bm{\theta}}\right)$ is defined as:

$\mathop{\mathbb{E}}_{\begin{subarray}{c}\overrightarrow{{\bm{x}}_{0}},\mathbf{\epsilon},t\end{subarray}}\left[\left\|\tilde{\mu}_{t}\left(\overrightarrow{{\bm{x}}_{t}},\overrightarrow{{\bm{x}}_{0}}\right)-\mu_{\theta}\left(\overrightarrow{{\bm{x}}_{t}},t\right)\right\|^{2}\right].$

For diffusion over $(T+1)$ states, if
$\mu_{\theta}\left(\overleftarrow{{\bm{x}}_{t}},t\right)={\bm{\theta}}\begin{bmatrix}\overleftarrow{{\bm{x}}_{t}}\\\
t\end{bmatrix}$, where ${\bm{\theta}}\in\mathbb{R}^{d\times(d+1)}$, then the
closed-form solution of DDPM Algorithm 1 is given by
${\bm{\theta}}^{*}=\left[\nu{\bm{A}}{\bm{A}}^{T}+\gamma{\bm{I}}_{d},\mathbf{0}\right]$.

###### Proof.

The proof is included in Appendix A.5. ∎

In general, $\\{\beta_{t}\\}_{t\geq 0}$ captures a list of variances ranging
from $\beta_{\min}$ to $\beta_{\max}$, such as a cosine or a linear schedule
[HJA20]. To capture the general framework, we provide $\nu$ and $\gamma$ as
explicit functions of $\beta_{t}$. One may wish to compute these values
exactly depending on the choice of variance schedule. For instance, in case of
diffusion over 2 states with $\beta_{0}=0$ and $\beta_{1}=\beta$, we have
$\nu=\left(\frac{1}{\left(1-\bar{\alpha}_{1}\right)}\right)\left(\bar{\alpha}_{0}\sqrt{\alpha_{1}}\right)-\left(\frac{\bar{\alpha}_{1}}{\left(1-\bar{\alpha}_{1}\right)}\right)\left(\sqrt{\alpha}_{1}\right)=\sqrt{1-\beta}$
and
$\gamma=\left(\frac{1}{\left(1-\bar{\alpha}_{1}\right)}\right)\left(\sqrt{\alpha}_{t}\left(1-\bar{\alpha}_{0}\right)\right)=0$.
This matches with the closed-form solution derived in Theorem 3.

For the sake of analysis, we present a simplified closed-form solution in
Corollary 10 under Assumption 9.

######  Assumption 9.

Let us denote by $\overrightarrow{\mathbf{\epsilon}}$ the noise used in the
dispersion of the forward SDE, i.e.,
$\overrightarrow{{\bm{x}}_{t-1}}=\sqrt{\bar{\alpha}_{t-1}}\overrightarrow{{\bm{x}}_{0}}+\sqrt{1-\bar{\alpha}_{t-1}}\overrightarrow{\mathbf{\epsilon}}$,
and $\overleftarrow{\mathbf{\epsilon}}$ in the reverse SDE, i.e.,
$\overleftarrow{{\bm{x}}_{t-1}}=\mu_{{\bm{\theta}}}\left(\overleftarrow{{\bm{x}}_{t}},t\right)+\sqrt{\beta_{t}}\overleftarrow{\mathbf{\epsilon}}$,
where
$\overleftarrow{{\bm{x}}_{t}}=\sqrt{\bar{\alpha}_{t}}\overrightarrow{{\bm{x}}_{0}}+\sqrt{1-\bar{\alpha}_{t}}\overleftarrow{\mathbf{\epsilon}}$.
Let $\overrightarrow{\mathbf{\epsilon}}$ and
$\overleftarrow{\mathbf{\epsilon}}$ be IID Gaussian random variables drawn
from ${\mathcal{N}}\left(\mathbf{0},{\bm{I}}_{d}\right)$.

######  Corollary 10.

Suppose Assumption 9 holds. For
$\bar{\nu}\coloneqq\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\frac{1}{\left(1-\bar{\alpha}_{t}\right)}\right)\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\bar{\alpha}_{t-1}\sqrt{\alpha_{t}}\right)-\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\frac{\bar{\alpha}_{t}}{\left(1-\bar{\alpha}_{t}\right)}\right)\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\bar{\alpha}_{t-1}\sqrt{\alpha_{t}}\right)$
in the same setting as Theorem 8, the closed-form solution of DDPM Algorithm 1
becomes
${\bm{\theta}}^{*}=\left[\bar{\nu}{\bm{A}}{\bm{A}}^{T},\mathbf{0}\right]$.

###### Proof.

The proof is included in Appendix A.6. ∎

For simplicity, let us use this closed-form solution for inpainting at
intermediate states. We prove in Theorem 11 that for diffusion over $(T+1)$
states, the generative model using DDPM learns to sample from the data
manifold.

######  Theorem 11.

For the coefficients of drift
$\omega_{t}=\omega=\frac{1}{\bar{\nu}\sqrt{2}^{1/T}}$ and dispersion
$\xi_{t}=\sqrt{\frac{2^{(t-1)/T}}{2\beta_{t}(T-1)}}$, if we choose the
Gaussian transition kernel using
${\bm{\theta}}^{*}=\left[\bar{\nu}{\bm{A}}{\bm{A}}^{T},\mathbf{0}\right]$ in
the reverse Markov process (3), then

$\displaystyle
TV\left(p_{{\bm{\theta}}^{*}}\left(\overleftarrow{{\bm{x}}_{0}}\right),q\left(\overrightarrow{{\bm{x}}_{0}}\right)\right)=0.$

###### Proof.

The proof is included in Appendix A.7. ∎

Since the total variation distance is zero, it ensures that DDPM learns the
target distribution $q\left(\overrightarrow{{\bm{x}}_{0}}\right)$. Now, let us
see whether inpainting known portions of $\overrightarrow{{\bm{x}}_{0}}$ at
intermediate states would recover the true sample. In this regard, an
important distinction from RePaint+ is that there is no resampling at
intermediate states.

######  Theorem 12.

For the coefficients of drift
$\omega_{t}=\omega=\frac{1}{\bar{\nu}\sqrt{2}^{1/T}}$ and dispersion
$\xi_{t}=\sqrt{\frac{2^{(t-1)/T}}{2\beta_{t}(T-1)}}$, following the reverse
Markov process (3) with inpainting at intermediate states, generates a sample
satisfying:

$\displaystyle\overleftarrow{{\bm{x}}_{0}}$
$\displaystyle=\frac{1}{\sqrt{2}}\left({\bm{A}}{\bm{A}}^{T}{\bm{D}}({\bm{m}})\right)^{T-1}{\bm{A}}{\bm{A}}^{T}\overleftarrow{{\bm{x}}_{T}}$
$\displaystyle+\left(\sum_{t=1}^{T-1}\left(\frac{1}{\sqrt{2}}\right)^{t/T}\sqrt{\bar{\alpha}_{t}}\left({\bm{A}}{\bm{A}}^{T}{\bm{D}}({\bm{m}})\right)^{t-1}\right){\bm{A}}{\bm{A}}^{T}{\bm{D}}(\bm{1}-{\bm{m}})\overrightarrow{{\bm{x}}_{0}}+\overrightarrow{\bar{\Sigma}}\overrightarrow{\mathbf{\epsilon}}+\overleftarrow{\bar{\Sigma}}\overleftarrow{\mathbf{\epsilon}},$

where $\overrightarrow{\bar{\Sigma}}$ and $\overleftarrow{\bar{\Sigma}}$
denote the covariance of collective dispersion in the forward and the reverse
processes, respectively.

###### Proof.

The proof is included in Appendix A.8. ∎

We note that $\overleftarrow{{\bm{x}}_{0}}$ is not necessarily an image
residing on the manifold of the data generative distribution
$q\left(\overrightarrow{{\bm{x}}_{0}}\right)$. This is partly due to the
residual bias arising from the noise injected in the forward and the reverse
processes. Another cause of concern is the lack of harmonization between
generated portions and rest of the image. Notably, the bias is persistent even
if we consider diffusion over infinitely many states, justifying the
importance of resampling over slowing down the diffusion process. This
phenomenon is also observed in practice (5, [LDRYTV22]).

On the contrary, diffusion over two states with many resampling steps recovers
the original sample as discussed in 3.2.2. In Algorithm 3, we generalize this
notion to $(T+1)$ states, where the generated image is harmonized at each
intermediate state. Thus, we gradually wipe out the contribution of noise
before it gets accumulated through the reverse transitions, leading to a
smooth recovery.

##  4 Conclusion

We studied mathematical underpinnings of sample recovery in image inpainting
using denoising diffusion probabilistic models. We presented a universal mask
principle to prove that diffusion based inpainting easily adapts to unseen
masks without retraining. Our theoretical analysis of a popular inpainting
algorithm in a linear setting revealed a previously not understood bias that
hampered perfect recovery. With realignment of drift and dispersion in the
reverse process, we proposed an algorithm called Repaint+ that eliminated this
bias, leading to perfect recovery in the limit. Further, We proved that
RePaint+ converged linearly to the $\epsilon$-neighborhood of the true
underlying sample. In the future, it would be interesting to study the
convergence properties in a more general non-linear manifold setting.

## References

  * [BHSN10] Waheed U. Bajwa, Jarvis Haupt, Akbar M. Sayeed and Robert Nowak  “Compressed Channel Sensing: A New Approach to Estimating Sparse Multipath Channels”  In _Proceedings of the IEEE_ 98.6, 2010, pp. 1058–1076  DOI: 10.1109/JPROC.2010.2042415
  * [BJPD17] Ashish Bora, Ajil Jalal, Eric Price and Alexandros G Dimakis  “Compressed sensing using generative models”  In _International Conference on Machine Learning_ , 2017, pp. 537–546  PMLR
  * [Bor22] Valentin De Bortoli  “Convergence of denoising diffusion models under the manifold hypothesis”  In _Transactions on Machine Learning Research_ , 2022  URL: https://openreview.net/forum?id=MhK5aXo3gB
  * [CLL22] Hongrui Chen, Holden Lee and Jianfeng Lu  “Improved Analysis of Score-based Generative Modeling: User-Friendly Bounds under Minimal Smoothness Assumptions”  In _arXiv preprint arXiv:2211.01916_ , 2022
  * [CCLLSZ22] 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 _arXiv preprint arXiv:2209.11215_ , 2022
  * [CKJGY21] Jooyoung Choi, Sungwon Kim, Yonghyun Jeong, Youngjune Gwon and Sungroh Yoon  “ILVR: Conditioning Method for Denoising Diffusion Probabilistic Models”  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , 2021, pp. 14367–14376
  * [DDDD22] Giannis Daras, Yuval Dagan, Alex Dimakis and Constantinos Daskalakis  “Score-Guided Intermediate Level Optimization: Fast Langevin Mixing for Inverse Problems”  In _International Conference on Machine Learning_ , 2022, pp. 4722–4753  PMLR
  * [DGE18] Manik Dhar, Aditya Grover and Stefano Ermon  “Modeling sparse deviations for compressed sensing using generative models”  In _International Conference on Machine Learning_ , 2018, pp. 1214–1223  PMLR
  * [DN21] Prafulla Dhariwal and Alexander Nichol  “Diffusion models beat gans on image synthesis”  In _Advances in Neural Information Processing Systems_ 34, 2021, pp. 8780–8794
  * [DSB17] Laurent Dinh, Jascha Sohl-Dickstein and Samy Bengio  “Density estimation using Real NVP”  In _International Conference on Learning Representations_ , 2017  URL: https://openreview.net/forum?id=HkpbnH9lx
  * [FMN16] Charles Fefferman, Sanjoy Mitter and Hariharan Narayanan  “Testing the manifold hypothesis”  In _Journal of the American Mathematical Society_ 29.4, 2016, pp. 983–1049
  * [FN03] M.A.T. Figueiredo and R.D. Nowak  “An EM algorithm for wavelet-based image restoration”  In _IEEE Transactions on Image Processing_ 12.8, 2003, pp. 906–916  DOI: 10.1109/TIP.2003.814255
  * [FNW07] MÁrio A.. Figueiredo, Robert D. Nowak and Stephen J. Wright  “Gradient Projection for Sparse Reconstruction: Application to Compressed Sensing and Other Inverse Problems”  In _IEEE Journal of Selected Topics in Signal Processing_ 1.4, 2007, pp. 586–597  DOI: 10.1109/JSTSP.2007.910281
  * [GPMXWOCB14] Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville and Yoshua Bengio  “Generative Adversarial Nets”  In _Advances in Neural Information Processing Systems_ 27, 2014
  * [Hag89] William W Hager  “Updating the inverse of a matrix”  In _SIAM review_ 31.2 SIAM, 1989, pp. 221–239
  * [HJA20] Jonathan Ho, Ajay Jain and Pieter Abbeel  “Denoising diffusion probabilistic models”  In _Advances in Neural Information Processing Systems_ 33, 2020, pp. 6840–6851
  * [HD05] Aapo Hyvärinen and Peter Dayan  “Estimation of non-normalized statistical models by score matching.”  In _Journal of Machine Learning Research_ 6.4, 2005
  * [JADPDT21] Ajil Jalal, Marius Arvinte, Giannis Daras, Eric Price, Alexandros G Dimakis and Jon Tamir  “Robust compressed sensing mri with deep generative priors”  In _Advances in Neural Information Processing Systems_ 34, 2021, pp. 14938–14954
  * [JLDC20] Ajil Jalal, Liu Liu, Alexandros G Dimakis and Constantine Caramanis  “Robust compressed sensing using generative models”  In _Advances in Neural Information Processing Systems_ 33, 2020, pp. 713–727
  * [Jar97] Christopher Jarzynski  “Equilibrium free-energy differences from nonequilibrium measurements: A master-equation approach”  In _Physical Review E_ 56.5 APS, 1997, pp. 5018
  * [KAAL22] Tero Karras, Miika Aittala, Timo Aila and Samuli Laine  “Elucidating the Design Space of Diffusion-Based Generative Models”  In _Advances in Neural Information Processing Systems_ , 2022  URL: https://openreview.net/forum?id=k7FuTOWMOc7
  * [KEES22] Bahjat Kawar, Michael Elad, Stefano Ermon and Jiaming Song  “Denoising Diffusion Restoration Models”  In _Advances in Neural Information Processing Systems_ , 2022
  * [KSSKM22] Dongjun Kim, Seungjae Shin, Kyungwoo Song, Wanmo Kang and Il-Chul Moon  “Soft Truncation: A Universal Training Technique of Score-based Diffusion Model for High Precision Score Estimation”  In _Proceedings of the 39th International Conference on Machine Learning_ 162, Proceedings of Machine Learning Research  PMLR, 2022, pp. 11201–11228  URL: https://proceedings.mlr.press/v162/kim22i.html
  * [KW13] Diederik P Kingma and Max Welling  “Auto-encoding variational bayes”  In _arXiv preprint arXiv:1312.6114_ , 2013
  * [KPHZC21] Zhifeng Kong, Wei Ping, Jiaji Huang, Kexin Zhao and Bryan Catanzaro  “DiffWave: A Versatile Diffusion Model for Audio Synthesis”  In _International Conference on Learning Representations_ , 2021  URL: https://openreview.net/forum?id=a-xFK8Ymz5J
  * [LLT22] Holden Lee, Jianfeng Lu and Yixin Tan  “Convergence for score-based generative modeling with polynomial complexity”  In _arXiv preprint arXiv:2206.06227_ , 2022
  * [LLT22a] Holden Lee, Jianfeng Lu and Yixin Tan  “Convergence of score-based generative modeling for general data distributions”  In _arXiv preprint arXiv:2209.12381_ , 2022
  * [LTGLH22] Xiang Lisa Li, John Thickstun, Ishaan Gulrajani, Percy Liang and Tatsunori Hashimoto  “Diffusion-LM Improves Controllable Text Generation”  In _Advances in Neural Information Processing Systems_ , 2022  URL: https://openreview.net/forum?id=3s9IrEsjLyk
  * [LDRYTV22] Andreas Lugmayr, Martin Danelljan, Andres Romero, Fisher Yu, Radu Timofte and Luc Van Gool  “Repaint: Inpainting using denoising diffusion probabilistic models”  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2022, pp. 11461–11471
  * [PLXL21] Jialun Peng, Dong Liu, Songcen Xu and Houqiang Li  “Generating diverse structure for image inpainting with hierarchical VQ-VAE”  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2021, pp. 10775–10784
  * [Pid22] Jakiw Pidstrigach  “Score-Based Generative Models Detect Manifolds”  In _Advances in Neural Information Processing Systems_ , 2022  URL: https://openreview.net/forum?id=AiNrnIrDfD9
  * [RDNCC22] Aditya Ramesh, Prafulla Dhariwal, Alex Nichol, Casey Chu and Mark Chen  “Hierarchical text-conditional image generation with clip latents”  In _arXiv preprint arXiv:2204.06125_ , 2022
  * [RPGGVRCS21] Aditya Ramesh, Mikhail Pavlov, Gabriel Goh, Scott Gray, Chelsea Voss, Alec Radford, Mark Chen and Ilya Sutskever  “Zero-shot text-to-image generation”  In _International Conference on Machine Learning_ , 2021, pp. 8821–8831  PMLR
  * [RVV19] Ali Razavi, Aaron Van den Oord and Oriol Vinyals  “Generating diverse high-fidelity images with vq-vae-2”  In _Advances in neural information processing systems_ 32, 2019
  * [RBLEO22] Robin Rombach, Andreas Blattmann, Dominik Lorenz, Patrick Esser and Björn Ommer  “High-Resolution Image Synthesis With Latent Diffusion Models”  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2022, pp. 10684–10695
  * [SCCLHSFN22] Chitwan Saharia, William Chan, Huiwen Chang, Chris Lee, Jonathan Ho, Tim Salimans, David Fleet and Mohammad Norouzi  “Palette: Image-to-image diffusion models”  In _ACM SIGGRAPH 2022 Conference Proceedings_ , 2022, pp. 1–10
  * [SCSLWDGAML+22] Chitwan Saharia, William Chan, Saurabh Saxena, Lala Li, Jay Whang, Emily Denton, Seyed Kamyar Seyed Ghasemipour, Burcu Karagol Ayan, S Sara Mahdavi and Rapha Gontijo Lopes  “Photorealistic Text-to-Image Diffusion Models with Deep Language Understanding”  In _arXiv preprint arXiv:2205.11487_ , 2022
  * [SWMG15] Jascha Sohl-Dickstein, Eric Weiss, Niru Maheswaranathan and Surya Ganguli  “Deep unsupervised learning using nonequilibrium thermodynamics”  In _International Conference on Machine Learning_ , 2015  PMLR
  * [SE19] Yang Song and Stefano Ermon  “Generative modeling by estimating gradients of the data distribution”  In _Advances in Neural Information Processing Systems_ 32, 2019
  * [SSXE21] Yang Song, Liyue Shen, Lei Xing and Stefano Ermon  “Solving Inverse Problems in Medical Imaging with Score-Based Generative Models”  In _International Conference on Learning Representations_ , 2021
  * [SSKKEP20] Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon and Ben Poole  “Score-Based Generative Modeling through Stochastic Differential Equations”  In _International Conference on Learning Representations_ , 2020
  * [TSSE21] Yusuke Tashiro, Jiaming Song, Yang Song and Stefano Ermon  “CSDI: Conditional score-based diffusion models for probabilistic time series imputation”  In _Advances in Neural Information Processing Systems_ 34, 2021, pp. 24804–24816
  * [TSL00] Joshua B Tenenbaum, Vin de Silva and John C Langford  “A global geometric framework for nonlinear dimensionality reduction”  In _science_ 290.5500 American Association for the Advancement of Science, 2000, pp. 2319–2323
  * [Tib96] Robert Tibshirani  “Regression shrinkage and selection via the lasso”  In _Journal of the Royal Statistical Society: Series B (Methodological)_ 58.1 Wiley Online Library, 1996, pp. 267–288
  * [Vin11] Pascal Vincent  “A connection between score matching and denoising autoencoders”  In _Neural computation_ 23.7 MIT Press, 2011, pp. 1661–1674
  * [WYWGLDQC18] Xintao Wang, Ke Yu, Shixiang Wu, Jinjin Gu, Yihao Liu, Chao Dong, Yu Qiao and Chen Change Loy  “Esrgan: Enhanced super-resolution generative adversarial networks”  In _Proceedings of the European conference on computer vision (ECCV) workshops_ , 2018
  * [WDTSDM22] Jay Whang, Mauricio Delbracio, Hossein Talebi, Chitwan Saharia, Alexandros G. Dimakis and Peyman Milanfar  “Deblurring via Stochastic Refinement”  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2022, pp. 16293–16303
  * [WRL19] Yan Wu, Mihaela Rosca and Timothy Lillicrap  “Deep compressed sensing”  In _International Conference on Machine Learning_ , 2019, pp. 6850–6860  PMLR
  * [YLYSLH18] Jiahui Yu, Zhe Lin, Jimei Yang, Xiaohui Shen, Xin Lu and Thomas S Huang  “Generative image inpainting with contextual attention”  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2018, pp. 5505–5514

##  Appendix A Technical Proofs

###  A.1 Proof of Theorem 3

For diffusion over two states, the posterior mean
$\tilde{\mu}_{t}\left(\overrightarrow{{\bm{x}}_{t}},\overrightarrow{{\bm{x}}_{0}}\right)\coloneqq\frac{\sqrt{\bar{\alpha}_{t-1}}\beta_{t}}{1-\bar{a}_{t}}\overrightarrow{{\bm{x}}_{0}}+\frac{\sqrt{\alpha}_{t}\left(1-\bar{\alpha}_{t-1}\right)}{1-\bar{\alpha}_{t}}\overrightarrow{{\bm{x}}_{t}}$
simplifies to
$\tilde{\mu}_{1}\left(\overrightarrow{{\bm{x}}_{1}},\overrightarrow{{\bm{x}}_{0}}\right)=\overrightarrow{{\bm{x}}_{0}}$.
Thus, the training loss becomes:

$\displaystyle\min_{{\bm{\theta}}}\mathbb{E}_{\overrightarrow{{\bm{x}}_{0}},\overrightarrow{\mathbf{\epsilon}}}\left[\left\|\tilde{\mu}_{1}(\overrightarrow{{\bm{x}}_{1}}\left(\overrightarrow{{\bm{x}}_{0}},\overrightarrow{\mathbf{\epsilon}}),\overrightarrow{{\bm{x}}_{0}}\right)-\mu_{\theta}\left(\overrightarrow{{\bm{x}}_{1}}\left(\overrightarrow{{\bm{x}}_{0}},\overrightarrow{\mathbf{\epsilon}}\right)\right)\right\|^{2}\right]=$
$\displaystyle\mathbb{E}_{\overrightarrow{{\bm{x}}_{0}},\overrightarrow{\mathbf{\epsilon}}}\left[\left\|\overrightarrow{{\bm{x}}_{0}}-\mu_{\bm{\theta}}\left(\overrightarrow{{\bm{x}}_{1}}\right)\right\|^{2}\right]=\mathbb{E}_{\overrightarrow{{\bm{x}}_{0}},\overrightarrow{\mathbf{\epsilon}}}\left[\left\|\overrightarrow{{\bm{x}}_{0}}-{\bm{\theta}}\overrightarrow{{\bm{x}}_{1}}\right\|^{2}\right]=$
$\displaystyle\mathop{\mathbb{E}}_{\begin{subarray}{c}\overrightarrow{{\bm{x}}_{0}}\sim
q\\\
\overrightarrow{\mathbf{\epsilon}}\sim\mathcal{N}\left(\mathbf{0},{\bm{I}}_{d}\right)\end{subarray}}\left[\left\|\overrightarrow{{\bm{x}}_{0}}-{\bm{\theta}}\left({\bm{x}}_{0}\sqrt{1-\beta}+\overrightarrow{\mathbf{\epsilon}}\sqrt{\beta}\right)\right\|^{2}\right]=$
$\displaystyle\mathop{\mathbb{E}}_{\begin{subarray}{c}\overrightarrow{{\bm{x}}_{0}}\sim
q\\\
\overrightarrow{\mathbf{\epsilon}}\sim\mathcal{N}\left(\mathbf{0},{\bm{I}}_{d}\right)\end{subarray}}\left[\sum_{i=1}^{d}\left(\overrightarrow{{\bm{x}}_{0,i}}-{\bm{\theta}}_{i}^{T}\left(\overrightarrow{{\bm{x}}_{0}}\sqrt{1-\beta}+\overrightarrow{\mathbf{\epsilon}}\sqrt{\beta}\right)\right)^{2}\right],$
(8)

where ${\bm{\theta}}_{i}^{T}$ denotes the $i^{th}$ row of ${\bm{\theta}}$.
Solving the MMSE problem (8) for ${\bm{\theta}}_{i}^{*}$ yields

$\displaystyle{\bm{\theta}}_{i}^{*}=\mathop{\mathbb{E}}_{\begin{subarray}{c}\overrightarrow{{\bm{x}}_{0}},\overrightarrow{\mathbf{\epsilon}}\end{subarray}}\left[\left(\overrightarrow{{\bm{x}}_{0}}\sqrt{1-\beta}+\overrightarrow{\mathbf{\epsilon}}\sqrt{\beta}\right)\left(\overrightarrow{{\bm{x}}_{0}}\sqrt{1-\beta}+\overrightarrow{\mathbf{\epsilon}}\sqrt{\beta}\right)^{T}\right]^{-1}\mathbb{E}_{\overrightarrow{{\bm{x}}_{0}},\overrightarrow{\mathbf{\epsilon}}}\left[\overrightarrow{{\bm{x}}_{0,i}}\left(\overrightarrow{{\bm{x}}_{0}}\sqrt{1-\beta}+\overrightarrow{\mathbf{\epsilon}}\sqrt{\beta}\right)\right]=$
$\displaystyle\mathop{\mathbb{E}}_{\begin{subarray}{c}\overrightarrow{{\bm{z}}_{0}},\overrightarrow{\mathbf{\epsilon}}\end{subarray}}\left[\left(\mathcal{A}{\bm{z}}_{0}\sqrt{1-\beta}+\overrightarrow{\mathbf{\epsilon}}\sqrt{\beta}\right)\left(\mathcal{A}{\bm{z}}_{0}\sqrt{1-\beta}+\overrightarrow{\mathbf{\epsilon}}\sqrt{\beta}\right)^{T}\right]^{-1}\mathbb{E}_{{\bm{z}}_{0},\overrightarrow{\mathbf{\epsilon}}}\left[\left({\bm{a}}_{i}^{T}{\bm{z}}_{0}\right)\left(\mathcal{A}{\bm{z}}_{0}\sqrt{1-\beta}+\overrightarrow{\mathbf{\epsilon}}\sqrt{\beta}\right)\right]=$
$\displaystyle\Big{[}(1-\beta)\mathcal{A}\mathbb{E}_{{\bm{z}}_{0}}\left[{\bm{z}}_{0}{\bm{z}}_{0}^{T}\right]\mathcal{A}^{T}+\beta\mathbb{E}_{\overrightarrow{\mathbf{\epsilon}}}\left[\overrightarrow{\mathbf{\epsilon}}\overrightarrow{\mathbf{\epsilon}}^{T}\right]\Big{]}^{-1}\left[\left(\sqrt{1-\beta}\right)\mathcal{A}\mathbb{E}_{\overrightarrow{{\bm{x}}_{0}}}\left[\overrightarrow{{\bm{x}}_{0}}\overrightarrow{{\bm{x}}_{0}}^{T}\right]{\bm{a}}_{i}\right].$
(9)

Using the fact that ${\bm{z}}_{0}$ and $\mathbf{\epsilon}$ are independent
Gaussian random variables with unit covariance, (9) gives

$\displaystyle{\bm{\theta}}_{i}^{*}=\left[(1-\beta)\mathcal{A}\mathcal{A}^{T}+\beta{\bm{I}}_{d}\right]^{-1}\left[\left(\sqrt{1-\beta}\right)\mathcal{A}{\bm{a}}_{i}\right]=$
$\displaystyle\frac{\sqrt{1-\beta}}{\beta}\left[{\bm{I}}_{d}+\left(\sqrt{\frac{1-\beta}{\beta}}\mathcal{A}\right)\left(\sqrt{\frac{1-\beta}{\beta}}\mathcal{A}\right)^{T}\right]^{-1}\mathcal{A}{\bm{a}}_{i}.$

Next, we use Woodbury matrix identity [Hag89] to get

$\displaystyle{\bm{\theta}}_{i}^{*}=$
$\displaystyle\frac{\sqrt{1-\beta}}{\beta}\Big{[}{\bm{I}}_{d}-\left(\sqrt{\frac{1-\beta}{\beta}}\mathcal{A}\right)\biggl{(}{\bm{I}}_{k}+\left(\sqrt{\frac{1-\beta}{\beta}}\mathcal{A}\right)^{T}\left(\sqrt{\frac{1-\beta}{\beta}}\mathcal{A}\right)\biggr{)}^{-1}\left(\sqrt{\frac{1-\beta}{\beta}}\mathcal{A}\right)^{T}\Big{]}\mathcal{A}{\bm{a}}_{i}=$
$\displaystyle\frac{\sqrt{1-\beta}}{\beta}\Big{[}{\bm{I}}_{d}-\left(\sqrt{\frac{1-\beta}{\beta}}{\bm{A}}\right)\biggl{(}{\bm{I}}_{k}\left(\frac{1-\beta}{\beta}\right)\mathcal{A}^{T}\mathcal{A}\biggr{)}^{-1}\left(\sqrt{\frac{1-\beta}{\beta}}\mathcal{A}\right)^{T}\Big{]}\mathcal{A}{\bm{a}}_{i}=$
$\displaystyle\frac{\sqrt{1-\beta}}{\beta}\left[{\bm{I}}_{d}-\left(\sqrt{\frac{1-\beta}{\beta}}{\bm{A}}\right)\beta{\bm{I}}_{k}\left(\sqrt{\frac{1-\beta}{\beta}}\mathcal{A}\right)^{T}\right]\mathcal{A}{\bm{a}}_{i}=$
$\displaystyle\frac{\sqrt{1-\beta}}{\beta}\left[{\bm{I}}_{d}-\left(1-\beta\right){\bm{A}}{\bm{A}}^{T}\right]\mathcal{A}{\bm{a}}_{i}=$
$\displaystyle\frac{\sqrt{1-\beta}}{\beta}\left[\mathcal{A}{\bm{a}}_{i}-\left(1-\beta\right){\bm{A}}{\bm{a}}_{i}\right]=\sqrt{1-\beta}{\bm{A}}{\bm{a}}_{i}.$

Stacking all the rows together, we get the optimal solution of the MMSE
problem (8), i.e., ${\bm{\theta}}^{*}=\sqrt{1-\beta}{\bm{A}}{\bm{A}}^{T}$.
This optimal solution characterizes the Gaussian transition kernel of the
reverse SDE (3). $\square$

###  A.2 Proof of Theorem 4

After $r$ resampling steps, we get the following output from the reverse
process:

$\displaystyle\overleftarrow{{\bm{x}}_{0}^{r}}=$
$\displaystyle{\bm{A}}{\bm{A}}^{T}\left({\bm{D}}\left({\bm{m}}\right)\overleftarrow{{\bm{x}}_{0}^{r-1}}+{\bm{D}}\left(\bm{1}-{\bm{m}}\right)\overrightarrow{{\bm{x}}_{0}}\right)=$
$\displaystyle{\bm{A}}{\bm{A}}^{T}{\bm{D}}\left({\bm{m}}\right)\overleftarrow{{\bm{x}}_{0}^{r-1}}+{\bm{A}}{\bm{A}}^{T}{\bm{D}}\left(\bm{1}-{\bm{m}}\right)\overrightarrow{{\bm{x}}_{0}}=$
$\displaystyle{\bm{A}}{\bm{A}}^{T}{\bm{D}}\left({\bm{m}}\right)\biggl{(}{\bm{A}}{\bm{A}}^{T}{\bm{D}}\left({\bm{m}}\right)\overleftarrow{{\bm{x}}_{0}^{r-2}}+{\bm{A}}{\bm{A}}^{T}{\bm{D}}\left(\bm{1}-{\bm{m}}\right)\overrightarrow{{\bm{x}}_{0}}\biggr{)}+{\bm{A}}{\bm{A}}^{T}{\bm{D}}\left(\bm{1}-{\bm{m}}\right)\overrightarrow{{\bm{x}}_{0}}=$
$\displaystyle\left({\bm{A}}{\bm{A}}^{T}{\bm{D}}\left({\bm{m}}\right)\right)^{r}\overleftarrow{{\bm{x}}_{0}}+\left({\bm{I}}+\sum_{i=1}^{r-1}\left({\bm{A}}{\bm{A}}^{T}{\bm{D}}({\bm{m}})\right)^{i}\right)\biggl{(}{\bm{A}}{\bm{A}}^{T}{\bm{D}}(\bm{1}-{\bm{m}}){\bm{A}}{\bm{z}}_{0}\biggr{)}.$
(10)

Using the truncation of Neumann series, (10) simplifies to

$\displaystyle\overleftarrow{{\bm{x}}_{0}^{r}}=\left({\bm{A}}{\bm{A}}^{T}{\bm{D}}\left({\bm{m}}\right)\right)^{r}\overleftarrow{{\bm{x}}_{0}}+\left({\bm{I}}-\left({\bm{A}}{\bm{A}}^{T}{\bm{D}}({\bm{m}})\right)^{r}\right)\left({\bm{I}}-{\bm{A}}{\bm{A}}^{T}{\bm{D}}({\bm{m}})\right)^{-1}{\bm{A}}{\bm{A}}^{T}{\bm{D}}(\bm{1}-{\bm{m}}){\bm{A}}{\bm{z}}_{0}=$
$\displaystyle\left({\bm{A}}{\bm{A}}^{T}{\bm{D}}\left({\bm{m}}\right)\right)^{r}\overleftarrow{{\bm{x}}_{0}}+\left({\bm{I}}-\left({\bm{A}}{\bm{A}}^{T}{\bm{D}}({\bm{m}})\right)^{r}\right)\left({\bm{I}}-{\bm{A}}{\bm{A}}^{T}{\bm{D}}({\bm{m}})\right)^{-1}{\bm{A}}{\bm{A}}^{T}\left({\bm{I}}-{\bm{D}}({\bm{m}})\right){\bm{A}}{\bm{z}}_{0}=$
$\displaystyle\left({\bm{A}}{\bm{A}}^{T}{\bm{D}}\left({\bm{m}}\right)\right)^{r}\overleftarrow{{\bm{x}}_{0}}+\left({\bm{I}}-\left({\bm{A}}{\bm{A}}^{T}{\bm{D}}({\bm{m}})\right)^{r}\right)\left({\bm{I}}-{\bm{A}}{\bm{A}}^{T}{\bm{D}}({\bm{m}})\right)^{-1}\left({\bm{A}}{\bm{A}}^{T}{\bm{A}}{\bm{z}}_{0}-{\bm{A}}{\bm{A}}^{T}{\bm{D}}({\bm{m}}){\bm{A}}{\bm{z}}_{0}\right)=$
$\displaystyle\left({\bm{A}}{\bm{A}}^{T}{\bm{D}}\left({\bm{m}}\right)\right)^{r}\overleftarrow{{\bm{x}}_{0}}+\left({\bm{I}}-\left({\bm{A}}{\bm{A}}^{T}{\bm{D}}({\bm{m}})\right)^{r}\right)\left({\bm{I}}-{\bm{A}}{\bm{A}}^{T}{\bm{D}}({\bm{m}})\right)^{-1}\left({\bm{A}}{\bm{z}}_{0}-{\bm{A}}{\bm{A}}^{T}{\bm{D}}({\bm{m}}){\bm{A}}{\bm{z}}_{0}\right)=$
$\displaystyle\left({\bm{A}}{\bm{A}}^{T}{\bm{D}}\left({\bm{m}}\right)\right)^{r}{\bm{A}}{\bm{A}}^{T}\overleftarrow{{\bm{x}}_{1}}+\left({\bm{I}}-\left({\bm{A}}{\bm{A}}^{T}{\bm{D}}({\bm{m}})\right)^{r}\right)\left({\bm{I}}-{\bm{A}}{\bm{A}}^{T}{\bm{D}}({\bm{m}})\right)^{-1}\left({\bm{I}}-{\bm{A}}{\bm{A}}^{T}{\bm{D}}({\bm{m}})\right){\bm{A}}{\bm{z}}_{0}=$
$\displaystyle\left({\bm{A}}{\bm{A}}^{T}{\bm{D}}\left({\bm{m}}\right)\right)^{r}\left({\bm{A}}{\bm{A}}^{T}\overleftarrow{{\bm{x}}_{1}}-{\bm{A}}{\bm{z}}_{0}\right)+{\bm{A}}{\bm{z}}_{0}.$
(11)

Next, we use the Singular Value Decomposition (SVD) of
${\bm{A}}{\bm{A}}^{T}{\bm{D}}({\bm{m}})$. Let us denote the left singular
vectors of ${\bm{A}}{\bm{A}}^{T}{\bm{D}}({\bm{m}})$ by
${\bm{U}}\in\mathbb{R}^{d\times d}$, the eigen values by the diagonal matrix
$\mathbf{\Sigma}\in\mathbb{R}^{d\times d}$, and the right singular vectors by
${\bm{V}}\in\mathbb{R}^{d\times d}$. Therefore, we simplify the expression
(11) to

$\displaystyle\overleftarrow{{\bm{x}}_{0}^{r}}$
$\displaystyle={\bm{U}}\mathbf{\Sigma}^{r}{\bm{V}}\left({\bm{A}}{\bm{A}}^{T}\overleftarrow{{\bm{x}}_{1}}-{\bm{A}}{\bm{z}}_{0}\right)+{\bm{A}}{\bm{z}}_{0}.$

Since the largest eigen value of ${\bm{A}}{\bm{A}}^{T}{\bm{D}}({\bm{m}})$ is
stricly less than $1$, we get the following result:
$\lim_{r\rightarrow\infty}\overleftarrow{{\bm{x}}_{0}^{r}}={\bm{A}}{\bm{z}}_{0}$.
This answers the question of recoverability by diffusion based inpainting.
Indeed, Algorithm 2 recovers the exact sample given infinitely many resampling
steps. For most practical applications, it is sufficient to get an
$\mathbf{\epsilon}$-close solution. Now, let us derive the rate at which
Algorithm 2 converges to the $\mathbf{\epsilon}$-neighborhood of
${\bm{A}}{\bm{z}}_{0}$.

Subtracting the original sample
$\overrightarrow{{\bm{x}}_{0}}={\bm{A}}{\bm{z}}_{0}$ from both sides of (11)
and taking the Euclidean norm, we arrive at

$\displaystyle\left\|\overleftarrow{{\bm{x}}_{0}^{r}}-\overrightarrow{{\bm{x}}_{0}}\right\|=$
$\displaystyle\left\|\left({\bm{A}}{\bm{A}}^{T}{\bm{D}}\left({\bm{m}}\right)\right)^{r}\left({\bm{A}}{\bm{A}}^{T}\overleftarrow{{\bm{x}}_{1}}-{\bm{A}}{\bm{z}}_{0}\right)\right\|\stackrel{{\scriptstyle(i)}}{{\leq}}$
$\displaystyle\left\|{\bm{A}}{\bm{A}}^{T}{\bm{D}}({\bm{m}})\right\|^{r}\left\|{\bm{A}}{\bm{A}}^{T}\overleftarrow{{\bm{x}}_{1}}-{\bm{A}}{\bm{z}}_{0}\right\|=$
$\displaystyle\left\|{\bm{A}}{\bm{A}}^{T}{\bm{D}}({\bm{m}})\right\|^{r}\left\|{\bm{A}}{\bm{A}}^{T}\overleftarrow{{\bm{x}}_{1}}-{\bm{A}}{\bm{A}}^{T}{\bm{A}}{\bm{z}}_{0}\right\|\leq$
$\displaystyle\lambda_{\max}^{r}\frac{\left\|{\bm{\theta}}^{*}\right\|}{\sqrt{1-\beta}}\left\|\overleftarrow{{\bm{x}}_{1}}-\overrightarrow{{\bm{x}}_{0}}\right\|,$

where (i) uses Cauchy-Schwarz inequality and the spectral norm identity. By
making
$\lambda_{\max}^{r}\frac{\left\|{\bm{\theta}}^{*}\right\|}{\sqrt{1-\beta}}\left\|\overleftarrow{{\bm{x}}_{1}}-\overrightarrow{{\bm{x}}_{0}}\right\|\leq\epsilon$,
we ensure that the algorithm converges to an $\epsilon$-accurate solution. We
achieve this convergence by running at least
$\frac{\log\left(\frac{\left\|\theta^{*}\right\|\left\|\overleftarrow{x_{1}}-\overrightarrow{{\bm{x}}_{0}}\right\|}{\epsilon\sqrt{1-\beta}}\right)}{\log\left(\frac{1}{\lambda_{\max}}\right)}$
resampling steps, which implies a linear rate of convergence. $\square$

###  A.3 Proof of Theorem 6

Similar to the proof of Theorem 4, we begin with the following expression in
the case of noisy generative model:

$\displaystyle\overleftarrow{{\bm{x}}_{0}^{r}}=\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right)\left({\bm{D}}\left({\bm{m}}\right)\overleftarrow{{\bm{x}}_{0}^{r-1}}+{\bm{D}}\left(\bm{1}-{\bm{m}}\right)\overrightarrow{{\bm{x}}_{0}}\right)=$
$\displaystyle\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}\left({\bm{m}}\right)\overleftarrow{{\bm{x}}_{0}^{r-1}}+\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}\left(\bm{1}-{\bm{m}}\right)\overrightarrow{{\bm{x}}_{0}}=$
$\displaystyle\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}\left({\bm{m}}\right)\biggl{(}\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}\left({\bm{m}}\right)\overleftarrow{{\bm{x}}_{0}^{r-2}}+\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}\left(\bm{1}-{\bm{m}}\right)\overrightarrow{{\bm{x}}_{0}}\biggr{)}+\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}\left(\bm{1}-{\bm{m}}\right)\overrightarrow{{\bm{x}}_{0}}=$
$\displaystyle\left(\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}\left({\bm{m}}\right)\right)^{r}\overleftarrow{{\bm{x}}_{0}}+\left(\sum_{i=0}^{r-1}\left(\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}({\bm{m}})\right)^{i}\right)\biggl{(}\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}(\bm{1}-{\bm{m}}){\bm{A}}{\bm{z}}_{0}\biggr{)}=$
$\displaystyle\left(\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}\left({\bm{m}}\right)\right)^{r}\overleftarrow{{\bm{x}}_{0}}+\left({\bm{I}}-\left(\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}({\bm{m}})\right)^{r}\right)\left({\bm{I}}-\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}({\bm{m}})\right)^{-1}$
$\displaystyle\times\left({\bm{A}}{\bm{A}}^{T}{\bm{A}}{\bm{z}}_{0}+\mathbf{\delta}{\bm{A}}{\bm{z}}_{0}-\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}({\bm{m}}){\bm{A}}{\bm{z}}_{0}\right)=$
$\displaystyle\left(\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}\left({\bm{m}}\right)\right)^{r}\overleftarrow{{\bm{x}}_{0}}+\left({\bm{I}}-\left(\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}({\bm{m}})\right)^{r}\right)\left({\bm{I}}-\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}({\bm{m}})\right)^{-1}$
$\displaystyle\times\left({\bm{A}}{\bm{z}}_{0}-\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}({\bm{m}}){\bm{A}}{\bm{z}}_{0}+\mathbf{\delta}{\bm{A}}{\bm{z}}_{0}\right)=$
$\displaystyle\left(\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}\left({\bm{m}}\right)\right)^{r}\overleftarrow{{\bm{x}}_{0}}+\left({\bm{I}}-\left(\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}({\bm{m}})\right)^{r}\right)\left({\bm{I}}-\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}({\bm{m}})\right)^{-1}$
$\displaystyle\times\left(\left({\bm{I}}-\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}({\bm{m}})\right){\bm{A}}{\bm{z}}_{0}+\mathbf{\delta}{\bm{A}}{\bm{z}}_{0}\right)=$
$\displaystyle\left(\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}\left({\bm{m}}\right)\right)^{r}\overleftarrow{{\bm{x}}_{0}}+\left({\bm{I}}-\left(\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}({\bm{m}})\right)^{r}\right)\left({\bm{A}}{\bm{z}}_{0}+\left({\bm{I}}-\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}({\bm{m}})\right)^{-1}\mathbf{\delta}{\bm{A}}{\bm{z}}_{0}\right)=$
$\displaystyle\left(\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}\left({\bm{m}}\right)\right)^{r}\overleftarrow{{\bm{x}}_{0}}+{\bm{A}}{\bm{z}}_{0}+\left({\bm{I}}-\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}({\bm{m}})\right)^{-1}\mathbf{\delta}{\bm{A}}{\bm{z}}_{0}-\left(\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}({\bm{m}})\right)^{r}{\bm{A}}{\bm{z}}_{0}$
$\displaystyle-\left(\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}({\bm{m}})\right)^{r}\left({\bm{I}}-\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}({\bm{m}})\right)^{-1}\mathbf{\delta}{\bm{A}}{\bm{z}}_{0}.$

Subtracting $\overrightarrow{{\bm{x}}_{0}}={\bm{A}}{\bm{z}}_{0}$ from both
sides of the above expression and taking the Euclidean norm, we get

$\displaystyle\left\|\overleftarrow{{\bm{x}}_{0}^{r}}-\overrightarrow{{\bm{x}}_{0}}\right\|=\Biggl{\|}\left(\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}\left({\bm{m}}\right)\right)^{r}\overleftarrow{{\bm{x}}_{0}}+\left({\bm{I}}-\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}({\bm{m}})\right)^{-1}\mathbf{\delta}{\bm{A}}{\bm{z}}_{0}-\left(\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}({\bm{m}})\right)^{r}{\bm{A}}{\bm{z}}_{0}$
$\displaystyle-\left(\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}({\bm{m}})\right)^{r}\left({\bm{I}}-\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}({\bm{m}})\right)^{-1}\mathbf{\delta}{\bm{A}}{\bm{z}}_{0}\Biggr{\|}.$

After rearranging the terms, we apply the triangle inequality followed by
Cauchy-Schwarz inequality to obtain

$\displaystyle\left\|\overleftarrow{{\bm{x}}_{0}^{r}}-\overrightarrow{{\bm{x}}_{0}}\right\|\leq$
$\displaystyle\left\|\left(\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}\left({\bm{m}}\right)\right)^{r}\right\|\left\|\left(\overleftarrow{{\bm{x}}_{0}}-{\bm{A}}{\bm{z}}_{0}-\left({\bm{I}}-\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}({\bm{m}})\right)^{-1}\mathbf{\delta}{\bm{A}}{\bm{z}}_{0}\right)\right\|$
$\displaystyle+\left\|\left({\bm{I}}-\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}({\bm{m}})\right)^{-1}\mathbf{\delta}{\bm{A}}{\bm{z}}_{0}\right\|\leq$
$\displaystyle\left\|\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}\left({\bm{m}}\right)\right\|^{r}\left(\left\|\overleftarrow{{\bm{x}}_{0}}-{\bm{A}}{\bm{z}}_{0}\right\|+\left\|\left({\bm{I}}-\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}({\bm{m}})\right)^{-1}\mathbf{\delta}{\bm{A}}{\bm{z}}_{0}\right\|\right)$
$\displaystyle+\left\|\left({\bm{I}}-\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}({\bm{m}})\right)^{-1}\mathbf{\delta}{\bm{A}}{\bm{z}}_{0}\right\|.$

Since $\hat{\lambda}_{\max}$ denotes the largest eigen value of
$\left({\bm{A}}{\bm{A}}^{T}+\mathbf{\delta}\right){\bm{D}}\left({\bm{m}}\right)$,
we further simplify as:

$\displaystyle\left\|\overleftarrow{{\bm{x}}_{0}^{r}}-\overrightarrow{{\bm{x}}_{0}}\right\|\leq$
$\displaystyle\hat{\lambda}_{\max}^{r}\left(\left\|\overleftarrow{{\bm{x}}_{0}}-{\bm{A}}{\bm{z}}_{0}\right\|+\left\|\left({\bm{I}}-\hat{\lambda}_{\max}{\bm{I}}\right)^{-1}\mathbf{\delta}{\bm{A}}{\bm{z}}_{0}\right\|\right)+\left\|\left({\bm{I}}-\hat{\lambda}_{\max}{\bm{I}}\right)^{-1}\mathbf{\delta}{\bm{A}}{\bm{z}}_{0}\right\|=$
$\displaystyle\hat{\lambda}_{\max}^{r}\left(\left\|\overleftarrow{{\bm{x}}_{0}}-{\bm{A}}{\bm{z}}_{0}\right\|+\frac{\left\|\mathbf{\delta}{\bm{A}}{\bm{z}}_{0}\right\|}{\left(1-\hat{\lambda}_{\max}\right)}\right)+\frac{\left\|\mathbf{\delta}{\bm{A}}{\bm{z}}_{0}\right\|}{\left(1-\hat{\lambda}_{\max}\right)}.$

Using the fact that
$\overleftarrow{{\bm{x}}_{0}}={\bm{A}}{\bm{A}}^{T}\overleftarrow{{\bm{x}}_{1}}$,
$\overrightarrow{{\bm{x}}_{0}}={\bm{A}}{\bm{z}}_{0}$, and
${\bm{A}}^{T}{\bm{A}}={\bm{I}}$, we get

$\displaystyle\left\|\overleftarrow{{\bm{x}}_{0}^{r}}-\overrightarrow{{\bm{x}}_{0}}\right\|\leq\hat{\lambda}_{\max}^{r}\left(\left\|{\bm{A}}{\bm{A}}^{T}\overleftarrow{{\bm{x}}_{1}}-{\bm{A}}({\bm{A}}^{T}{\bm{A}}){\bm{z}}_{0}\right\|+\frac{\left\|\mathbf{\delta}{\bm{A}}{\bm{z}}_{0}\right\|}{\left(1-\hat{\lambda}_{\max}\right)}\right)+\frac{\left\|\mathbf{\delta}{\bm{A}}{\bm{z}}_{0}\right\|}{\left(1-\hat{\lambda}_{\max}\right)}\leq$
$\displaystyle\hat{\lambda}_{\max}^{r}\left(\left\|{\bm{A}}{\bm{A}}^{T}\right\|\left\|\overleftarrow{{\bm{x}}_{1}}-{\bm{A}}{\bm{z}}_{0}\right\|+\frac{\left\|\mathbf{\delta}\right\|\left\|{\bm{A}}{\bm{z}}_{0}\right\|}{\left(1-\hat{\lambda}_{\max}\right)}\right)+\frac{\left\|\mathbf{\delta}\right\|\left\|{\bm{A}}{\bm{z}}_{0}\right\|}{\left(1-\hat{\lambda}_{\max}\right)}=$
$\displaystyle\hat{\lambda}_{\max}^{r}\left(\frac{\left\|{\bm{\theta}}^{*}\right\|\left\|\overleftarrow{{\bm{x}}_{1}}-\overrightarrow{{\bm{x}}_{0}}\right\|}{\sqrt{1-\beta}}+\frac{\left\|{\bm{\theta}}-{\bm{\theta}}^{*}\right\|\left\|\overrightarrow{{\bm{x}}_{0}}\right\|}{\left(1-\hat{\lambda}_{\max}\right)\sqrt{1-\beta}}\right)+\frac{\left\|{\bm{\theta}}-{\bm{\theta}}^{*}\right\|\left\|\overrightarrow{{\bm{x}}_{0}}\right\|}{\left(1-\hat{\lambda}_{\max}\right)\sqrt{1-\beta}},$

which completes the statement of the theorem. $\square$

###  A.4 Proof of Corollary 7

We begin with the statement of Theorem 6,

$\displaystyle\left\|\overleftarrow{{\bm{x}}_{0}^{r}}-\overrightarrow{{\bm{x}}_{0}}\right\|\leq\hat{\lambda}_{\max}^{r}\left(\frac{\left\|{\bm{\theta}}^{*}\right\|\left\|\overleftarrow{{\bm{x}}_{1}}-\overrightarrow{{\bm{x}}_{0}}\right\|}{\sqrt{1-\beta}}+\frac{\left\|{\bm{\theta}}-{\bm{\theta}}^{*}\right\|\left\|\overrightarrow{{\bm{x}}_{0}}\right\|}{\left(1-\hat{\lambda}_{\max}\right)\sqrt{1-\beta}}\right)+\frac{\left\|{\bm{\theta}}-{\bm{\theta}}^{*}\right\|\left\|\overrightarrow{{\bm{x}}_{0}}\right\|}{\left(1-\hat{\lambda}_{\max}\right)\sqrt{1-\beta}}.$

Since $\hat{\lambda}_{\max}<1$, the first term vanishes in the limit. Using
the fact that $q(\overrightarrow{{\bm{x}}_{0}})$ has a compact support, we get

$\displaystyle\mathop{\sup}_{\begin{subarray}{c}\overrightarrow{{\bm{x}}_{0}}\sim
q(\overrightarrow{{\bm{x}}_{0}})\end{subarray}}\left\|\overleftarrow{{\bm{x}}_{0}^{r}}-\overrightarrow{{\bm{x}}_{0}}\right\|$
$\displaystyle\leq\mathop{\sup}_{\begin{subarray}{c}\overrightarrow{{\bm{x}}_{0}}\sim
q(\overrightarrow{{\bm{x}}_{0}})\end{subarray}}\frac{\left\|{\bm{\theta}}-{\bm{\theta}}^{*}\right\|\left\|\overrightarrow{{\bm{x}}_{0}}\right\|}{\left(1-\hat{\lambda}_{\max}\right)\sqrt{1-\beta}}$
$\displaystyle\leq\frac{\kappa\left\|{\bm{\theta}}-{\bm{\theta}}^{*}\right\|}{\left(1-\hat{\lambda}_{\max}\right)\sqrt{1-\beta}}$
$\displaystyle=\frac{\kappa\delta}{\left(1-\hat{\lambda}_{\max}\right)}.$

Substituting
$\delta={\mathcal{O}}\left(\epsilon\left(1-\hat{\lambda}_{\max}\right)/\kappa\right)$
for a $\delta$-approximate generator, we have
$\mathop{\sup}_{\begin{subarray}{c}\overrightarrow{{\bm{x}}_{0}}\sim
q(\overrightarrow{{\bm{x}}_{0}})\end{subarray}}\left\|\overleftarrow{{\bm{x}}_{0}^{r}}-\overrightarrow{{\bm{x}}_{0}}\right\|\leq\epsilon$,
which finishes the proof. $\square$

###  A.5 Proof of Theorem 8

For diffusion over $T+1$ states, the training objective to minimize the
difference in posterior means is given by (6):

$\displaystyle\min_{{\bm{\theta}}}\mathop{\mathbb{E}}_{\begin{subarray}{c}\overrightarrow{{\bm{x}}_{0}}\sim
q(\overrightarrow{{\bm{x}}_{0}})\\\
\mathbf{\epsilon}\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}\right)\\\
t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\left[\left\|\tilde{\mu}_{t}\left(\overrightarrow{{\bm{x}}_{t}},\overrightarrow{{\bm{x}}_{0}}\right)-\mu_{\theta}\left(\overleftarrow{{\bm{x}}_{t}},t\right)\right\|^{2}\right]\coloneqq\left[\left\|\frac{\sqrt{\bar{\alpha}_{t-1}}\beta_{t}}{1-\bar{a}_{t}}\overrightarrow{{\bm{x}}_{0}}+\frac{\sqrt{\alpha}_{t}\left(1-\bar{\alpha}_{t-1}\right)}{1-\bar{\alpha}_{t}}\overrightarrow{{\bm{x}}_{t}}-\mu_{\bm{\theta}}\left(\overleftarrow{{\bm{x}}_{t}},t\right)\right\|^{2}\right]=$
$\displaystyle\mathop{\mathbb{E}}_{\begin{subarray}{c}\overrightarrow{{\bm{x}}_{0}}\sim
q(\overrightarrow{{\bm{x}}_{0}})\\\
\mathbf{\epsilon}\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}\right)\\\
t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\left[\left\|\frac{\sqrt{\bar{\alpha}_{t-1}}\beta_{t}}{1-\bar{a}_{t}}\overrightarrow{{\bm{x}}_{0}}+\frac{\sqrt{\alpha}_{t}\left(1-\bar{\alpha}_{t-1}\right)}{1-\bar{\alpha}_{t}}\overrightarrow{{\bm{x}}_{t}}-{\bm{\theta}}\begin{bmatrix}\overleftarrow{{\bm{x}}_{t}}\\\
t\end{bmatrix}\right\|^{2}\right].$

We know from our previous analysis that the closed-form solution of the above
MMSE problem is expressed as:

$\displaystyle{\bm{\theta}}_{i}^{*}=\left[\mathop{\mathbb{E}}_{\begin{subarray}{c}\overrightarrow{{\bm{x}}_{0}}\sim
q(\overrightarrow{{\bm{x}}_{0}})\\\
\mathbf{\epsilon}\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}\right)\\\
t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\begin{bmatrix}\overleftarrow{{\bm{x}}_{t}}\\\
t\end{bmatrix}\begin{bmatrix}\overleftarrow{{\bm{x}}_{t}}\\\
t\end{bmatrix}^{T}\right]^{-1}\times\mathop{\mathbb{E}}_{\begin{subarray}{c}\overrightarrow{{\bm{x}}_{0}}\sim
q(\overrightarrow{{\bm{x}}_{0}})\\\
\mathbf{\epsilon}\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}\right)\\\
t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\left[\left(\frac{\sqrt{\bar{\alpha}_{t-1}}\beta_{t}}{1-\bar{a}_{t}}\overrightarrow{{\bm{x}}_{0}}+\frac{\sqrt{\alpha}_{t}\left(1-\bar{\alpha}_{t-1}\right)}{1-\bar{\alpha}_{t}}\overrightarrow{{\bm{x}}_{t}}\right)_{i}\begin{bmatrix}\overleftarrow{{\bm{x}}_{t}}\\\
t\end{bmatrix}\right].$

By substituting the expression for samples generated at intermediate states,
we get the following result:

$\displaystyle{\bm{\theta}}_{i}^{*}=\left[\mathop{\mathbb{E}}_{\begin{subarray}{c}\overrightarrow{{\bm{x}}_{0}}\sim
q(\overrightarrow{{\bm{x}}_{0}})\\\
\mathbf{\epsilon}\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}\right)\\\
t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\begin{bmatrix}\overrightarrow{{\bm{x}}_{0}}\sqrt{\bar{\alpha}_{t}}+\mathbf{\epsilon}\sqrt{1-\bar{\alpha}_{t}}\\\
t\end{bmatrix}\begin{bmatrix}\overrightarrow{{\bm{x}}_{0}}\sqrt{\bar{\alpha}_{t}}+\mathbf{\epsilon}\sqrt{1-\bar{\alpha}_{t}}\\\
t\end{bmatrix}^{T}\right]^{-1}$
$\displaystyle\times\mathop{\mathbb{E}}_{\begin{subarray}{c}\overrightarrow{{\bm{x}}_{0}}\sim
q(\overrightarrow{{\bm{x}}_{0}})\\\
\mathbf{\epsilon}\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}\right)\\\
t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\left[\left(\frac{\sqrt{\bar{\alpha}_{t-1}}\beta_{t}}{1-\bar{a}_{t}}\overrightarrow{{\bm{x}}_{0}}+\frac{\sqrt{\alpha}_{t}\left(1-\bar{\alpha}_{t-1}\right)}{1-\bar{\alpha}_{t}}\left(\overrightarrow{{\bm{x}}_{0}}\sqrt{\bar{\alpha}_{t}}+\mathbf{\epsilon}\sqrt{1-\bar{\alpha}_{t}}\right)\right)_{i}\begin{bmatrix}\overrightarrow{{\bm{x}}_{0}}\sqrt{\bar{\alpha}_{t}}+\mathbf{\epsilon}\sqrt{1-\bar{\alpha}_{t}}\\\
t\end{bmatrix}\right].$

Taking the outer product and expanding further, we arrive at

$\displaystyle{\bm{\theta}}_{i}^{*}=\left[\mathop{\mathbb{E}}_{\begin{subarray}{c}\overrightarrow{{\bm{x}}_{0}}\sim
q(\overrightarrow{{\bm{x}}_{0}})\\\
\mathbf{\epsilon}\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}\right)\\\
t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\begin{bmatrix}\left(\overrightarrow{{\bm{x}}_{0}}\sqrt{\bar{\alpha}_{t}}+\mathbf{\epsilon}\sqrt{1-\bar{\alpha}_{t}}\right)\left(\overrightarrow{{\bm{x}}_{0}}\sqrt{\bar{\alpha}_{t}}+\mathbf{\epsilon}\sqrt{1-\bar{\alpha}_{t}}\right)^{T}&\left(\overrightarrow{{\bm{x}}_{0}}\sqrt{\bar{\alpha}_{t}}+\mathbf{\epsilon}\sqrt{1-\bar{\alpha}_{t}}\right)t\\\
t\left(\overrightarrow{{\bm{x}}_{0}}\sqrt{\bar{\alpha}_{t}}+\mathbf{\epsilon}\sqrt{1-\bar{\alpha}_{t}}\right)&t^{2}\end{bmatrix}\right]^{-1}$
$\displaystyle\times\mathop{\mathbb{E}}_{\begin{subarray}{c}\overrightarrow{{\bm{x}}_{0}}\sim
q(\overrightarrow{{\bm{x}}_{0}})\\\
\mathbf{\epsilon}\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}\right)\\\
t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\left[\left(\overrightarrow{{\bm{x}}_{0}}\left(\frac{\sqrt{\bar{\alpha}_{t-1}}\beta_{t}}{1-\bar{a}_{t}}+\frac{\sqrt{\bar{\alpha}_{t}}\sqrt{\alpha}_{t}\left(1-\bar{\alpha}_{t-1}\right)}{1-\bar{\alpha}_{t}}\right)+\mathbf{\epsilon}\left(\frac{\sqrt{\alpha}_{t}\left(1-\bar{\alpha}_{t-1}\right)\sqrt{1-\bar{\alpha}_{t}}}{1-\bar{\alpha}_{t}}\right)\right)_{i}\begin{bmatrix}\overrightarrow{{\bm{x}}_{0}}\sqrt{\bar{\alpha}_{t}}+\mathbf{\epsilon}\sqrt{1-\bar{\alpha}_{t}}\\\
t\end{bmatrix}\right],$

which upon rearrangement gives the following expression:

$\displaystyle{\bm{\theta}}_{i}^{*}=$
$\displaystyle\left[\mathop{\mathbb{E}}_{\begin{subarray}{c}\overrightarrow{{\bm{x}}_{0}}\sim
q(\overrightarrow{{\bm{x}}_{0}})\\\
\mathbf{\epsilon}\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}\right)\\\
t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\begin{bmatrix}\left(\overrightarrow{{\bm{x}}_{0}}\overrightarrow{{\bm{x}}_{0}}^{T}\bar{\alpha}_{t}+\overrightarrow{{\bm{x}}_{0}}~{}\mathbf{\epsilon}^{T}~{}\sqrt{\bar{\alpha}_{t}}\sqrt{1-\bar{\alpha}_{t}}+\mathbf{\epsilon}~{}\overrightarrow{{\bm{x}}_{0}}^{T}~{}\sqrt{\bar{\alpha}_{t}}\sqrt{1-\bar{\alpha}_{t}}+\mathbf{\epsilon}~{}\mathbf{\epsilon}^{T}\left(1-\bar{\alpha}_{t}\right)\right)&\left(\overrightarrow{{\bm{x}}_{0}}\sqrt{\bar{\alpha}_{t}}+\mathbf{\epsilon}\sqrt{1-\bar{\alpha}_{t}}\right)t\\\
t\left(\overrightarrow{{\bm{x}}_{0}}\sqrt{\bar{\alpha}_{t}}+\mathbf{\epsilon}\sqrt{1-\bar{\alpha}_{t}}\right)&t^{2}\end{bmatrix}\right]^{-1}$
$\displaystyle\times\mathop{\mathbb{E}}_{\begin{subarray}{c}\overrightarrow{{\bm{x}}_{0}}\sim
q(\overrightarrow{{\bm{x}}_{0}})\\\
\mathbf{\epsilon}\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}\right)\\\
t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\left[\left(\overrightarrow{{\bm{x}}_{0}}\left(\frac{\sqrt{\bar{\alpha}_{t-1}}\beta_{t}}{1-\bar{a}_{t}}+\frac{\sqrt{\bar{\alpha}_{t}}\sqrt{\alpha}_{t}\left(1-\bar{\alpha}_{t-1}\right)}{1-\bar{\alpha}_{t}}\right)+\mathbf{\epsilon}\left(\frac{\sqrt{\alpha}_{t}\left(1-\bar{\alpha}_{t-1}\right)\sqrt{1-\bar{\alpha}_{t}}}{1-\bar{\alpha}_{t}}\right)\right)_{i}\begin{bmatrix}\overrightarrow{{\bm{x}}_{0}}\sqrt{\bar{\alpha}_{t}}+\mathbf{\epsilon}\sqrt{1-\bar{\alpha}_{t}}\\\
t\end{bmatrix}\right].$

Recall that $\overrightarrow{{\bm{x}}_{0}}$ and $\mathbf{\epsilon}$ have zero
means. Now, let us use the fact that $\overrightarrow{{\bm{x}}_{0}}$,
$\mathbf{\epsilon}$, and $t$ are all independent random variables.

$\displaystyle{\bm{\theta}}_{i}^{*}=\left[\mathop{\mathbb{E}}_{\begin{subarray}{c}{\bm{z}}_{0}\sim{\mathcal{N}}(\mathbf{0},{\bm{I}})\\\
\mathbf{\epsilon}\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}\right)\\\
t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\begin{bmatrix}\left({\bm{A}}{\bm{z}}_{0}\left({\bm{A}}{\bm{z}}_{0}\right)^{T}\bar{\alpha}_{t}+\mathbf{\epsilon}~{}\mathbf{\epsilon}^{T}\left(1-\bar{\alpha}_{t}\right)\right)&\mathbf{0}\\\
\mathbf{0}&t^{2}\end{bmatrix}\right]^{-1}$
$\displaystyle\times\mathop{\mathbb{E}}_{\begin{subarray}{c}\overrightarrow{{\bm{x}}_{0}}\sim
q(\overrightarrow{{\bm{x}}_{0}})\\\
\mathbf{\epsilon}\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}\right)\\\
t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\left[\left(\overrightarrow{{\bm{x}}_{0}}\left(\frac{\sqrt{\bar{\alpha}_{t-1}}\beta_{t}}{1-\bar{a}_{t}}+\frac{\sqrt{\bar{\alpha}_{t}}\sqrt{\alpha}_{t}\left(1-\bar{\alpha}_{t-1}\right)}{1-\bar{\alpha}_{t}}\right)+\mathbf{\epsilon}\left(\frac{\sqrt{\alpha}_{t}\left(1-\bar{\alpha}_{t-1}\right)\sqrt{1-\bar{\alpha}_{t}}}{1-\bar{\alpha}_{t}}\right)\right)_{i}\begin{bmatrix}\overrightarrow{{\bm{x}}_{0}}\sqrt{\bar{\alpha}_{t}}+\mathbf{\epsilon}\sqrt{1-\bar{\alpha}_{t}}\\\
t\end{bmatrix}\right].$

Since $\alpha_{t}=1-\beta_{t}$ and
$\sqrt{\bar{\alpha}_{t}}=\sqrt{\bar{\alpha}_{t-1}}\sqrt{\alpha_{t}}$, we
further simplify the closed-form solution as follows:

$\displaystyle{\bm{\theta}}_{i}^{*}=$
$\displaystyle\left[\mathop{\mathbb{E}}_{\begin{subarray}{c}{\bm{z}}_{0}\sim{\mathcal{N}}(\mathbf{0},{\bm{I}})\\\
\mathbf{\epsilon}\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}\right)\\\
t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\begin{bmatrix}\left({\bm{A}}{\bm{z}}_{0}{\bm{z}}_{0}^{T}{\bm{A}}^{T}\bar{\alpha}_{t}+\mathbf{\epsilon}~{}\mathbf{\epsilon}^{T}\left(1-\bar{\alpha}_{t}\right)\right)&\mathbf{0}\\\
\mathbf{0}&t^{2}\end{bmatrix}\right]^{-1}$
$\displaystyle\times\mathop{\mathbb{E}}_{\begin{subarray}{c}\overrightarrow{{\bm{x}}_{0}}\sim
q(\overrightarrow{{\bm{x}}_{0}})\\\
\mathbf{\epsilon}\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}\right)\\\
t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\left[\left(\overrightarrow{{\bm{x}}_{0}}\left(\frac{\sqrt{\bar{\alpha}_{t-1}}\left(\beta_{t}+\alpha_{t}\left(1-\bar{\alpha}_{t-1}\right)\right)}{1-\bar{a}_{t}}\right)+\mathbf{\epsilon}\left(\frac{\sqrt{\alpha}_{t}\left(1-\bar{\alpha}_{t-1}\right)}{\sqrt{1-\bar{\alpha}_{t}}}\right)\right)_{i}\begin{bmatrix}\overrightarrow{{\bm{x}}_{0}}\sqrt{\bar{\alpha}_{t}}+\mathbf{\epsilon}\sqrt{1-\bar{\alpha}_{t}}\\\
t\end{bmatrix}\right]=$
$\displaystyle\left[\mathop{\mathbb{E}}_{\begin{subarray}{c}{\bm{z}}_{0}\sim{\mathcal{N}}(\mathbf{0},{\bm{I}})\\\
\mathbf{\epsilon}\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}\right)\\\
t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\begin{bmatrix}\left({\bm{A}}{\bm{A}}^{T}\bar{\alpha}_{t}+\mathbf{\epsilon}~{}\mathbf{\epsilon}^{T}\left(1-\bar{\alpha}_{t}\right)\right)&\mathbf{0}\\\
\mathbf{0}&t^{2}\end{bmatrix}\right]^{-1}$
$\displaystyle\times\mathop{\mathbb{E}}_{\begin{subarray}{c}\overrightarrow{{\bm{x}}_{0}}\sim
q(\overrightarrow{{\bm{x}}_{0}})\\\
\mathbf{\epsilon}\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}\right)\\\
t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\left[\left(\overrightarrow{{\bm{x}}_{0}}\left(\frac{\sqrt{\bar{\alpha}_{t-1}}\left(\beta_{t}+\alpha_{t}-\bar{\alpha}_{t}\right)}{1-\bar{a}_{t}}\right)+\mathbf{\epsilon}\left(\frac{\sqrt{\alpha}_{t}\left(1-\bar{\alpha}_{t-1}\right)}{\sqrt{1-\bar{\alpha}_{t}}}\right)\right)_{i}\begin{bmatrix}\overrightarrow{{\bm{x}}_{0}}\sqrt{\bar{\alpha}_{t}}+\mathbf{\epsilon}\sqrt{1-\bar{\alpha}_{t}}\\\
t\end{bmatrix}\right]=$
$\displaystyle\mathop{\mathbb{E}}_{\begin{subarray}{c}t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\begin{bmatrix}\left({\bm{A}}{\bm{A}}^{T}\bar{\alpha}_{t}+{\bm{I}}\left(1-\bar{\alpha}_{t}\right)\right)&\mathbf{0}\\\
\mathbf{0}&\frac{T^{3}-1}{3T}\end{bmatrix}^{-1}$
$\displaystyle\times\mathop{\mathbb{E}}_{\begin{subarray}{c}\overrightarrow{{\bm{x}}_{0}}\sim
q(\overrightarrow{{\bm{x}}_{0}})\\\
\mathbf{\epsilon}\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}\right)\\\
t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\left[\left(\overrightarrow{{\bm{x}}_{0}}\left(\frac{\sqrt{\bar{\alpha}_{t-1}}\left(\beta_{t}+1-\beta_{t}-\bar{\alpha}_{t}\right)}{1-\bar{a}_{t}}\right)+\mathbf{\epsilon}\left(\frac{\sqrt{\alpha}_{t}\left(1-\bar{\alpha}_{t-1}\right)}{\sqrt{1-\bar{\alpha}_{t}}}\right)\right)_{i}\begin{bmatrix}\overrightarrow{{\bm{x}}_{0}}\sqrt{\bar{\alpha}_{t}}+\mathbf{\epsilon}\sqrt{1-\bar{\alpha}_{t}}\\\
t\end{bmatrix}\right]=$
$\displaystyle\mathop{\mathbb{E}}_{\begin{subarray}{c}t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\begin{bmatrix}\left[{\bm{A}}{\bm{A}}^{T}\bar{\alpha}_{t}+{\bm{I}}\left(1-\bar{\alpha}_{t}\right)\right]^{-1}&\mathbf{0}\\\
\mathbf{0}&\frac{3T}{T^{3}-1}\end{bmatrix}$
$\displaystyle\times\mathop{\mathbb{E}}_{\begin{subarray}{c}\overrightarrow{{\bm{x}}_{0}}\sim
q(\overrightarrow{{\bm{x}}_{0}})\\\
\mathbf{\epsilon}\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}\right)\\\
t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\left[\left(\overrightarrow{{\bm{x}}_{0}}\sqrt{\bar{\alpha}_{t-1}}+\mathbf{\epsilon}\left(\frac{\sqrt{\alpha}_{t}\left(1-\bar{\alpha}_{t-1}\right)}{\sqrt{1-\bar{\alpha}_{t}}}\right)\right)_{i}\begin{bmatrix}\overrightarrow{{\bm{x}}_{0}}\sqrt{\bar{\alpha}_{t}}+\mathbf{\epsilon}\sqrt{1-\bar{\alpha}_{t}}\\\
t\end{bmatrix}\right].$

Using Woodburry matrix inverse [Hag89] in the above expression, we get

$\displaystyle{\bm{\theta}}_{i}^{*}=\mathop{\mathbb{E}}_{\begin{subarray}{c}t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\begin{bmatrix}\left(1-\bar{\alpha}_{t}\right)^{-1}\left[{\bm{I}}+\left({\bm{A}}\sqrt{\frac{\bar{\alpha}_{t}}{\left(1-\bar{\alpha}_{t}\right)}}\right)\left({\bm{A}}\sqrt{\frac{\bar{\alpha}_{t}}{\left(1-\bar{\alpha}_{t}\right)}}\right)^{T}\right]^{-1}&\mathbf{0}\\\
\mathbf{0}&\frac{3T}{T^{3}-1}\end{bmatrix}$
$\displaystyle\times\mathop{\mathbb{E}}_{\begin{subarray}{c}{\bm{z}}_{0}\sim{\mathcal{N}}(\mathbf{0},{\bm{I}})\\\
\mathbf{\epsilon}\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}\right)\\\
t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\left[\left({\bm{A}}{\bm{z}}_{0}\sqrt{\bar{\alpha}_{t-1}}+\mathbf{\epsilon}\left(\frac{\sqrt{\alpha}_{t}\left(1-\bar{\alpha}_{t-1}\right)}{\sqrt{1-\bar{\alpha}_{t}}}\right)\right)_{i}\begin{bmatrix}{\bm{A}}{\bm{z}}_{0}\sqrt{\bar{\alpha}_{t}}+\mathbf{\epsilon}\sqrt{1-\bar{\alpha}_{t}}\\\
t\end{bmatrix}\right]=$
$\displaystyle\mathop{\mathbb{E}}_{\begin{subarray}{c}t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\begin{bmatrix}\left(1-\bar{\alpha}_{t}\right)^{-1}\left[{\bm{I}}-\left({\bm{A}}\sqrt{\frac{\bar{\alpha}_{t}}{\left(1-\bar{\alpha}_{t}\right)}}\right)\left({\bm{I}}+\left({\bm{A}}\sqrt{\frac{\bar{\alpha}_{t}}{\left(1-\bar{\alpha}_{t}\right)}}\right)^{T}\left({\bm{A}}\sqrt{\frac{\bar{\alpha}_{t}}{\left(1-\bar{\alpha}_{t}\right)}}\right)\right)^{-1}\left({\bm{A}}\sqrt{\frac{\bar{\alpha}_{t}}{\left(1-\bar{\alpha}_{t}\right)}}\right)^{T}\right]^{-1}&\mathbf{0}\\\
\mathbf{0}&\frac{3T}{T^{3}-1}\end{bmatrix}$
$\displaystyle\times\mathop{\mathbb{E}}_{\begin{subarray}{c}{\bm{z}}_{0}\sim{\mathcal{N}}(\mathbf{0},{\bm{I}})\\\
\mathbf{\epsilon}\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}\right)\\\
t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\left[\left({\bm{A}}{\bm{z}}_{0}\sqrt{\bar{\alpha}_{t-1}}+\mathbf{\epsilon}\left(\frac{\sqrt{\alpha}_{t}\left(1-\bar{\alpha}_{t-1}\right)}{\sqrt{1-\bar{\alpha}_{t}}}\right)\right)_{i}\begin{bmatrix}{\bm{A}}{\bm{z}}_{0}\sqrt{\bar{\alpha}_{t}}+\mathbf{\epsilon}\sqrt{1-\bar{\alpha}_{t}}\\\
t\end{bmatrix}\right].$

Next, we simplify the second term further by extracting the $i^{th}$ component
$\left({\bm{A}}{\bm{z}}_{0}\sqrt{\bar{\alpha}_{t-1}}+\mathbf{\epsilon}\left(\frac{\sqrt{\alpha}_{t}\left(1-\bar{\alpha}_{t-1}\right)}{\sqrt{1-\bar{\alpha}_{t}}}\right)\right)_{i}$
as:

$\displaystyle{\bm{\theta}}_{i}^{*}=$
$\displaystyle\mathop{\mathbb{E}}_{\begin{subarray}{c}t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\begin{bmatrix}\left(1-\bar{\alpha}_{t}\right)^{-1}\left[{\bm{I}}-\left({\bm{A}}\sqrt{\frac{\bar{\alpha}_{t}}{\left(1-\bar{\alpha}_{t}\right)}}\right)\left({\bm{I}}+\left({\bm{A}}\sqrt{\frac{\bar{\alpha}_{t}}{\left(1-\bar{\alpha}_{t}\right)}}\right)^{T}\left({\bm{A}}\sqrt{\frac{\bar{\alpha}_{t}}{\left(1-\bar{\alpha}_{t}\right)}}\right)\right)^{-1}\left({\bm{A}}\sqrt{\frac{\bar{\alpha}_{t}}{\left(1-\bar{\alpha}_{t}\right)}}\right)^{T}\right]^{-1}&\mathbf{0}\\\
\mathbf{0}&\frac{3T}{T^{3}-1}\end{bmatrix}$
$\displaystyle\times\mathop{\mathbb{E}}_{\begin{subarray}{c}{\bm{z}}_{0}\sim{\mathcal{N}}(\mathbf{0},{\bm{I}})\\\
\mathbf{\epsilon}\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}\right)\\\
t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\left[\left({\bm{a}}_{i}^{T}{\bm{z}}_{0}\sqrt{\bar{\alpha}_{t-1}}+\mathbf{\epsilon}_{i}\left(\frac{\sqrt{\alpha}_{t}\left(1-\bar{\alpha}_{t-1}\right)}{\sqrt{1-\bar{\alpha}_{t}}}\right)\right)\begin{bmatrix}{\bm{A}}{\bm{z}}_{0}\sqrt{\bar{\alpha}_{t}}+\mathbf{\epsilon}\sqrt{1-\bar{\alpha}_{t}}\\\
t\end{bmatrix}\right]=$
$\displaystyle\mathop{\mathbb{E}}_{\begin{subarray}{c}t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\begin{bmatrix}\left(1-\bar{\alpha}_{t}\right)^{-1}\left[{\bm{I}}-\left({\bm{A}}\sqrt{\frac{\bar{\alpha}_{t}}{\left(1-\bar{\alpha}_{t}\right)}}\right)\left({\bm{I}}+{\bm{I}}\left(\frac{\bar{\alpha}_{t}}{1-\bar{\alpha}_{t}}\right)\right)^{-1}\left({\bm{A}}\sqrt{\frac{\bar{\alpha}_{t}}{\left(1-\bar{\alpha}_{t}\right)}}\right)^{T}\right]^{-1}&\mathbf{0}\\\
\mathbf{0}&\frac{3T}{T^{3}-1}\end{bmatrix}$
$\displaystyle\times\mathop{\mathbb{E}}_{\begin{subarray}{c}{\bm{z}}_{0}\sim{\mathcal{N}}(\mathbf{0},{\bm{I}})\\\
\mathbf{\epsilon}\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}\right)\\\
t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\begin{bmatrix}\left({\bm{a}}_{i}^{T}{\bm{z}}_{0}\sqrt{\bar{\alpha}_{t-1}}+\mathbf{\epsilon}_{i}\left(\frac{\sqrt{\alpha}_{t}\left(1-\bar{\alpha}_{t-1}\right)}{\sqrt{1-\bar{\alpha}_{t}}}\right)\right)\left({\bm{A}}{\bm{z}}_{0}\sqrt{\bar{\alpha}_{t}}+\mathbf{\epsilon}\sqrt{1-\bar{\alpha}_{t}}\right)\\\
\left({\bm{a}}_{i}^{T}{\bm{z}}_{0}\sqrt{\bar{\alpha}_{t-1}}+\mathbf{\epsilon}_{i}\left(\frac{\sqrt{\alpha}_{t}\left(1-\bar{\alpha}_{t-1}\right)}{\sqrt{1-\bar{\alpha}_{t}}}\right)\right)t\end{bmatrix}=$
$\displaystyle\mathop{\mathbb{E}}_{\begin{subarray}{c}t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\begin{bmatrix}\left(1-\bar{\alpha}_{t}\right)^{-1}\left[{\bm{I}}-\left({\bm{A}}\sqrt{\frac{\bar{\alpha}_{t}}{\left(1-\bar{\alpha}_{t}\right)}}\right){\bm{I}}\left(1-\bar{\alpha}_{t}\right)\left({\bm{A}}\sqrt{\frac{\bar{\alpha}_{t}}{\left(1-\bar{\alpha}_{t}\right)}}\right)^{T}\right]&\mathbf{0}\\\
\mathbf{0}&\frac{3T}{T^{3}-1}\end{bmatrix}$
$\displaystyle\times\mathop{\mathbb{E}}_{\begin{subarray}{c}{\bm{z}}_{0}\sim{\mathcal{N}}(\mathbf{0},{\bm{I}})\\\
\mathbf{\epsilon}\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}\right)\\\
t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\begin{bmatrix}\left({\bm{A}}{\bm{z}}_{0}({\bm{a}}_{i}^{T}{\bm{z}}_{0})\sqrt{\bar{\alpha}_{t-1}}\sqrt{\bar{\alpha}_{t}}+\mathbf{\epsilon}_{i}\left(\frac{\sqrt{\alpha}_{t}\left(1-\bar{\alpha}_{t-1}\right)}{\sqrt{1-\bar{\alpha}_{t}}}\right)\mathbf{\epsilon}\sqrt{1-\bar{\alpha}_{t}}\right)\\\
0\end{bmatrix}=$
$\displaystyle\mathop{\mathbb{E}}_{\begin{subarray}{c}t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\begin{bmatrix}\left(1-\bar{\alpha}_{t}\right)^{-1}\left[{\bm{I}}-{\bm{A}}{\bm{A}}^{T}\bar{\alpha}_{t}\right]&\mathbf{0}\\\
\mathbf{0}&\frac{3T}{T^{3}-1}\end{bmatrix}$
$\displaystyle\times\mathop{\mathbb{E}}_{\begin{subarray}{c}\mathbf{\epsilon}\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}\right)\\\
t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\begin{bmatrix}\left({\bm{A}}{\bm{a}}_{i}\bar{\alpha}_{t-1}\sqrt{\alpha_{t}}+\mathbf{\epsilon}_{i}\left(\sqrt{\alpha}_{t}\left(1-\bar{\alpha}_{t-1}\right)\right)\mathbf{\epsilon}\right)\\\
0\end{bmatrix}=$
$\displaystyle\mathop{\mathbb{E}}_{\begin{subarray}{c}t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\begin{bmatrix}\frac{1}{\left(1-\bar{\alpha}_{t}\right)}\left[{\bm{I}}-{\bm{A}}{\bm{A}}^{T}\bar{\alpha}_{t}\right]&\mathbf{0}\\\
\mathbf{0}&\frac{3T}{T^{3}-1}\end{bmatrix}$
$\displaystyle\times\mathop{\mathbb{E}}_{\begin{subarray}{c}t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\begin{bmatrix}{\bm{A}}{\bm{a}}_{i}\bar{\alpha}_{t-1}\sqrt{\alpha_{t}}+{\bm{e}}_{i}\left(\sqrt{\alpha}_{t}\left(1-\bar{\alpha}_{t-1}\right)\right)\\\
0\end{bmatrix},$

where ${\bm{e}}_{i}\in\\{0,1\\}^{d}$ is a one-hot encoded vector with only
$i^{th}$ coordinate taking the value $1$. Bringing the expectation inside and
using matrix products, we obtain

$\displaystyle{\bm{\theta}}_{i}^{*}=\begin{bmatrix}\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\frac{1}{\left(1-\bar{\alpha}_{t}\right)}\right){\bm{I}}-\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\frac{\bar{\alpha}_{t}}{\left(1-\bar{\alpha}_{t}\right)}\right){\bm{A}}{\bm{A}}^{T}&\mathbf{0}\\\
\mathbf{0}&\frac{3T}{T^{3}-1}\end{bmatrix}\begin{bmatrix}\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\bar{\alpha}_{t-1}\sqrt{\alpha_{t}}\right){\bm{A}}{\bm{a}}_{i}+\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\sqrt{\alpha}_{t}\left(1-\bar{\alpha}_{t-1}\right)\right){\bm{e}}_{i}\\\
0\end{bmatrix}=$
$\displaystyle\begin{bmatrix}\left(\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\frac{1}{\left(1-\bar{\alpha}_{t}\right)}\right){\bm{I}}-\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\frac{\bar{\alpha}_{t}}{\left(1-\bar{\alpha}_{t}\right)}\right){\bm{A}}{\bm{A}}^{T}\right)\left(\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\bar{\alpha}_{t-1}\sqrt{\alpha_{t}}\right){\bm{A}}{\bm{a}}_{i}+\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\sqrt{\alpha}_{t}\left(1-\bar{\alpha}_{t-1}\right)\right){\bm{e}}_{i}\right)\\\
0\end{bmatrix}=$
$\displaystyle\begin{bmatrix}\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\frac{1}{\left(1-\bar{\alpha}_{t}\right)}\right)\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\bar{\alpha}_{t-1}\sqrt{\alpha_{t}}\right){\bm{A}}{\bm{a}}_{i}-\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\frac{\bar{\alpha}_{t}}{\left(1-\bar{\alpha}_{t}\right)}\right)\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\bar{\alpha}_{t-1}\sqrt{\alpha_{t}}\right){\bm{A}}{\bm{A}}^{T}{\bm{A}}{\bm{a}}_{i}\\\
+\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\frac{1}{\left(1-\bar{\alpha}_{t}\right)}\right)\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\sqrt{\alpha}_{t}\left(1-\bar{\alpha}_{t-1}\right)\right){\bm{e}}_{i}-\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\frac{\bar{\alpha}_{t}}{\left(1-\bar{\alpha}_{t}\right)}\right)\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\sqrt{\alpha}_{t}\left(1-\bar{\alpha}_{t-1}\right)\right)\left({\bm{A}}{\bm{A}}^{T}\right){\bm{e}}_{i}\\\
0\end{bmatrix}=$
$\displaystyle\begin{bmatrix}\left(\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\frac{1}{\left(1-\bar{\alpha}_{t}\right)}\right)\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\bar{\alpha}_{t-1}\sqrt{\alpha_{t}}\right)-\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\frac{\bar{\alpha}_{t}}{\left(1-\bar{\alpha}_{t}\right)}\right)\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\bar{\alpha}_{t-1}\sqrt{\alpha_{t}}\right)\par\right){\bm{A}}{\bm{a}}_{i}\\\
+\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\frac{1}{\left(1-\bar{\alpha}_{t}\right)}\right)\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\sqrt{\alpha}_{t}\left(1-\bar{\alpha}_{t-1}\right)\right){\bm{e}}_{i}-\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\frac{\bar{\alpha}_{t}}{\left(1-\bar{\alpha}_{t}\right)}\right)\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\sqrt{\alpha}_{t}\left(1-\bar{\alpha}_{t-1}\right)\right)\left({\bm{A}}{\bm{A}}^{T}\right){\bm{e}}_{i}\\\
0\end{bmatrix}.$

Stacking all the rows together, we obtain the closed-form solution for
generative modeling:

$\displaystyle{\bm{\theta}}^{*}=\begin{bmatrix}\left(\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\frac{1}{\left(1-\bar{\alpha}_{t}\right)}\right)\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\bar{\alpha}_{t-1}\sqrt{\alpha_{t}}\right)-\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\frac{\bar{\alpha}_{t}}{\left(1-\bar{\alpha}_{t}\right)}\right)\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\bar{\alpha}_{t-1}\sqrt{\alpha_{t}}\right)\right){\bm{A}}{\bm{A}}^{T}\\\
+\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\frac{1}{\left(1-\bar{\alpha}_{t}\right)}\right)\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\sqrt{\alpha}_{t}\left(1-\bar{\alpha}_{t-1}\right)\right){\bm{I}}\\\
-\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\frac{\bar{\alpha}_{t}}{\left(1-\bar{\alpha}_{t}\right)}\right)\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\sqrt{\alpha}_{t}\left(1-\bar{\alpha}_{t-1}\right)\right){\bm{A}}{\bm{A}}^{T}&\mathbf{0}\end{bmatrix}=$
$\displaystyle\begin{bmatrix}\left(\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\frac{1}{\left(1-\bar{\alpha}_{t}\right)}\right)\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\bar{\alpha}_{t-1}\sqrt{\alpha_{t}}\right)-\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\frac{\bar{\alpha}_{t}}{\left(1-\bar{\alpha}_{t}\right)}\right)\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\bar{\alpha}_{t-1}\sqrt{\alpha_{t}}\right)\right){\bm{A}}{\bm{A}}^{T}\\\
+\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\frac{1}{\left(1-\bar{\alpha}_{t}\right)}\right)\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\sqrt{\alpha}_{t}\left(1-\bar{\alpha}_{t-1}\right)\right){\bm{I}}\\\
-\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\frac{\bar{\alpha}_{t}}{\left(1-\bar{\alpha}_{t}\right)}\right)\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\sqrt{\alpha}_{t}\right){\bm{A}}{\bm{A}}^{T}+\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\frac{\bar{\alpha}_{t}}{\left(1-\bar{\alpha}_{t}\right)}\right)\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\sqrt{\alpha}_{t}\bar{\alpha}_{t-1}\right){\bm{A}}{\bm{A}}^{T}&\mathbf{0}\end{bmatrix}=$
$\displaystyle\begin{bmatrix}\left(\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\frac{1}{\left(1-\bar{\alpha}_{t}\right)}\right)\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\bar{\alpha}_{t-1}\sqrt{\alpha_{t}}\right)-\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\frac{\bar{\alpha}_{t}}{\left(1-\bar{\alpha}_{t}\right)}\right)\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\sqrt{\alpha}_{t}\right)\right){\bm{A}}{\bm{A}}^{T}\\\
+\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\frac{1}{\left(1-\bar{\alpha}_{t}\right)}\right)\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\sqrt{\alpha}_{t}\left(1-\bar{\alpha}_{t-1}\right)\right){\bm{I}}&\mathbf{0}\end{bmatrix}=$
$\displaystyle\begin{bmatrix}\nu{\bm{A}}{\bm{A}}^{T}+\gamma{\bm{I}}&\mathbf{0}\end{bmatrix}.$

This completes the proof of the theorem. $\square$

###  A.6 Proof of Corollary 10

Using the arguments from Appendix A.5, we express the closed-form solution of
the MMSE problem (A.5) as:

$\displaystyle{\bm{\theta}}_{i}^{*}=\left[\mathop{\mathbb{E}}_{\begin{subarray}{c}\overrightarrow{{\bm{x}}_{0}}\sim
q(\overrightarrow{{\bm{x}}_{0}})\\\
\mathbf{\epsilon}\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}\right)\\\
t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\begin{bmatrix}\overleftarrow{{\bm{x}}_{t}}\\\
t\end{bmatrix}\begin{bmatrix}\overleftarrow{{\bm{x}}_{t}}\\\
t\end{bmatrix}^{T}\right]^{-1}\times\mathop{\mathbb{E}}_{\begin{subarray}{c}\overrightarrow{{\bm{x}}_{0}}\sim
q(\overrightarrow{{\bm{x}}_{0}})\\\
\mathbf{\epsilon}\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}\right)\\\
t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\left[\left(\frac{\sqrt{\bar{\alpha}_{t-1}}\beta_{t}}{1-\bar{a}_{t}}\overrightarrow{{\bm{x}}_{0}}+\frac{\sqrt{\alpha}_{t}\left(1-\bar{\alpha}_{t-1}\right)}{1-\bar{\alpha}_{t}}\overrightarrow{{\bm{x}}_{t}}\right)_{i}\begin{bmatrix}\overleftarrow{{\bm{x}}_{t}}\\\
t\end{bmatrix}\right]=$
$\displaystyle\mathop{\mathbb{E}}_{\begin{subarray}{c}t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\begin{bmatrix}\left[{\bm{A}}{\bm{A}}^{T}\bar{\alpha}_{t}+{\bm{I}}\left(1-\bar{\alpha}_{t}\right)\right]^{-1}&\mathbf{0}\\\
\mathbf{0}&\frac{3T}{T^{3}-1}\end{bmatrix}$
$\displaystyle\times\mathop{\mathbb{E}}_{\begin{subarray}{c}\overrightarrow{{\bm{x}}_{0}}\sim
q(\overrightarrow{{\bm{x}}_{0}})\\\
\mathbf{\epsilon}\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}\right)\\\
t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\left[\left(\overrightarrow{{\bm{x}}_{0}}\sqrt{\bar{\alpha}_{t-1}}+\overrightarrow{\mathbf{\epsilon}}\left(\frac{\sqrt{\alpha}_{t}\left(1-\bar{\alpha}_{t-1}\right)}{\sqrt{1-\bar{\alpha}_{t}}}\right)\right)_{i}\begin{bmatrix}\overrightarrow{{\bm{x}}_{0}}\sqrt{\bar{\alpha}_{t}}+\overleftarrow{\mathbf{\epsilon}}\sqrt{1-\bar{\alpha}_{t}}\\\
t\end{bmatrix}\right]=$
$\displaystyle\mathop{\mathbb{E}}_{\begin{subarray}{c}t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\begin{bmatrix}\left(1-\bar{\alpha}_{t}\right)^{-1}\left[{\bm{I}}-\left({\bm{A}}\sqrt{\frac{\bar{\alpha}_{t}}{\left(1-\bar{\alpha}_{t}\right)}}\right){\bm{I}}\left(1-\bar{\alpha}_{t}\right)\left({\bm{A}}\sqrt{\frac{\bar{\alpha}_{t}}{\left(1-\bar{\alpha}_{t}\right)}}\right)^{T}\right]&\mathbf{0}\\\
\mathbf{0}&\frac{3T}{T^{3}-1}\end{bmatrix}$
$\displaystyle\times\mathop{\mathbb{E}}_{\begin{subarray}{c}{\bm{z}}_{0}\sim{\mathcal{N}}(\mathbf{0},{\bm{I}})\\\
\mathbf{\epsilon}\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}\right)\\\
t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\begin{bmatrix}\left({\bm{A}}{\bm{z}}_{0}({\bm{a}}_{i}^{T}{\bm{z}}_{0})\sqrt{\bar{\alpha}_{t-1}}\sqrt{\bar{\alpha}_{t}}+\overrightarrow{\mathbf{\epsilon}_{i}}\left(\frac{\sqrt{\alpha}_{t}\left(1-\bar{\alpha}_{t-1}\right)}{\sqrt{1-\bar{\alpha}_{t}}}\right)\overleftarrow{\mathbf{\epsilon}}\sqrt{1-\bar{\alpha}_{t}}\right)\\\
0\end{bmatrix}=$
$\displaystyle\mathop{\mathbb{E}}_{\begin{subarray}{c}t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\begin{bmatrix}\left(1-\bar{\alpha}_{t}\right)^{-1}\left[{\bm{I}}-{\bm{A}}{\bm{A}}^{T}\bar{\alpha}_{t}\right]&\mathbf{0}\\\
\mathbf{0}&\frac{3T}{T^{3}-1}\end{bmatrix}\mathop{\mathbb{E}}_{\begin{subarray}{c}t\sim{\mathcal{U}}\left(1,\dots,T\right)\end{subarray}}\begin{bmatrix}{\bm{A}}{\bm{a}}_{i}\bar{\alpha}_{t-1}\sqrt{\alpha_{t}}\\\
0\end{bmatrix}=$
$\displaystyle\begin{bmatrix}\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\frac{1}{\left(1-\bar{\alpha}_{t}\right)}\right){\bm{I}}-\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\frac{\bar{\alpha}_{t}}{\left(1-\bar{\alpha}_{t}\right)}\right){\bm{A}}{\bm{A}}^{T}&\mathbf{0}\\\
\mathbf{0}&\frac{3T}{T^{3}-1}\end{bmatrix}\begin{bmatrix}\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\bar{\alpha}_{t-1}\sqrt{\alpha_{t}}\right){\bm{A}}{\bm{a}}_{i}\\\
0\end{bmatrix}=$
$\displaystyle\begin{bmatrix}\left(\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\frac{1}{\left(1-\bar{\alpha}_{t}\right)}\right)\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\bar{\alpha}_{t-1}\sqrt{\alpha_{t}}\right)-\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\frac{\bar{\alpha}_{t}}{\left(1-\bar{\alpha}_{t}\right)}\right)\mathop{\mathbb{E}}_{\begin{subarray}{c}t\end{subarray}}\left(\bar{\alpha}_{t-1}\sqrt{\alpha_{t}}\right)\par\right){\bm{A}}{\bm{a}}_{i}\\\
0\end{bmatrix}$

Now, stacking all the rows together as in Appendix A.5, we obtain the closed-
form solution for DDPM generator as
${\bm{\theta}}^{*}=\begin{bmatrix}\bar{\nu}{\bm{A}}{\bm{A}}^{T}&\mathbf{0}\end{bmatrix}$,
which finishes the proof of the corollary. $\square$

###  A.7 Proof of Theorem 11

To prove DDPM learns the distribution
$q\left(\overrightarrow{{\bm{x}}_{0}}\right)$, it is sufficient to show that
$\overleftarrow{{\bm{x}}_{0}}\sim p\left(\overleftarrow{{\bm{x}}_{0}}\right)$
has the structure ${\bm{A}}\overleftarrow{{\bm{z}}_{0}}$, where
$\overleftarrow{{\bm{z}}_{0}}\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}\right)$.
Using the closed-form solution from Corollary 10, we run one step of the
reverse process starting from
$\overleftarrow{{\bm{x}}_{T}}\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}\right)$.
Then, we rectify the misalignment in drift and dispersion of the reverse
process. After recursively applying the reverse step for all the states, we
obtain

$\displaystyle\overleftarrow{{\bm{x}}_{0}}=\left(\prod_{t=1}^{T}\left(\omega_{t}\bar{\nu}{\bm{A}}{\bm{A}}^{T}\right)\right)\overleftarrow{{\bm{x}}_{T}}+\left(\prod_{t=1}^{T-1}\left(\omega_{t}\bar{\nu}{\bm{A}}{\bm{A}}^{T}\right)\right)\xi_{T}\sqrt{\beta_{T}}\overleftarrow{\epsilon_{T}}+\left(\prod_{t=1}^{T-2}\left(\omega_{t}\bar{\nu}{\bm{A}}{\bm{A}}^{T}\right)\right)\xi_{T}\sqrt{\beta_{T-1}}\overleftarrow{\epsilon_{T-1}}$
$\displaystyle+\cdots+\left(\omega_{1}\bar{\nu}{\bm{A}}{\bm{A}}^{T}\right)\left(\omega_{2}\bar{\nu}{\bm{A}}{\bm{A}}^{T}\right)\xi_{3}\sqrt{\beta_{3}}\overleftarrow{\mathbf{\epsilon}_{3}}+\left(\omega_{1}\bar{\nu}{\bm{A}}{\bm{A}}^{T}\right)\xi_{2}\sqrt{\beta_{2}}\overleftarrow{\mathbf{\epsilon}_{2}}.$

Since $\\{\overleftarrow{\mathbf{\epsilon}_{t}}\\}_{t>1}$ are IID Gaussian
random variables, we rewrite the above epxression as:

$\displaystyle\overleftarrow{{\bm{x}}_{0}}=\left(\prod_{t=1}^{T}\left(\omega_{t}\bar{\nu}{\bm{A}}{\bm{A}}^{T}\right)\right)\overleftarrow{{\bm{x}}_{T}}+\bar{\Sigma}^{1/2}\overleftarrow{\mathbf{\epsilon}},$

where
$\overleftarrow{\mathbf{\epsilon}}\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}\right)$
and $\bar{\Sigma}$ captures the sum of variances of the terms containing
$\\{\overleftarrow{\mathbf{\epsilon}_{t}}\\}_{t>1}$. Recall that
${\bm{A}}^{T}{\bm{A}}={\bm{I}}$, the coefficients of drift
$\omega_{t}=\omega=\frac{1}{\bar{\nu}\sqrt{2}^{1/T}}$ and dispersion
$\xi_{t}=\sqrt{\frac{2^{(t-1)/T}}{2\beta_{t}(T-1)}}$. Now, let us compute the
variances of each term:

$\displaystyle\mathbb{E}\left[\left(\left(\omega_{1}\bar{\nu}{\bm{A}}{\bm{A}}^{T}\right)\xi_{2}\sqrt{\beta_{2}}\overleftarrow{\mathbf{\epsilon}_{2}}\right)\left(\left(\omega_{1}\bar{\nu}{\bm{A}}{\bm{A}}^{T}\right)\xi_{2}\sqrt{\beta_{2}}\overleftarrow{\mathbf{\epsilon}_{2}}\right)^{T}\right]=\omega_{1}^{2}\bar{\nu}^{2}\xi_{2}^{2}\beta_{2}{\bm{A}}{\bm{A}}^{T}\mathbb{E}\left[\overleftarrow{\mathbf{\epsilon}_{2}}\overleftarrow{\mathbf{\epsilon}_{2}}^{T}\right]{\bm{A}}{\bm{A}}^{T}=\omega_{1}^{2}\bar{\nu}^{2}\xi_{2}^{2}\beta_{2}{\bm{A}}{\bm{A}}^{T}{\bm{A}}{\bm{A}}^{T}$
$\displaystyle=\omega_{1}^{2}\bar{\nu}^{2}\xi_{2}^{2}\beta_{2}{\bm{A}}{\bm{A}}^{T}=\left(\frac{1}{\bar{\nu}\sqrt{2}^{1/T}}\right)^{2}\bar{\nu}^{2}\left(\frac{2^{1/T}}{2\beta_{2}(T-1)}\right)\beta_{2}{\bm{A}}{\bm{A}}^{T}=\frac{{\bm{A}}{\bm{A}}^{T}}{2(T-1)};$
$\displaystyle\mathbb{E}\left[\left(\left(\omega_{1}\bar{\nu}{\bm{A}}{\bm{A}}^{T}\right)\left(\omega_{2}\bar{\nu}{\bm{A}}{\bm{A}}^{T}\right)\xi_{3}\sqrt{\beta_{3}}\overleftarrow{\mathbf{\epsilon}_{3}}\right)\left(\left(\omega_{1}\bar{\nu}{\bm{A}}{\bm{A}}^{T}\right)\left(\omega_{2}\bar{\nu}{\bm{A}}{\bm{A}}^{T}\right)\xi_{3}\sqrt{\beta_{3}}\overleftarrow{\mathbf{\epsilon}_{3}}\right)^{T}\right]=\frac{{\bm{A}}{\bm{A}}^{T}}{2(T-1)};\cdots;$
$\displaystyle\mathbb{E}\left[\left(\left(\prod_{t=1}^{T-1}\left(\omega_{t}\bar{\nu}{\bm{A}}{\bm{A}}^{T}\right)\right)\xi_{T}\sqrt{\beta_{T}}\overleftarrow{\epsilon_{T}}\right)\left(\left(\prod_{t=1}^{T-1}\left(\omega_{t}\bar{\nu}{\bm{A}}{\bm{A}}^{T}\right)\right)\xi_{T}\sqrt{\beta_{T}}\overleftarrow{\epsilon_{T}}\right)^{T}\right]=\frac{{\bm{A}}{\bm{A}}^{T}}{2(T-1)}.$

Note that
$\left({\bm{A}}{\bm{A}}^{T}\right)^{2}\coloneqq{\bm{A}}{\bm{A}}^{T}{\bm{A}}{\bm{A}}^{T}={\bm{A}}{\bm{A}}^{T}$.
Since there are $(T-1)$ terms, the total covariance
$\bar{\Sigma}=\frac{{\bm{A}}{\bm{A}}^{T}}{2}$, which is equivalent to the
following sample generated at the end of the reverse process:

$\displaystyle\overleftarrow{{\bm{x}}_{0}}=\left(\prod_{t=1}^{T}\left(\omega_{t}\bar{\nu}{\bm{A}}{\bm{A}}^{T}\right)\right)\overleftarrow{{\bm{x}}_{T}}+\frac{{\bm{A}}{\bm{A}}^{T}}{\sqrt{2}}\overleftarrow{\mathbf{\epsilon}}={\bm{A}}{\bm{A}}^{T}\left(\prod_{t=1}^{T}\left(\frac{1}{\bar{\nu}\sqrt{2}^{1/T}}\right)\bar{\nu}\right)\overleftarrow{{\bm{x}}_{T}}+\frac{{\bm{A}}{\bm{A}}^{T}}{\sqrt{2}}\overleftarrow{\mathbf{\epsilon}}$
$\displaystyle=\frac{{\bm{A}}{\bm{A}}^{T}}{\sqrt{2}}\overleftarrow{{\bm{x}}_{T}}+\frac{{\bm{A}}{\bm{A}}^{T}}{\sqrt{2}}\overleftarrow{\mathbf{\epsilon}}={\bm{A}}\left(\frac{{\bm{A}}^{T}\overleftarrow{{\bm{x}}_{T}}}{\sqrt{2}}+\frac{{\bm{A}}^{T}\overleftarrow{\mathbf{\epsilon}}}{\sqrt{2}}\right)={\bm{A}}\overleftarrow{{\bm{z}}_{0}}.$

Now, let us verify if $\overleftarrow{{\bm{x}}_{0}}$ lies on the data
manifold. To see this, it is sufficient to show $\overleftarrow{{\bm{z}}_{0}}$
has zero mean and unit covariance. Therefore, we compute

$\displaystyle\mathbb{E}\left[\overleftarrow{{\bm{z}}_{0}}\right]=\mathbb{E}\left[\frac{{\bm{A}}^{T}\overleftarrow{{\bm{x}}_{T}}}{\sqrt{2}}+\frac{{\bm{A}}^{T}\overleftarrow{\mathbf{\epsilon}}}{\sqrt{2}}\right]=\frac{{\bm{A}}^{T}\mathbb{E}\left[\overleftarrow{{\bm{x}}_{T}}\right]}{\sqrt{2}}+\frac{{\bm{A}}^{T}\mathbb{E}\left[\overleftarrow{\mathbf{\epsilon}}\right]}{\sqrt{2}}=\mathbf{0}$

and

$\displaystyle\mathbb{E}\left[\overleftarrow{{\bm{z}}_{0}}\overleftarrow{{\bm{z}}_{0}}^{T}\right]=\mathbb{E}\left[\left(\frac{{\bm{A}}^{T}\overleftarrow{{\bm{x}}_{T}}}{\sqrt{2}}+\frac{{\bm{A}}^{T}\overleftarrow{\mathbf{\epsilon}}}{\sqrt{2}}\right)\left(\frac{{\bm{A}}^{T}\overleftarrow{{\bm{x}}_{T}}}{\sqrt{2}}+\frac{{\bm{A}}^{T}\overleftarrow{\mathbf{\epsilon}}}{\sqrt{2}}\right)^{T}\right]\stackrel{{\scriptstyle(i)}}{{=}}\frac{{\bm{A}}^{T}\mathbb{E}\left[\overleftarrow{{\bm{x}}_{T}}\overleftarrow{{\bm{x}}_{T}}^{T}\right]{\bm{A}}}{2}+\frac{{\bm{A}}^{T}\mathbb{E}\left[\overleftarrow{\mathbf{\epsilon}}\overleftarrow{\mathbf{\epsilon}}^{T}\right]{\bm{A}}}{2}\stackrel{{\scriptstyle(ii)}}{{=}}{\bm{I}},$

where (i) uses the fact that $\overleftarrow{{\bm{x}}_{T}}$ and
$\overleftarrow{\mathbf{\epsilon}}$ are independent zero-mean random
variables, and (ii) relies on unit covariance of
$\overleftarrow{{\bm{x}}_{T}}$, $\overleftarrow{\mathbf{\epsilon}}$, and the
property that ${\bm{A}}^{T}{\bm{A}}={\bm{I}}$. Thus, we finish the proof.
$\square$

###  A.8 Proof of Theorem 12

Besides the inpainting step at every intermediate state, the proof follows
from the arguments in Appendix A.7. Starting from the Gaussian prior
$\overleftarrow{{\bm{x}}_{T}}\sim{\mathcal{N}}\left(\mathbf{0},{\bm{I}}\right)$,
We first compute

$\displaystyle\overleftarrow{{\bm{x}}_{T-1}}=\omega_{T}\bar{\nu}{\bm{A}}{\bm{A}}^{T}\overleftarrow{{\bm{x}}_{T}}+\xi_{T}\sqrt{\beta_{T}}\overleftarrow{\mathbf{\epsilon}_{T}}.$

In the next step, we run the forward SDE to compute the known portions of the
given image with the noise level at $(T-1)$:

$\displaystyle\overrightarrow{{\bm{x}}_{T-1}}=\sqrt{\bar{\alpha}_{T-1}}\overrightarrow{{\bm{x}}_{0}}+\sqrt{1-\bar{\beta}_{T-1}}\overrightarrow{\mathbf{\epsilon}_{T}}.$

Then, we update $\overleftarrow{{\bm{x}}_{T-1}}$ by combining known portions
with the generated unknown portions of the image:

$\displaystyle\overleftarrow{{\bm{x}}_{T-1}}={\bm{D}}({\bm{m}})~{}\overleftarrow{{\bm{x}}_{T-1}}+{\bm{D}}(1-{\bm{m}})~{}\overrightarrow{{\bm{x}}_{T-1}}$

This completes inpainting at the state $(T-1)$ of the reverse process.
Recursively applying this over all the states, we obtain