# Improving Diffusion Models for Inverse Problems
Using Optimal Posterior Covariance

Xinyu Peng  Ziyang Zheng  Wenrui Dai  Nuoqian Xiao  Chenglin Li  Junni Zou
Hongkai Xiong

###### Abstract

Recent diffusion models provide a promising zero-shot solution to noisy linear
inverse problems without retraining for specific inverse problems. In this
paper, we propose the first unified interpretation for existing zero-shot
methods from the perspective of approximating the conditional posterior mean
for the reverse diffusion process of conditional sampling. We reveal that
recent methods are equivalent to making isotropic Gaussian approximations to
intractable posterior distributions over clean images given diffused noisy
images, with the only difference in the handcrafted design of isotropic
posterior covariances. Inspired by this finding, we propose a general plug-
and-play posterior covariance optimization based on maximum likelihood
estimation to improve recent methods. To achieve optimal posterior covariance
without retraining, we provide general solutions based on two approaches
specifically designed to leverage pre-trained models with and without reverse
covariances. Experimental results demonstrate that the proposed methods
significantly enhance the overall performance or robustness to hyperparameters
of recent methods. Code is available at
https://github.com/xypeng9903/k-diffusion-inverse-problems

Machine Learning, ICML

##  1 Introduction

Noisy linear inverse problems are widely studied for a variety of image
processing tasks, including denoising, inpainting, deblurring, and super-
resolution. The noisy linear inverse problems are formulated to accommodate
the widely adopted degradation model where images are measured with a linear
projection under noise corruption.

Recently, diffusion models have been emerging as promising methods for solving
inverse problems. According to the training strategies, these methods can be
categorized into two groups: 1) supervised methods that aim to learn a
conditional diffusion model using datasets consisting of pairs of degraded and
clean images (Saharia et al., 2022; Whang et al., 2022; Luo et al., 2023; Chan
et al., 2023), and 2) zero-shot methods that leverage pre-trained
unconditional diffusion models for conditional sampling in various scenarios
of inverse problems without the requirement of retraining. In this paper, we
focus on the zero-shot methods that can accommodate various tasks without
retraining. To ensure the consistency of optimizing fidelity under conditional
sampling, existing zero-shot methods adopt projection onto the measurement
subspace (Choi et al., 2021; Lugmayr et al., 2022; Song et al., 2022; Wang et
al., 2023; Zhu et al., 2023), leverage the similar idea as classifier guidance
(Song et al., 2021b; Dhariwal & Nichol, 2021) to modify the sampling process
with the likelihood score (Song et al., 2023; Chung et al., 2023), or resort
to variational inference (Feng et al., 2023; Mardani et al., 2023).

Diffusion models initially establish a forward process that introduces noise
to the original data ${\bm{x}}_{0}$ to generate noisy data ${\bm{x}}_{t}$ at
time $t$, and then implement the reverse process to generate ${\bm{x}}_{0}$
obeying the distribution of original data. The key to realize the reverse
process is the posterior mean
$\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]=\int{\bm{x}}_{0}p_{t}({\bm{x}}_{0}|{\bm{x}}_{t})\mathrm{d}{\bm{x}}_{0}$
that could be the optimal estimation of ${\bm{x}}_{0}$ given ${\bm{x}}_{t}$ in
the sense of minimum mean square error (MMSE). In the reverse process for
conditional sampling in inverse problems, estimating ${\bm{x}}_{0}$ is based
on both ${\bm{x}}_{t}$ and measurement ${\bm{y}}$. Thus, the optimal
estimation is determined by the conditional posterior mean
$\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}}]$.

In this paper, we provide the first unified interpretation for zero-shot
methods from the view of approximating the conditional posterior mean for the
reverse diffusion process of conditional sampling. We interpret recent zero-
shot methods (Song et al., 2023; Chung et al., 2023; Wang et al., 2023; Zhu et
al., 2023) with two types of guidance that uniformly employ an isotropic
Gaussian approximation for the intractable posterior distribution
$p_{t}({\bm{x}}_{0}|{\bm{x}}_{t})$ to approximate the conditional posterior
mean. Motivated by the insight, we reveal their limitation in handcrafted
design of posterior covariance and propose a generalized method for improving
recent methods that optimize the posterior covariance based on maximum
likelihood estimation. We achieve optimal posterior covariance using pre-
trained unconditional diffusion models in a plug-and-play fashion by
converting (available) reverse covariances or via Monte Carlo estimation
without reverse covariances. Remarkably, this is the first attempt that
successfully predicts posterior covariance in addition to posterior mean from
a general class of existing unconditional diffusion models. Experimental
results demonstrate that the proposed method significantly outperforms
existing zero-shot methods in a wide range of tasks like inpainting,
deblurring, and super-resolution, and enhances the robustness to
hyperparameters.

##  2 Background

###  2.1 Bayesian Framework for Solving Inverse Problems

In inverse problems, we wish to recover the unknown signal
${\bm{x}}_{0}\in\mathbb{R}^{d}$ from noisy measurements
${\bm{y}}\in\mathbb{R}^{m}$, where

${\bm{y}}={\bm{A}}{\bm{x}}_{0}+{\bm{n}}$ (1)

It is usually assumed that by the expert knowledge of the sensing device,
${\bm{A}}\in\mathbb{R}^{m\times d}$ is known and
${\bm{n}}\sim\mathcal{N}(\bm{0},\sigma^{2}{\bm{I}})$ is an i.i.d. additive
Gaussian noise with a known standard deviation of $\sigma$. However,
recovering ${\bm{x}}_{0}$ from the degraded measurements ${\bm{y}}$ becomes
ill-posed due to the under-determined ${\bm{A}}$ and the existence of noise.
To solve the ill-posed inverse problem, a prior is required to constrain the
solution according to desired image statistics (Feng et al., 2023). From a
Bayesian perspective, the prior is described by assuming ${\bm{x}}_{0}$ obeys
an unknown prior distribution $p({\bm{x}}_{0})$, and we have a likelihood
function
$p({\bm{y}}|{\bm{x}}_{0})=\mathcal{N}({\bm{y}}|{\bm{A}}{\bm{x}}_{0},\sigma^{2}{\bm{I}})$
determined by equation 1. We solve inverse problems by formulating a posterior
distribution $p({\bm{x}}_{0}|{\bm{y}})$ over clean images given noisy
measurements, which in principle is determined given $p({\bm{x}}_{0})$ and
$p({\bm{y}}|{\bm{x}}_{0})$ by Bayes’ theorem:
$p({\bm{x}}_{0}|{\bm{y}})=p({\bm{x}}_{0})p({\bm{y}}|{\bm{x}}_{0})/\int
p({\bm{x}}_{0})p({\bm{y}}|{\bm{x}}_{0})\mathrm{d}{\bm{x}}_{0}$.

###  2.2 Diffusion Models and Conditioning

We use diffusion models to model the complex posterior distribution
$p({\bm{x}}_{0}|{\bm{y}})$ for inverse problems. Let us define a family of
Gaussian perturbation kernels $p_{t}({\bm{x}}_{t}|{\bm{x}}_{0})$ of
${\bm{x}}_{0}\sim p({\bm{x}}_{0})$ by injecting _i.i.d._ Gaussian noise of
standard deviation $\sigma_{t}$ to ${\bm{x}}_{0}$ and scaling by the factor of
$s_{t}$, _i.e._ ,
$p_{t}({\bm{x}}_{t}|{\bm{x}}_{0})=\mathcal{N}({\bm{x}}_{t}|s_{t}{\bm{x}}_{0},s_{t}^{2}\sigma_{t}^{2}{\bm{I}})$,
where $\sigma_{t}$ is monotonically increasing with respect to time
$t=0,1,\cdots,T$. We start from $\sigma_{0}=0$ and reach a value of
$\sigma_{T}$ that is much larger than the standard deviation of
$p({\bm{x}}_{0})$, to ensure that samples from ${\bm{x}}_{T}\sim
p({\bm{x}}_{T})$ are indistinguishable to samples from
$\mathcal{N}(\bm{0},s_{T}^{2}\sigma_{T}^{2}{\bm{I}})$. Since ${\bm{x}}_{t}$ is
independent of ${\bm{y}}$ once ${\bm{x}}_{0}$ is known, we characterize the
joint distribution between ${\bm{x}}_{0},{\bm{y}}$, and ${\bm{x}}_{t}$ as
$p_{t}({\bm{x}}_{0},{\bm{y}},{\bm{x}}_{t})=p({\bm{x}}_{0})p({\bm{y}}|{\bm{x}}_{0})p_{t}({\bm{x}}_{t}|{\bm{x}}_{0})$,
which can also be represented by a probabilistic graphical model
${\bm{y}}\leftarrow{\bm{x}}_{0}\rightarrow{\bm{x}}_{t}$. There exist multiple
formulations of diffusion models in the literature (Song & Ermon, 2019; Ho et
al., 2020; Song et al., 2021a; Kingma et al., 2021; Song et al., 2021b; Karras
et al., 2022). Here, we use the ordinary differential equation (ODE)
formulation and select $s_{t}=1,\sigma_{t}=t$ as suggested in (Karras et al.,
2022) for simplicity. Let us consider the following ODE:

$\mathrm{d}{\bm{x}}_{t}=\frac{{\bm{x}}_{t}-\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]}{t}\mathrm{d}t,\
\ {\bm{x}}_{T}\sim p_{T}({\bm{x}}_{T}).$ (2)

In equation 2, the only source of randomness is the initial sample
${\bm{x}}_{T}\sim p_{T}({\bm{x}}_{T})$. The ODE possesses an important
property that ${\bm{x}}_{t}$ generated by the ODE maintains the exact same
marginals to ${\bm{x}}_{t}$ obtained by injecting Gaussian noise to
${\bm{x}}_{0}$, _i.e._ , $p_{t}({\bm{x}}_{t})$. Generally, diffusion models
approximate $\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]$ with a time-dependent
denoiser $D_{t}({\bm{x}}_{t})$ (referred to as the unconditional diffusion
model throughout the paper) that is trained by minimizing the $L_{2}$ loss for
all $t\in[0,T]$:

$\min_{D_{t}}\mathbb{E}_{p_{t}({\bm{x}}_{0},{\bm{x}}_{t})}\left[\lVert{\bm{x}}_{0}-D_{t}({\bm{x}}_{t})\rVert_{2}^{2}\right].$
(3)

With sufficient data and model capacity, the optimal $D_{t}({\bm{x}}_{t})$ is
the MMSE estimator of ${\bm{x}}_{0}$ given ${\bm{x}}_{t}$ and equals to
$\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]$. Thus, samples from $p({\bm{x}}_{0})$
can be obtained by first sampling ${\bm{x}}_{T}$ from
$\mathcal{N}(\bm{0},s_{T}^{2}\sigma_{T}^{2}{\bm{I}})$ and then simulating
equation 2 from $t=T$ to $t=0$ using black box ODE solver that replaces
$\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]$ with a well-trained
$D_{t}({\bm{x}}_{t})$.

To solve inverse problems, we are interested in $p({\bm{x}}_{0}|{\bm{y}})$. We
condition the diffusion models and formulate an ODE whose marginals are
$p_{t}({\bm{x}}_{t}|{\bm{y}})$ such that we can simulate the ODE to sample
from $p({\bm{x}}_{0}|{\bm{y}})$. The desired ODE is

$\mathrm{d}{\bm{x}}_{t}=\frac{{\bm{x}}_{t}-\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}}]}{t}\mathrm{d}t,\
\ {\bm{x}}_{T}\sim p_{T}({\bm{x}}_{T}|{\bm{y}}).$ (4)

For sufficiently large $\sigma_{T}$, samples ${\bm{x}}_{T}\sim
p({\bm{x}}_{T}|{\bm{y}})$ are indistinguishable to samples from
$\mathcal{N}(\bm{0},s_{T}^{2}\sigma_{T}^{2}{\bm{I}})$ (see Appendix D.2 in
(Dhariwal & Nichol, 2021)). The unconditioned sampling procedure in equation 2
can be used to achieve conditioned sampling in equation 4, by substituting
$\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]$ with the conditional posterior mean
$\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}}]$111Note that this conclusion
can be extended to DDPM, DDIM, and SDE-formalized diffusion models..

##  3 Unified Interpretation of Diffusion-based Solvers to Inverse Problems

To simulate equation 4 for solving inverse problems, a standard approach is
training a conditional diffusion model for approximating
$\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}}]$ using supervised learning.
However, this approach can be computationally demanding by training separate
models for different scenarios. In this paper, we aim to leverage pre-trained
unconditional diffusion models for conditional sampling in inverse problems to
avoid additional training. We summarize three steps below:
1) Obtain the estimation of $\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]$ through
the unconditional diffusion model $D_{t}({\bm{x}}_{t})$;
2) Obtain the estimation of the conditional posterior mean
$\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}}]$ based on the estimation of
$\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]$;
3) Substitute the estimation of
$\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}}]$ for
$\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]$ in the sampling process (see equation
2 and equation 4).

Recent zero-shot methods, including DPS (Chung et al., 2023), $\Pi\text{GDM}$
(Song et al., 2023), DDNM (Wang et al., 2023), and DiffPIR (Zhu et al., 2023),
explicitly or implicitly follow the three steps but only differ in the
estimation of $\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}}]$ from
$\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]$ in Step 2). Specifically, these
methods can be unified with the interpretation of estimating
$\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}}]$ with isotropic Gaussian
approximations
$\mathcal{N}({\bm{x}}_{0}|\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}],r_{t}^{2}{\bm{I}})$
to the intractable posterior $p_{t}({\bm{x}}_{0}|{\bm{x}}_{t})$ with different
$r_{t}$, as summarized in Table 1. We classify them into two categories,
_i.e._ , Type I and Type II guidance, according to approximation paradigms, as
elaborated below.

###  3.1 Type I Guidance: Approximating the Likelihood Score Function

We classify DPS and $\Pi\text{GDM}$ into one category, referred to as Type I
guidance. Type I guidance is closely related to classifier guidance (Dhariwal
& Nichol, 2021), where the conditional posterior mean
$\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}}]$ is approximated based on the
following proposition.

######  Proposition 1.

The conditional posterior mean equals to the posterior mean drifted by scaled
likelihood score function. Formally,

$\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}}]=\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]+s_{t}\sigma^{2}_{t}\nabla_{{\bm{x}}_{t}}\log
p_{t}({\bm{y}}|{\bm{x}}_{t}).$ (5)

###### Proof.

Please refer to Appendix A.2. ∎

According to Proposition 1222In addition to realizing guidance for diffusion
models, Proposition 1 can be used for refining the estimation of
${\bm{x}}_{0}$ for learning sampling patterns in MRI (Ravula et al., 2023) and
realizing guidance for flow-based generative models to solve inverse problems
(Pokle et al., 2023)., we can approximate
$\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]$ using the pre-trained unconditional
diffusion model $D_{t}({\bm{x}}_{t})$; however, the likelihood score
$\nabla_{{\bm{x}}_{t}}\log p_{t}({\bm{y}}|{\bm{x}}_{t})$ is computationally
intractable. In inverse problems, only the likelihood
$p({\bm{y}}|{\bm{x}}_{0})$ at $t=0$ is known, and
$p_{t}({\bm{y}}|{\bm{x}}_{t})$ for any $t>0$ is obtained by an intractable
integral over all possible ${\bm{x}}_{0}$ as

$p_{t}({\bm{y}}|{\bm{x}}_{t})=\int
p({\bm{y}}|{\bm{x}}_{0})p_{t}({\bm{x}}_{0}|{\bm{x}}_{t})\mathrm{d}{\bm{x}}_{0}.$
(6)

Thus, Type I guidance like DPS (Chung et al., 2023) and $\Pi\text{GDM}$ (Song
et al., 2023) in nature approximately computes $p_{t}({\bm{y}}|{\bm{x}}_{t})$
based on isotropic Gaussian approximations to the true posterior
$p_{t}({\bm{x}}_{0}|{\bm{x}}_{t})$.

Methods | Guidance | $r_{t}$
---|---|---
DPS (Chung et al., 2023) | I | $\text{approach}~{}0$
$\Pi\text{GDM}$ (Song et al., 2023) | I | $\sqrt{\sigma_{t}^{2}/(\sigma_{t}^{2}+1)}$
DDNM (Wang et al., 2023) | II | any fixed value
DiffPIR (Zhu et al., 2023) | II | $\sigma_{t}/\sqrt{\lambda}$
Table 1: Unified interpretation of diffusion-based solvers to inverse
problems. Recent methods can be regarded as making isotropic Gaussian
approximations to the posterior distributions over clean images given noisy
images.

DPS (Chung et al., 2023) DPS can be viewed as approximating
$p_{t}({\bm{x}}_{0}|{\bm{x}}_{t})$ using a delta distribution
$\delta({\bm{x}}_{0}-\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}])$ centered at the
posterior mean $\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]$, which can be regarded
as the limit of the Gaussian
$\mathcal{N}({\bm{x}}_{0}|\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}],r_{t}^{2}{\bm{I}})$
when the variance $r_{t}^{2}$ approaches zero. The likelihood
$p_{t}({\bm{y}}|{\bm{x}}_{t})$ is approximated by

$\displaystyle p_{t}({\bm{y}}|{\bm{x}}_{t})$ $\displaystyle\approx\int
p({\bm{y}}|{\bm{x}}_{0})\delta({\bm{x}}_{0}-\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}])\mathrm{d}{\bm{x}}_{0}$
$\displaystyle=p({\bm{y}}|{\bm{x}}_{0}=\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]).$
(7)

However, $p({\bm{y}}|{\bm{x}}_{0}=\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}])$ is
not a good approximation to compute the likelihood score
$\nabla_{{\bm{x}}_{t}}\log p_{t}({\bm{y}}|{\bm{x}}_{t})$ in practice (Chung et
al., 2023). To address this issue, DPS empirically adjusts the strength of
guidance by replacing the likelihood score $\nabla_{{\bm{x}}_{t}}\log
p_{t}({\bm{y}}|{\bm{x}}_{t})$ with
$-\zeta_{t}\nabla_{{\bm{x}}_{t}}\lVert{\bm{y}}-{\bm{A}}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]\rVert_{2}^{2}$,
where
$\zeta_{t}=\zeta/\lVert{\bm{y}}-{\bm{A}}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]\rVert$
with a hyper-parameter $\zeta$.

$\Pi\text{GDM}$ (Song et al., 2023) The delta distribution used in DPS is a
very rough approximation to $p_{t}({\bm{x}}_{0}|{\bm{x}}_{t})$ as it
completely ignores the uncertainty of ${\bm{x}}_{0}$ given ${\bm{x}}_{t}$. As
$t$ increases, the uncertainty in $p_{t}({\bm{x}}_{0}|{\bm{x}}_{t})$ becomes
larger and is close to the original data distribution $p({\bm{x}}_{0})$. Thus,
it is more reasonable to choose a positive $r_{t}$. In $\Pi\text{GDM}$,
$r_{t}$ is heuristically selected as
$\sqrt{\sigma_{t}^{2}/(1+\sigma_{t}^{2})}$ under the assumption that
$p({\bm{x}}_{0})$ is the standard normal distribution
$\mathcal{N}(\bm{0},{\bm{I}})$. In such case, the likelihood
$p_{t}({\bm{y}}|{\bm{x}}_{t})$ is approximated by

$\displaystyle p_{t}({\bm{y}}|{\bm{x}}_{t})$
$\displaystyle\approx\int\mathcal{N}({\bm{y}}|{\bm{A}}{\bm{x}}_{0},\sigma^{2}{\bm{I}})\mathcal{N}({\bm{x}}_{0}|\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}],r_{t}^{2}{\bm{I}})\mathrm{d}{\bm{x}}_{0}$
$\displaystyle=\mathcal{N}({\bm{y}}|{\bm{A}}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}],\sigma^{2}{\bm{I}}+r_{t}^{2}{\bm{A}}{\bm{A}}^{T}).$
(8)

According to equation 8, the likelihood score is a vector-Jacobian product
that can be computed using back-propagation (Song et al., 2023).

###  3.2 Type II Guidance: Approximating the Conditional Posterior Mean Using
Proximal Solution

We classify DiffPIR and DDNM into the category of Type II guidance. Type II
guidance approximates the conditional posterior mean
$\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}}]$ by finding a solution
$\hat{{\bm{x}}}_{0}^{(t)}$ that is close to
$\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]$ and consistent with the measurement
${\bm{y}}$ at the same time. Compared to Type I guidance, Type II guidance is
more efficient by circumventing the expensive back-propagation of the
unconditional diffusion model $D_{t}({\bm{x}}_{t})$ for evaluating gradients.
Although Type II guidance seems to be developed from a different perspective,
we show that it can also be regarded as an isotropic Gaussian approximation to
the true posterior $p_{t}({\bm{x}}_{0}|{\bm{x}}_{t})$.

DiffPIR (Zhu et al., 2023) The core step in DiffPIR is replacing
$\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]$ in the unconditional sampling process
with the solution to the following proximal problem:

$\hat{{\bm{x}}}_{0}^{(t)}=\arg\min_{{\bm{x}}_{0}}\lVert{\bm{y}}-{\bm{A}}{\bm{x}}_{0}\rVert^{2}+\rho_{t}\lVert{\bm{x}}_{0}-\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]\rVert^{2},$
(9)

where $\rho_{t}=\lambda\sigma^{2}/\sigma_{t}^{2}$ and $\lambda$ is a hyper-
parameter. Our key insight is that equation 9 for DiffPIR can be interpreted
as approximating $\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}}]$ with the
mean of an approximate distribution
$q_{t}({\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}})$ derived from the isotropic
Gaussian approximation $q_{t}({\bm{x}}_{0}|{\bm{x}}_{t})$ to the intractable
posterior distribution $p_{t}({\bm{x}}_{0}|{\bm{x}}_{t})$. Let us define the
isotropic Gaussian approximation
$q_{t}({\bm{x}}_{0}|{\bm{x}}_{t})=\mathcal{N}({\bm{x}}_{0}|\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}],r_{t}^{2}{\bm{I}})$.
We can obtain that the approximate distribution
$q_{t}({\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}})\propto
p({\bm{y}}|{\bm{x}}_{0})q_{t}({\bm{x}}_{0}|{\bm{x}}_{t})$ for
$p_{t}({\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}})\propto
p({\bm{y}}|{\bm{x}}_{0})p_{t}({\bm{x}}_{0}|{\bm{x}}_{t})$333$p_{t}({\bm{y}}|{\bm{x}}_{0},{\bm{x}}_{t})=p({\bm{y}}|{\bm{x}}_{0})$
since we leverage the conditional independence between ${\bm{y}}$ and
${\bm{x}}_{t}$ given ${\bm{x}}_{0}$. is also a Gaussian and its mean
$\mathbb{E}_{q}[{\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}}]$ can be obtained by
solving the optimization problem:

$\displaystyle\mathbb{E}_{q}[{\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}}]$
$\displaystyle\ =\arg\max_{{\bm{x}}_{0}}\log
q_{t}({\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}})$ $\displaystyle\
=\arg\max_{{\bm{x}}_{0}}\left[\log p({\bm{y}}|{\bm{x}}_{0})\\!+\\!\log
q_{t}({\bm{x}}_{0}|{\bm{x}}_{t})\right]$ $\displaystyle\
=\arg\min_{{\bm{x}}_{0}}\left[\lVert{\bm{y}}\\!-\\!{\bm{A}}{\bm{x}}_{0}\rVert^{2}\\!+\\!\frac{\sigma^{2}}{r_{t}^{2}}\lVert{\bm{x}}_{0}\\!-\\!\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]\rVert^{2}\right].$
(10)

Note that equation 9 used in DiffPIR is a special case of equation 3.2 with
$r_{t}=\sigma_{t}/\sqrt{\lambda}$. Therefore, DiffPIR can be viewed as using
isotropic Gaussian approximation
$q_{t}({\bm{x}}_{0}|{\bm{x}}_{t})=\mathcal{N}({\bm{x}}_{0}|\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}],(\sigma_{t}^{2}/\lambda){\bm{I}})$
for $p_{t}({\bm{x}}_{0}|{\bm{x}}_{t})$.

DDNM (Wang et al., 2023) The core step in DDNM resorts to the range-null space
decomposition of $\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]$. Specifically, to
ensure the measurement consistency ${\bm{y}}={\bm{A}}\hat{{\bm{x}}}_{0}^{(t)}$
of the solution $\hat{{\bm{x}}}_{0}^{(t)}$, DDNM formulates
$\hat{{\bm{x}}}_{0}^{(t)}={\bm{A}}^{\dagger}{\bm{y}}+({\bm{I}}-{\bm{A}}^{\dagger}{\bm{A}})\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]$
by replacing the range space component of
$\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]$ with ${\bm{A}}^{\dagger}{\bm{y}}$ but
keeping the null space component unchanged. We find that, when the measurement
noise $\sigma$ vanishes, replacing the range space component
${\bm{A}}^{\dagger}{\bm{A}}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]$ with
${\bm{A}}^{\dagger}{\bm{y}}$ is equivalent to replacing
$\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]$ with
$\mathbb{E}_{q}[{\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}}]$. Therefore, DDNM is
equivalent to approximating $p_{t}({\bm{x}}_{0}|{\bm{x}}_{t})$ using an
isotropic Gaussian, as formalized in Proposition 2.

######  Proposition 2.

For any $r_{t}>0$, $\mathbb{E}_{q}[{\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}}]$
approaches $\hat{{\bm{x}}}_{0}^{(t)}$ used in DDNM as the variance of
measurement noise approaches zero. Formally, we have

$\lim_{\sigma\rightarrow
0}\mathbb{E}_{q}[{\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}}]={\bm{A}}^{\dagger}{\bm{y}}+({\bm{I}}-{\bm{A}}^{\dagger}{\bm{A}})\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}].$
(11)

###### Proof.

Please refer to Appendix A.3. ∎

###  3.3 Solving Inverse Problems With Optimal Posterior Covariance

In Sections 3.1 and 3.2, we unify existing zero-shot diffusion-based solvers
with the interpretation of isotropic Gaussian approximations to the
intractable posterior distributions $p_{t}({\bm{x}}_{0}|{\bm{x}}_{t})$. This
unified interpretation motivates us to further optimize the posterior
covariance in the Gaussian approximation for enhancing existing diffusion-
based solvers. Specifically, we consider to approximate
$p_{t}({\bm{x}}_{0}|{\bm{x}}_{t})$ using variational Gaussian
$q_{t}({\bm{x}}_{0}|{\bm{x}}_{t})=\mathcal{N}({\bm{x}}_{0}|{\bm{\mu}}_{t}({\bm{x}}_{t}),{\bm{\Sigma}}_{t}({\bm{x}}_{t}))$
with learnable covariance ${\bm{\Sigma}}_{t}({\bm{x}}_{t})$. In the pre-
training stage, $q_{t}({\bm{x}}_{0}|{\bm{x}}_{t})$ is learned by minimizing
the weighted integral of expected forward KL divergence between
$p_{t}({\bm{x}}_{0}|{\bm{x}}_{t})$ and $q_{t}({\bm{x}}_{0}|{\bm{x}}_{t})$ as

$\min_{q}\int\omega_{t}\mathbb{E}_{p_{t}({\bm{x}}_{t})}[D_{KL}(p_{t}({\bm{x}}_{0}|{\bm{x}}_{t})\|q_{t}({\bm{x}}_{0}|{\bm{x}}_{t}))]\mathrm{d}t.$
(12)

Equation 12 can be achieved in a tractable way by maximizing the log-
likelihood of $q_{t}({\bm{x}}_{0}|{\bm{x}}_{t})$. When
${\bm{\Sigma}}_{t}({\bm{x}}_{t})$ is full rank, the optimal
${\bm{\mu}}_{t}({\bm{x}}_{t})$ is exactly the MMSE estimator
$\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]$. Once
$q_{t}({\bm{x}}_{0}|{\bm{x}}_{t})$ is learned, we can leverage
${\bm{\mu}}_{t}$ and ${\bm{\Sigma}}_{t}$ to solve inverse problems and
accommodate to existing methods by plugging ${\bm{\mu}}_{t}({\bm{x}}_{t})$ in
equation 2 to achieve unconditional sampling. The solutions to Type I and Type
II guidance are summarized as below.

Type I guidance. The likelihood is approximated in a similar way to equation 8
using ${\bm{\mu}}_{t}$ and ${\bm{\Sigma}}_{t}$:

$p_{t}({\bm{y}}|{\bm{x}}_{t})\approx\mathcal{N}({\bm{y}}|{\bm{A}}{\bm{\mu}}_{t}({\bm{x}}_{t}),\sigma^{2}{\bm{I}}+{\bm{A}}{\bm{\Sigma}}_{t}({\bm{x}}_{t}){\bm{A}}^{T}).$
(13)

Type II guidance. $\hat{{\bm{x}}}_{0}^{(t)}$ is solved from the auto-weighted
proximal problem associated with ${\bm{\mu}}_{t}$ and ${\bm{\Sigma}}_{t}$:

$\hat{{\bm{x}}}_{0}^{(t)}=\arg\min_{{\bm{x}}_{0}}\lVert{\bm{y}}-{\bm{A}}{\bm{x}}_{0}\rVert^{2}+\sigma^{2}\lVert{\bm{x}}_{0}-{\bm{\mu}}_{t}({\bm{x}}_{t})\rVert^{2}_{{\bm{\Sigma}}_{t}^{-1}({\bm{x}}_{t})},$
(14)

where
$\lVert{\bm{x}}\rVert_{{\bm{\Lambda}}}^{2}={\bm{x}}^{T}{\bm{\Lambda}}{\bm{x}}$
with positive definite matrix ${\bm{\Lambda}}$.

To avoid the inversion of high-dimensional matrices for realizing guidances in
practice, we develop closed-form solutions for Type I guidance in equation 13
and Type II guidance in equation 14 under the isotropic posterior covariance
${\bm{\Sigma}}_{t}({\bm{x}}_{t})=r_{t}^{2}({\bm{x}}_{t}){\bm{I}}$, as well as
numerical solutions based on conjugate gradient method (CG) for more general
covariance. Please refer to Appendices B and C for details.

##  4 Plug-and-Play Posterior Covariance Optimization

Since training diffusion models from scratch using equation 12 is expensive,
directly leveraging available checkpoints of unconditional diffusion models is
prevailing in prior works. To achieve plug-and-play posterior covariance
optimization without retraining, it is non-trivial to obtain posterior
covariance prediction from unconditional diffusion models, which has yet to be
explored in existing methods. In this section, we develop two alternative ways
to achieve optimal diagonal posterior covariance without retraining for two
common cases: 1) reverse covariance prediction is available from the given
unconditional diffusion model (Section 4.1), and 2) reverse covariance
prediction is not available (Section 4.2). We also discuss the relation
between this paper and the concurrent work (Boys et al., 2023). It is worth
noting that, this paper focuses on optimizing diagonal posterior covariance to
accommodate existing unconditional diffusion models. Nevertheless, we
demonstrate in Appendix D.5 that the proposed method has the potential for
generalizing to correlated (non-diagonal) posterior covariance. We will
further explore this problem in the future.

###  4.1 Converting Optimal Reverse Variances

Recent pre-trained diffusion models often predict the optimal reverse
variances444Variances are the diagonal elements of the covariance matrix.
(Nichol & Dhariwal, 2021) for improving performance when using ancestral
sampling of discrete diffusion models (or denoising diffusion probabilistic
models (DDPM) (Ho et al., 2020; Sohl-Dickstein et al., 2015). Here, we develop
a solution that can directly leverage off-the-shelf reverse covariance
prediction based on the fixed-point solution of DDPM to obtain posterior
covariance prediction solved from equation 12 and avoid retraining. We first
develop the fixed-point solutions to equation 12 under diagonal posterior
covariance case in Proposition 3, and then establish the relation between the
fixed-point solutions to equation 12 and DDPM in equation 19 in Theorem 4.

######  Proposition 3 (Fixed-point solutions of variational Gaussian
posterior).

The optimal mean ${\bm{\mu}}_{t}^{*}({\bm{x}}_{t})$ and the optimal diagnol
posterior covariances
${\bm{\Sigma}}_{t}^{*}({\bm{x}}_{t})=\mathrm{diag}[{\bm{r}}_{t}^{*2}({\bm{x}}_{t})]$
to equation 12 are obtained by

$\displaystyle{\bm{\mu}}_{t}^{*}({\bm{x}}_{t})=\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}],$
(15)
$\displaystyle{\bm{r}}_{t}^{*2}({\bm{x}}_{t})=\mathbb{E}_{p_{t}({\bm{x}}_{0}|{\bm{x}}_{t})}[({\bm{x}}_{0}-\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}])^{2}].$
(16)

###### Proof.

Please refer to Appendix A.5. ∎

Subsequently, we consider the fixed-point solutions of DDPM. Unlike continuous
ODE formulation introduced in Section 2.2, diffusion models pre-trained under
the DDPM framework are latent variable models defined by $p({\bm{x}}_{0})=\int
p({\bm{x}}_{0:T})\mathrm{d}{\bm{x}}_{1:T}$. The joint distribution
$p({\bm{x}}_{0:T})$ is referred to as the reverse process defined as a Markov
chain of learnable Gaussian transitions starting at
$p({\bm{x}}_{T})=\mathcal{N}(\bm{0},{\bm{I}})$ (Ho et al., 2020). The Gaussian
transition at time $t$ is characterized by the mean ${\bm{m}}_{t}$ and
covariance ${\bm{C}}_{t}$ dependent on the condition ${\bm{x}}_{t}$:

$\displaystyle
p({\bm{x}}_{0:T})=p({\bm{x}}_{T})\prod_{t=1}^{T}p({\bm{x}}_{t-1}|{\bm{x}}_{t}),$
(17) $\displaystyle
p({\bm{x}}_{t-1}|{\bm{x}}_{t})=\mathcal{N}({\bm{x}}_{t-1}|{\bm{m}}_{t}({\bm{x}}_{t}),{\bm{C}}_{t}({\bm{x}}_{t})).$
(18)

DDPM defines a forward process $q({\bm{x}}_{1:T}|{\bm{x}}_{0})$ by gradually
injecting noise to the data. Please refer to Appendix A.4 for the detailed
definition of $q({\bm{x}}_{1:T}|{\bm{x}}_{0})$. We fit $p({\bm{x}}_{0})$ to
the original data distribution $q({\bm{x}}_{0})$ by minimizing the KL
divergence between the forward and reverse processes:

$\min_{p}D_{KL}(q({\bm{x}}_{0:T})||p({\bm{x}}_{0:T})),$ (19)

where $q({\bm{x}}_{0:T})=q({\bm{x}}_{0})q({\bm{x}}_{1:T}|{\bm{x}}_{0})$.

In Theorem 4, we present the fixed-point solution of DDPM under signal-
dependent diagonal covariance
${\bm{C}}_{t}({\bm{x}}_{t})=\mathrm{diag}[{\bm{v}}_{t}^{2}({\bm{x}}_{t})]$ as
used in (Nichol & Dhariwal, 2021).

######  Theorem 4 (Fixed-point solutions of DDPM).

Let ${\bm{C}}_{t}({\bm{x}}_{t})=\mathrm{diag}[{\bm{v}}_{t}^{2}({\bm{x}}_{t})]$
be a signal-dependent diagonal covariance for the reverse covariance. When
$\tilde{{\bm{\mu}}}_{t},\beta_{t},\bar{\alpha}_{t},\bar{\beta}_{t}$ are
determined by the forward process $q({\bm{x}}_{1:T}|{\bm{x}}_{0})$, the
optimal solutions ${\bm{m}}^{*}_{t}({\bm{x}}_{t})$ and
${\bm{v}}_{t}^{*2}({\bm{x}}_{t})$ to equation 19 are

$\displaystyle{\bm{m}}^{*}_{t}({\bm{x}}_{t})=\tilde{{\bm{\mu}}}_{t}({\bm{x}}_{t},\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]),$
(20)
$\displaystyle{\bm{v}}_{t}^{*2}({\bm{x}}_{t})=\tilde{\beta}_{t}+(\frac{\sqrt{\bar{\alpha}_{t-1}}\beta_{t}}{1-\bar{\alpha}_{t}})^{2}\cdot{\bm{r}}_{t}^{*2}({\bm{x}}_{t})$
(21)

where ${\bm{r}}_{t}^{*2}({\bm{x}}_{t})$ is the optimal posterior variances
determined by equation 16 under $s_{t}=\sqrt{\bar{\alpha}_{t}}$ and
$\sigma_{t}=\sqrt{\bar{\beta}_{t}/\bar{\alpha}_{t}}$, and
$\tilde{\beta}_{t}=(\bar{\beta}_{t-1}/\bar{\beta}_{t})\beta_{t}$.

###### Proof.

Please refer to Appendix A.4. ∎

Theorem 4 implies that, given the reverse variances
$\hat{{\bm{v}}}_{t}^{2}({\bm{x}}_{t})$ predicted by a pre-trained DDPM model
at time $t$, the posterior variances are obtained according to equation 21 by

$\hat{{\bm{r}}}_{t}^{2}({\bm{x}}_{t})=(\hat{{\bm{v}}}_{t}^{2}({\bm{x}}_{t})-\tilde{\beta}_{t})\cdot\left(\frac{\sqrt{\bar{\alpha}_{t-1}}\beta_{t}}{1-\bar{\alpha}_{t}}\right)^{-2}.$
(22)

In fact, similar fixed-point solutions of DDPM have been proposed in (Bao et
al., 2022b) and (Bao et al., 2022a) for determining optimal reverse variances.
Our approach can be viewed as opposite direction to (Bao et al., 2022b) and
(Bao et al., 2022a) for applying the fixed-point solutions.

###  4.2 Monte Carlo Estimation of Posterior Variances

We further develop the optimal posterior variances using a pre-trained
unconditional diffusion model without providing the reverse variances
prediction $\hat{{\bm{v}}}_{t}^{2}({\bm{x}}_{t})$ (_e.g._ , (Ho et al.,
2020)). In this case, we suppose signal-independent
${\bm{\Sigma}}_{t}({\bm{x}}_{t})=r_{t}^{2}{\bm{I}}$ is an isotropic covariance
with a time-dependent standard deviation $r_{t}$. We directly optimize $r_{t}$
by forcing the derivative of equation 12 with an optimal ${\bm{\mu}}_{t}$ with
regard to $r_{t}$ to zero. Thus, we obtain that

$r_{t}^{*2}=\frac{1}{d}\mathbb{E}_{p_{t}({\bm{x}}_{0},{\bm{x}}_{t})}[\lVert{\bm{x}}_{0}-\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]\rVert^{2}].$
(23)

Note that we only estimate the expected reconstruction error of the MMSE
estimator $\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]$. Similar to (Bao et al.,
2022b), given a pre-trained unconditional diffusion model $D_{t}$,
$r_{t}^{*2}$ is estimated using Monte Carlo samples as

$\hat{r}_{t}^{2}=\frac{1}{dM}\sum_{m=1}^{M}\lVert{\bm{x}}_{0}^{(m)}-D_{t}({\bm{x}}_{t}^{(m)})\rVert^{2}_{2},$
(24)

where ${\bm{x}}_{0}^{(m)},{\bm{x}}_{t}^{(m)}\sim
p({\bm{x}}_{0})p_{t}({\bm{x}}_{t}|{\bm{x}}_{0})$ and $M$ is the number of
Monte Carlo samples. In the discrete-time case, $\hat{r}_{t}^{2}$ is pre-
computed for any $t$ and reused in subsequent computations. In the continuous-
time case, we pre-compute $\hat{r}_{t}^{2}$ for 1000 discrete time steps, and
for any $t$ used in sampling, we use the pre-computed result with the nearest
time step to $t$.

###  4.3 Comparison With Tweedie Moment Projection

Concurrent to our work, Boys et al. (2023) also leverage the idea of
optimizing posterior covariance for solving inverse problems. Specifically,
they leverage second-order Tweedie’s formula to achieve posterior covariance
optimization:
${\bm{\Sigma}}^{*}({\bm{x}}_{t})=\sigma_{t}^{2}\nabla_{{\bm{x}}_{t}}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]$,
where ${\bm{\Sigma}}_{t}^{*}({\bm{x}}_{t})$ is the optimal
${\bm{\Sigma}}_{t}({\bm{x}}_{t})$ for equation 12 and
$\nabla_{{\bm{x}}_{t}}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]$ is the Jacobian
of $\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]$. However, without strong
approximation, their approach cannot scale to high dimensional data like
images. In particular, they leverage row sum approximation for the Jacobian:
$\nabla_{{\bm{x}}_{t}}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]\approx\mathrm{diag}[\nabla_{{\bm{x}}_{t}}\mathbf{1}^{T}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]]$
for image experiments. In contrast, we achieve posterior covariance
optimization without introducing additional computational cost, as elaborated
in Sections 4.1 and 4.2.

##  5 Experiments

We implement our techniques and re-implement prior works on a newly written
codebase based on an open-source diffusion codebase
k-diffusion555https://github.com/crowsonkb/k-diffusion. This allows us to
minimize the effect of different implementations for fair comparison and
allows prior works to use a more advanced sampler that is not supported by the
original codebase. Following (Chung et al., 2023; Wang et al., 2023), we
perform experiments on the FFHQ 256$\times$256 and ImageNet 256$\times$256
datasets to compare different methods. For each dataset, we report the average
results over 100 validation images. Following (Chung et al., 2023; Zhu et al.,
2023), we use the pre-trained unconditional diffusion models from (Dhariwal &
Nichol, 2021) and (Chung et al., 2023) for ImageNet and FFHQ datasets,
respectively.

The degradation models are specified mostly following (Zhu et al., 2023): (i)
For inpainting, we randomly mask 50 percent of the total pixels. (ii) For
Gaussian debluring and motion debluring, we use the same setup of the bluring
kernels to (Chung et al., 2023). (iii) For super resolution (SR), we consider
bicubic downsampling. All measurements are corrupted by Gaussian noise with
$\sigma=0.05$. See Appendix D for more experimental details and results.

Task | Avg. MAD
---|---
Inpaint (Random) | $1.6975\times 10^{-6}$
Deblur (Gaussian) | $6.4629\times 10^{-4}$
Deblur (Motion) | $1.5461\times 10^{-3}$
Super resolution ($4\times$) | $1.1937\times 10^{-2}$
Table 2: DDNM v.s. DiffPIR in noiseless inverse problems. We report the
averaged MAD between their ${\bm{x}}_{0}^{(t)}$ averaged over all sampling
steps and test images.

Figure 1: Averaged values of ${\bm{e}}$ (black line) and ${\bm{r}}_{t}^{2}({\bm{x}}_{t})$ (blue line) on FFHQ (Left) and ImageNet (Right). Figure 2: Visualization of ${\bm{e}}$ and $\hat{{\bm{r}}}_{t}^{2}({\bm{x}}_{t})$ of an example image ${\bm{x}}_{0}$ at different $t$. Top: ${\bm{e}}$; Bottom: $\hat{{\bm{r}}}_{t}^{2}({\bm{x}}_{t})$. The results are averaged over RGB channels for better visualization. Dataset | Method | Inpaint (Random) | Deblur (Gaussian) | Deblur (Motion) | Super resolution ($4\times$)
---|---|---|---|---|---
SSIM $\uparrow$ |  LPIPS $\downarrow$ |  FID $\downarrow$ |  SSIM $\uparrow$ |  LPIPS $\downarrow$ |  FID $\downarrow$ |  SSIM $\uparrow$ |  LPIPS $\downarrow$ |  FID $\downarrow$ |  SSIM $\uparrow$ |  LPIPS $\downarrow$ |  FID $\downarrow$
FFHQ |  Convert (Ours) | 0.9279 | 0.0794 | 25.90 | 0.7905 | 0.1836 | 52.42 | 0.7584 | 0.2156 | 62.88 | 0.7878 | 0.1962 | 58.37
Analytic (Ours) | 0.9272 | 0.0845 | 28.83 | 0.7926 | 0.1850 | 53.09 | 0.7579 | 0.2183 | 64.77 | 0.7878 | 0.1968 | 59.83
TMPD | 0.8224 | 0.1924 | 70.92 | 0.7289 | 0.2523 | 76.52 | 0.7014 | 0.2718 | 83.49 | 0.7085 | 0.2701 | 79.58
DPS | 0.8891 | 0.1323 | 49.46 | 0.6284 | 0.3652 | 136.12 | 0.4904 | 0.4924 | 212.48 | 0.7719 | 0.2054 | 61.36
$\Pi\text{GDM}$ | 0.8784 | 0.1422 | 49.89 | 0.7890 | 0.1910 | 59.93 | 0.7543 | 0.2209 | 66.14 | 0.7850 | 0.2005 | 61.46
ImageNet |  Convert (Ours) | 0.8559 | 0.1329 | 29.14 | 0.6007 | 0.3327 | 95.23 | 0.5634 | 0.3656 | 109.61 | 0.5869 | 0.3477 | 96.76
Analytic (Ours) | 0.8481 | 0.1446 | 35.51 | 0.6009 | 0.3334 | 93.21 | 0.5611 | 0.3668 | 113.39 | 0.5958 | 0.3495 | 95.33
TMPD | 0.7011 | 0.2892 | 293.80 | 0.5430 | 0.4114 | 291.29 | 0.4773 | 0.4567 | 302.40 | 0.5186 | 0.4298 | 296.73
DPS | 0.8623 | 0.1490 | 36.58 | 0.4603 | 0.4630 | 173.77 | 0.3582 | 0.5554 | 282.21 | 0.5860 | 0.3231 | 92.89
$\Pi\text{GDM}$ | 0.7658 | 0.2328 | 64.96 | 0.5946 | 0.3429 | 102.89 | 0.5534 | 0.3781 | 113.89 | 0.5925 | 0.3552 | 100.36
Table 3: Quantitative results on FFHQ and ImageNet dataset for Type I
guidance. We use bold and underline for the best and second best,
respectively. Figure 3: LPIPS comparisons on FFHQ for Type II guidance. We
report LPIPS under different $\lambda$.

###  5.1 Validation of Theoretical Results

DDNM as noiseless DiffPIR. To validate Proposition 2, we re-implement DDNM
under DiffPIR codebase. DiffPIR deals with noiseless inverse problems by
setting $\sigma$ to a relatively low value for $\rho_{t}$ in equation 9
(_e.g._ , 0.001 in DiffPIR codebase). We demonstrate that directly using the
DDNM solution in equation 11 produces similar results in the noiseless case.
Table 2 reports the Mean Absolute Difference (MAD)666MAD between ${\bm{x}}$
and ${\bm{y}}$ is defined by $\lVert{\bm{x}}-{\bm{y}}\rVert_{1}/d$. between
the conditional posterior means in equations 9 and 11 by averaging over all
the sampling steps and test images. Overall, the MAD is negligible in
comparison to the data range of $[-1,1]$. Debluring and super resolution yield
larger MAD than inpainting, since several approximations are made for
computing ${\bm{A}}^{\dagger}$, while for inpainting ${\bm{A}}^{\dagger}$ is
exact (see Appendix D.3 for details).

Sanity check to equation 22. To validate the effectiveness of equation 22 in
practice, we compare the ground-truth square errors made by the denoiser
${\bm{e}}=({\bm{x}}_{0}-D_{t}({\bm{x}}_{t}))^{2}$ with the posterior variance
prediction $\hat{{\bm{r}}}_{t}^{2}({\bm{x}}_{t})$. We compare ${\bm{e}}$ with
$\hat{{\bm{r}}}_{t}^{2}({\bm{x}}_{t})$ since the posterior variance
${\bm{r}}_{t}^{*2}({\bm{x}}_{t})$ is the MMSE estimator of ${\bm{e}}$ given
${\bm{x}}_{t}$ (assuming
$D_{t}({\bm{x}}_{t})=\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]$), as suggested by
equation 16. This implies that a good $\hat{{\bm{r}}}_{t}^{2}({\bm{x}}_{t})$
should be a reliable predictor of ${\bm{e}}$. Thus, we compare ${\bm{e}}$ with
$\hat{{\bm{r}}}_{t}^{2}({\bm{x}}_{t})$ as a sanity check to gauge the
effectiveness of equation 22 in practice. In Figure 1, we plot their averaged
values over all the pixels and test images. The posterior variances prediction
obtained via equation 22 is accurate only in the regions with low noise level.
This is reasonable since the optimal reverse variance
${\bm{v}}^{*2}({\bm{x}}_{t})$ is bounded by the upper bound $\beta_{t}$ and
the lower bound $\tilde{\beta}_{t}$. The two bounds are almost equal in the
regions with high noise level (Nichol & Dhariwal, 2021). When $t$ is large,
${\bm{v}}^{*2}({\bm{x}}_{t})\\!-\\!\tilde{\beta}_{t}\\!\approx\\!\bm{0}$ and
equation 22 becomes a $0/0$ limit that possesses high numerical instability.
In Figure 2, we visualize ${\bm{e}}$ and
$\hat{{\bm{r}}}_{t}^{2}({\bm{x}}_{t})$ of an image ${\bm{x}}_{0}$ at different
$t$ and show that ${\bm{e}}$ can be well predicted by
$\hat{{\bm{r}}}_{t}^{2}({\bm{x}}_{t})$ when $t$ is small.

###  5.2 Quantitative Results

Experimental setup. To evaluate different methods, we follow Chung et al.
(2023) to use three metrics: Structure Similarity Index Measure (SSIM),
Learned Perceptual Image Patch Similarity (LPIPS, Zhang et al. (2018)) and
Frechet Inception Distance (FID, Heusel et al. (2017)). All Type I guidance
methods use the same Heun’s 2nd deterministic sampler suggested in Karras et
al. (2022) with 50 sampling steps, and all Type II guidance methods use the
same Heun’s 2nd stochastic sampler
($S_{\text{churn}}=80,S_{\text{tmin}}=0.05,S_{\text{tmax}}=1,S_{\text{noise}}=1.007$,
definition see Karras et al. (2022)) with 50 sampling steps since we found
that Type II guidance does not perform well using deterministic samplers777To
understand how to leverage pre-trained DDPM model to perform sampling under
perturbation kernels given in Section 2.2, please refer to Appendix D.2..

Dataset | Method | Inpaint (Random) | Deblur (Gaussian) | Deblur (Motion) | Super resolution ($4\times$)
---|---|---|---|---|---
SSIM $\uparrow$ |  LPIPS $\downarrow$ |  FID $\downarrow$ |  SSIM $\uparrow$ |  LPIPS $\downarrow$ |  FID $\downarrow$ |  SSIM $\uparrow$ |  LPIPS $\downarrow$ |  FID $\downarrow$ |  SSIM $\uparrow$ |  LPIPS $\downarrow$ |  FID $\downarrow$
FFHQ |  Convert (Ours) | 0.9241 | 0.0822 | 27.50 | 0.7783 | 0.1969 | 59.31 | 0.7329 | 0.2324 | 66.18 | 0.7632 | 0.2183 | 67.17
Analytic (Ours) | 0.9232 | 0.0852 | 28.63 | 0.7790 | 0.1971 | 59.80 | 0.7331 | 0.2336 | 69.82 | 0.7622 | 0.2195 | 68.84
$\Pi\text{GDM}$ | 0.7078 | 0.2605 | 77.46 | 0.7221 | 0.2421 | 71.19 | 0.6977 | 0.2607 | 75.15 | 0.7205 | 0.2442 | 72.41
ImageNet |  Convert (Ours) | 0.8492 | 0.1394 | 33.46 | 0.5770 | 0.3568 | 100.50 | 0.5341 | 0.3944 | 131.48 | 0.5613 | 0.3846 | 116.27
Analytic (Ours) | 0.8417 | 0.1484 | 37.98 | 0.5722 | 0.3583 | 103.67 | 0.5341 | 0.3947 | 126.39 | 0.5537 | 0.3867 | 116.87
$\Pi\text{GDM}$ | 0.5102 | 0.4293 | 141.80 | 0.5071 | 0.4095 | 127.41 | 0.4887 | 0.4267 | 132.38 | 0.5150 | 0.4098 | 127.40
Table 4: Results for $\Pi\text{GDM}$ with adaptive weight on FFHQ and ImageNet
dataset. We use bold and underline for the best and second best, respectively.
Figure 4: LPIPS comparisons on FFHQ for DPS with heuristic guidance strength.
We report LPIPS under different $\zeta$.

In initial experiments, using our optimal posterior covariance for all
sampling steps results in poor performance. This may due to Gaussian
variational $q_{t}({\bm{x}}_{0}|{\bm{x}}_{t})$ is an effective approximation
to posterior $p_{t}({\bm{x}}_{0}|{\bm{x}}_{t})$ only for small noise level
(Xiao et al., 2022). To address this, we heuristically use the optimal
posterior covariance only at the last several sampling steps with low noise
level but use $\Pi\text{GDM}$ and DiffPIR covariances for Type I and Type II
guidance for high noise level, respectively. We empirically find that sampling
with our optimal covariance yields highly stable results when $\sigma_{t}<0.2$
(_i.e._ , 12 out of 50 steps).

Results for pure covariance comparisons. Table 3 summarizes the results for
Type I guidance. DPS and $\Pi\text{GDM}$ refer to the posterior covariance
type discussed in Section 3.1; Convert, Analytic, and TMPD refer to the
posterior covariance type obtained using approaches presented in Sections 4.1,
4.2, and 4.3, respectively. Note that for Type I guidance, our methods achieve
the best results in almost all tasks. Although DPS outperforms us in several
cases, we observe that its performance is very unstable (see DPS performance
in debluring tasks). For Type II guidance, since we use DiffPIR covariance in
high noise level region, the performance of our method and DiffPIR baseline
are both influenced by the hyper-parameter $\lambda$. Therefore, we report the
performance under different $\lambda$ in Figure 3. Our method is more robust
than DiffPIR to the hyper-parameter $\lambda$, which demonstrates that the
last sampling steps are near optimal when using optimal posterior covariance.

Results with heuristics implementations. Prior works propose several
heuristics to enhance performance, in addition to using Gaussian for posterior
approximation. The results in Table 3 are based on re-implementations that
exclude these heuristics for pure covariance comparisons. We further
investigate the effectiveness of our techniques by applying them to the last
sampling steps of these re-implementations.

$\Pi\text{GDM}$ introduces an adaptive weight to the likelihood score
$\nabla_{{\bm{x}}_{t}}\log p_{t}({\bm{y}}|{\bm{x}}_{t})$ to adjust the
guidance strength according to the timestep. Specifically, the likelihood
score is substituted with $r_{t}^{2}\nabla_{{\bm{x}}_{t}}\log
p_{t}({\bm{y}}|{\bm{x}}_{t})$ with
$r_{t}^{2}=\sigma_{t}^{2}/(1+\sigma_{t}^{2})$. Table 4 summarizes the results
for $\Pi\text{GDM}$ with adaptive weight. Obviously, $\Pi\text{GDM}$ achieves
significant performance gain using Type I guidance with optimal variances at
the last sampling steps. As for DPS, it introduces a trick to heuristically
determine the strength of guidance as mentioned in Section 3.1. Since this
trick involves a hyper-parameter $\zeta$, we report the performance under
different $\zeta$. Figure 4 summarizes the quantitative results for DPS with
heuristic guidance strength. DPS with the last sampling steps replaced by our
techniques can achieve the best results across all hyper-parameters.

##  6 Conclusions

We show that recent diffusion-based inverse problem solvers can be interpreted
as isotropic Gaussian approximations to the intractable posterior
distributions. Inspired by this, we propose to enhance these existing methods
by optimizing the posterior covariances using two novel plug-and-play
strategies. Experimental results demonstrate that the proposed methods
significantly outperform existing zero-shot methods and enhance the robustness
to hyperparameters.

Limitations Despite the reasonable improvements our approaches offer over
prior works, it is worth noting that the Gaussian assumption with a diagonal
covariance might be an oversimplification to approximate the true pixel-
correlated image posterior. This could potentially limit the performance.
Future research could explore better principles for designing the covariance,
such as leveraging orthogonal transformation for correlation reduction (see
Appendix D.5 for discussion). Additionally, methods like efficient
approximation for Tweedie moment projection could be investigated. These
represent intriguing directions for future work.

## References

  * Ballé et al. (2020) Ballé, J., Chou, P. A., Minnen, D., Singh, S., Johnston, N., Agustsson, E., Hwang, S. J., and Toderici, G.  Nonlinear transform coding.  _IEEE Journal of Selected Topics in Signal Processing_ , 15(2):339–353, 2020.
  * Bao et al. (2022a) Bao, F., Li, C., Sun, J., Zhu, J., and Zhang, B.  Estimating the optimal covariance with imperfect mean in diffusion probabilistic models.  In _Proceedings of the 39th International Conference on Machine Learning_ , pp. 1555–1584, 2022a.
  * Bao et al. (2022b) Bao, F., Li, C., Zhu, J., and Zhang, B.  Analytic-DPM: an analytic estimate of the optimal reverse variance in diffusion probabilistic models.  In _The 10th International Conference on Learning Representations_ , 2022b.  URL https://openreview.net/forum?id=0xiJLKH-ufZ.
  * Bishop (2006) Bishop, C. M.  _Pattern Recognition and Machine Learning_.  Springer, 2006.
  * Boys et al. (2023) Boys, B., Girolami, M., Pidstrigach, J., Reich, S., Mosca, A., and Akyildiz, O. D.  Tweedie moment projected diffusions for inverse problems.  _arXiv preprint arXiv:2310.06721_ , 2023.
  * Chan et al. (2023) Chan, M. A., Young, S. I., and Metzler, C. A.  SUD2: Supervision by denoising diffusion models for image reconstruction.  In _NeurIPS 2023 Deep Inverse Workshop_ , 2023.
  * Choi et al. (2021) Choi, J., Kim, S., Jeong, Y., Gwon, Y., and Yoon, S.  Ilvr: Conditioning method for denoising diffusion probabilistic models.  _arXiv preprint arXiv:2108.02938_ , 2021.
  * Chung et al. (2023) Chung, H., Kim, J., Mccann, M. T., Klasky, M. L., and Ye, J. C.  Diffusion posterior sampling for general noisy inverse problems.  In _The 11th International Conference on Learning Representations_ , 2023.  URL https://openreview.net/forum?id=OnD9zGAGT0k.
  * Dhariwal & Nichol (2021) Dhariwal, P. and Nichol, A.  Diffusion models beat GANs on image synthesis.  In _Advances in Neural Information Processing Systems 34_ , pp. 8780–8794, 2021.
  * Feng et al. (2023) Feng, B. T., Smith, J., Rubinstein, M., Chang, H., Bouman, K. L., and Freeman, W. T.  Score-based diffusion models as principled priors for inverse imaging.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pp. 10520–10531, 2023.
  * Hestenes et al. (1952) Hestenes, M. R., Stiefel, E., et al.  Methods of conjugate gradients for solving linear systems.  _Journal of research of the National Bureau of Standards_ , 49(6):409–436, 1952.
  * Heusel et al. (2017) Heusel, M., Ramsauer, H., Unterthiner, T., Nessler, B., and Hochreiter, S.  Gans trained by a two time-scale update rule converge to a local nash equilibrium.  In _Advances in Neural Information Processing Systems 30_ , 2017.
  * Ho et al. (2020) Ho, J., Jain, A., and Abbeel, P.  Denoising diffusion probabilistic models.  In _Advances in Neural Information Processing Systems 33_ , pp. 6840–6851, 2020.
  * Karras et al. (2022) Karras, T., Aittala, M., Aila, T., and Laine, S.  Elucidating the design space of diffusion-based generative models.  In _Advances in Neural Information Processing Systems 35_ , pp. 26565–26577, 2022.
  * Kingma et al. (2021) Kingma, D., Salimans, T., Poole, B., and Ho, J.  Variational diffusion models.  In _Advances in Neural Information Processing Systems 34_ , pp. 21696–21707, 2021.
  * Lugmayr et al. (2022) Lugmayr, A., Danelljan, M., Romero, A., Yu, F., Timofte, R., and Van Gool, L.  Repaint: Inpainting using denoising diffusion probabilistic models.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pp. 11461–11471, 2022.
  * Luo et al. (2023) Luo, Z., Gustafsson, F. K., Zhao, Z., Sjölund, J., and Schön, T. B.  Refusion: Enabling large-size realistic image restoration with latent-space diffusion models.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pp. 1680–1691, 2023.
  * Mardani et al. (2023) Mardani, M., Song, J., Kautz, J., and Vahdat, A.  A variational perspective on solving inverse problems with diffusion models.  _arXiv preprint arXiv:2305.04391_ , 2023.
  * Nichol & Dhariwal (2021) Nichol, A. Q. and Dhariwal, P.  Improved denoising diffusion probabilistic models.  In _Proceedings of the 38th International Conference on Machine Learning_ , pp. 8162–8171, 2021.
  * Pokle et al. (2023) Pokle, A., Muckley, M. J., Chen, R. T., and Karrer, B.  Training-free linear image inversion via flows.  _arXiv preprint arXiv:2310.04432_ , 2023.
  * Ravula et al. (2023) Ravula, S., Levac, B., Jalal, A., Tamir, J. I., and Dimakis, A. G.  Optimizing sampling patterns for compressed sensing MRI with diffusion generative models.  _arXiv preprint arXiv:2306.03284_ , 2023.
  * Rezende & Viola (2018) Rezende, D. J. and Viola, F.  Taming VAEs.  _arXiv preprint arXiv:1810.00597_ , 2018.
  * Saharia et al. (2022) Saharia, C., Ho, J., Chan, W., Salimans, T., Fleet, D. J., and Norouzi, M.  Image super-resolution via iterative refinement.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 45(4):4713–4726, 2022.
  * Sohl-Dickstein et al. (2015) Sohl-Dickstein, J., Weiss, E., Maheswaranathan, N., and Ganguli, S.  Deep unsupervised learning using nonequilibrium thermodynamics.  In _Proceeding of the 35th International Conference on Machine Learning_ , pp. 2256–2265, 2015.
  * Song et al. (2021a) Song, J., Meng, C., and Ermon, S.  Denoising diffusion implicit models.  In _The 9th International Conference on Learning Representations_ , 2021a.  URL https://openreview.net/forum?id=St1giarCHLP.
  * Song et al. (2023) Song, J., Vahdat, A., Mardani, M., and Kautz, J.  Pseudoinverse-guided diffusion models for inverse problems.  In _The 11th International Conference on Learning Representations_ , 2023.  URL https://openreview.net/forum?id=9_gsMA8MRKQ.
  * Song & Ermon (2019) Song, Y. and Ermon, S.  Generative modeling by estimating gradients of the data distribution.  In _Advances in Neural Information Processing Systems 32_ , 2019.
  * Song et al. (2021b) Song, Y., Sohl-Dickstein, J., Kingma, D. P., Kumar, A., Ermon, S., and Poole, B.  Score-based generative modeling through stochastic differential equations.  In _The 9th International Conference on Learning Representations_ , 2021b.  URL https://openreview.net/forum?id=PxTIG12RRHS.
  * Song et al. (2022) Song, Y., Shen, L., Xing, L., and Ermon, S.  Solving inverse problems in medical imaging with score-based generative models.  In _The 10th International Conference on Learning Representations_ , 2022.  URL https://openreview.net/forum?id=vaRCHVj0uGI.
  * Wang et al. (2023) Wang, Y., Yu, J., and Zhang, J.  Zero-shot image restoration using denoising diffusion null-space model.  In _The 11th International Conference on Learning Representations_ , 2023.  URL https://openreview.net/forum?id=mRieQgMtNTQ.
  * Whang et al. (2022) Whang, J., Delbracio, M., Talebi, H., Saharia, C., Dimakis, A. G., and Milanfar, P.  Deblurring via stochastic refinement.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pp. 16293–16303, 2022.
  * Xiao et al. (2022) Xiao, Z., Kreis, K., and Vahdat, A.  Tackling the generative learning trilemma with denoising diffusion GANs.  In _The 10th International Conference on Learning Representations_ , 2022.  URL https://openreview.net/forum?id=JprM0p-q0Co.
  * Zhang et al. (2020) Zhang, K., Gool, L. V., and Timofte, R.  Deep unfolding network for image super-resolution.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pp. 3217–3226, 2020.
  * Zhang et al. (2018) Zhang, R., Isola, P., Efros, A. A., Shechtman, E., and Wang, O.  The unreasonable effectiveness of deep features as a perceptual metric.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , pp. 586–595, 2018.
  * Zhu et al. (2023) Zhu, Y., Zhang, K., Liang, J., Cao, J., Wen, B., Timofte, R., and Van Gool, L.  Denoising diffusion models for plug-and-play image restoration.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pp. 1219–1229, 2023.

##  Appendix A Proofs to Theoretical Results

######  Lemma 5 (Tweedie’s formula).

If the joint distribution between ${\bm{x}}_{0},{\bm{x}}_{t}$ is given by
$p_{t}({\bm{x}}_{0},{\bm{x}}_{t})=p({\bm{x}}_{0})p_{t}({\bm{x}}_{t}|{\bm{x}}_{0})$
with
$p_{t}({\bm{x}}_{t}|{\bm{x}}_{0})=\mathcal{N}({\bm{x}}_{t}|s_{t}{\bm{x}}_{0},s_{t}^{2}\sigma_{t}^{2}{\bm{I}})$
, then $\nabla_{{\bm{x}}_{t}}\log
p_{t}({\bm{x}}_{t})=\frac{1}{s_{t}^{2}\sigma_{t}^{2}}(s_{t}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]-{\bm{x}}_{t})$.

###### Proof.

$\displaystyle\nabla_{{\bm{x}}_{t}}\log p_{t}({\bm{x}}_{t})$
$\displaystyle=\frac{\nabla_{{\bm{x}}_{t}}p_{t}({\bm{x}}_{t})}{p_{t}({\bm{x}}_{t})}$
(25) $\displaystyle=\frac{1}{p_{t}({\bm{x}}_{t})}\nabla_{{\bm{x}}_{t}}\int
p({\bm{x}}_{0})p_{t}({\bm{x}}_{t}|{\bm{x}}_{0})\mathrm{d}{\bm{x}}_{0}$ (26)
$\displaystyle=\frac{1}{p_{t}({\bm{x}}_{t})}\int
p({\bm{x}}_{0})\nabla_{{\bm{x}}_{t}}p_{t}({\bm{x}}_{t}|{\bm{x}}_{0})\mathrm{d}{\bm{x}}_{0}$
(27) $\displaystyle=\frac{1}{p_{t}({\bm{x}}_{t})}\int
p({\bm{x}}_{0})p_{t}({\bm{x}}_{t}|{\bm{x}}_{0})\nabla_{{\bm{x}}_{t}}\log
p_{t}({\bm{x}}_{t}|{\bm{x}}_{0})\mathrm{d}{\bm{x}}_{0}$ (28)
$\displaystyle=\int p_{t}({\bm{x}}_{0}|{\bm{x}}_{t})\nabla_{{\bm{x}}_{t}}\log
p_{t}({\bm{x}}_{t}|{\bm{x}}_{0})\mathrm{d}{\bm{x}}_{0}$ (29)
$\displaystyle=\mathbb{E}_{p_{t}({\bm{x}}_{0}|{\bm{x}}_{t})}[\nabla_{{\bm{x}}_{t}}\log
p_{t}({\bm{x}}_{t}|{\bm{x}}_{0})]$ (30)

For Gaussian perturbation kernel
$p_{t}({\bm{x}}_{t}|{\bm{x}}_{0})=\mathcal{N}({\bm{x}}_{t}|s_{t}{\bm{x}}_{0},s_{t}^{2}\sigma_{t}^{2}{\bm{I}})$,
we have $\nabla_{{\bm{x}}_{t}}\log
p_{t}({\bm{x}}_{t}|{\bm{x}}_{0})=\frac{1}{s_{t}^{2}\sigma_{t}^{2}}(s_{t}{\bm{x}}_{0}-{\bm{x}}_{t})$.
Plug it into equation 30, we conclude the proof. ∎

######  Lemma 6 (Conditional Tweedie’s formula).

If the joint distribution between ${\bm{x}}_{0},{\bm{y}},{\bm{x}}_{t}$ is
given by
$p_{t}({\bm{x}}_{0},{\bm{y}},{\bm{x}}_{t})=p({\bm{x}}_{0})p({\bm{y}}|{\bm{x}}_{0})p_{t}({\bm{x}}_{t}|{\bm{x}}_{0})$
with
$p_{t}({\bm{x}}_{t}|{\bm{x}}_{0})=\mathcal{N}({\bm{x}}_{t}|s_{t}{\bm{x}}_{0},s_{t}^{2}\sigma_{t}^{2}{\bm{I}})$
, then $\nabla_{{\bm{x}}_{t}}\log
p_{t}({\bm{x}}_{t}|{\bm{y}})=\frac{1}{s_{t}^{2}\sigma_{t}^{2}}(s_{t}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}}]-{\bm{x}}_{t})$

###### Proof.

$\displaystyle\nabla_{{\bm{x}}_{t}}\log p_{t}({\bm{x}}_{t}|{\bm{y}})$
$\displaystyle=\frac{\nabla_{{\bm{x}}_{t}}p_{t}({\bm{x}}_{t}|{\bm{y}})}{p_{t}({\bm{x}}_{t}|{\bm{y}})}$
(31)
$\displaystyle=\frac{1}{p_{t}({\bm{x}}_{t}|{\bm{y}})}\nabla_{{\bm{x}}_{t}}\int
p_{t}({\bm{x}}_{t}|{\bm{x}}_{0},{\bm{y}})p({\bm{x}}_{0}|{\bm{y}})\mathrm{d}{\bm{x}}_{0}$
(32)
$\displaystyle=\frac{1}{p_{t}({\bm{x}}_{t}|{\bm{y}})}\nabla_{{\bm{x}}_{t}}\int
p_{t}({\bm{x}}_{t}|{\bm{x}}_{0})p({\bm{x}}_{0}|{\bm{y}})\mathrm{d}{\bm{x}}_{0}$
(33) $\displaystyle=\frac{1}{p_{t}({\bm{x}}_{t}|{\bm{y}})}\int
p({\bm{x}}_{0}|{\bm{y}})\nabla_{{\bm{x}}_{t}}p_{t}({\bm{x}}_{t}|{\bm{x}}_{0})\mathrm{d}{\bm{x}}_{0}$
(34) $\displaystyle=\frac{1}{p_{t}({\bm{x}}_{t}|{\bm{y}})}\int
p({\bm{x}}_{0}|{\bm{y}})p_{t}({\bm{x}}_{t}|{\bm{x}}_{0},{\bm{y}})\nabla_{{\bm{x}}_{t}}\log
p_{t}({\bm{x}}_{t}|{\bm{x}}_{0})\mathrm{d}{\bm{x}}_{0}$ (35)
$\displaystyle=\int
p_{t}({\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}})\nabla_{{\bm{x}}_{t}}\log
p_{t}({\bm{x}}_{t}|{\bm{x}}_{0})\mathrm{d}{\bm{x}}_{0}$ (36)
$\displaystyle=\mathbb{E}_{p_{t}({\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}})}[\nabla_{{\bm{x}}_{t}}\log
p_{t}({\bm{x}}_{t}|{\bm{x}}_{0})]$ (37)

where equation 33 and equation 35 are due to the conditional independence
between ${\bm{x}}_{t}$ and ${\bm{y}}$ given ${\bm{x}}_{0}$, such that
$p_{t}({\bm{x}}_{t}|{\bm{x}}_{0},{\bm{y}})=p_{t}({\bm{x}}_{t}|{\bm{x}}_{0})$.
For Gaussian perturbation kernel
$p_{t}({\bm{x}}_{t}|{\bm{x}}_{0})=\mathcal{N}({\bm{x}}_{t}|s_{t}{\bm{x}}_{0},s_{t}^{2}\sigma_{t}^{2}{\bm{I}})$,
we have $\nabla_{{\bm{x}}_{t}}\log
p_{t}({\bm{x}}_{t}|{\bm{x}}_{0})=\frac{1}{s_{t}^{2}\sigma_{t}^{2}}(s_{t}{\bm{x}}_{0}-{\bm{x}}_{t})$.
Plug it into equation 37, we conclude the proof. ∎

###  A.1 Derivation of the Marginal Preserving Property of Diffusion ODEs

For the sake of completeness, here we prove that the ODEs given in equation 2
and equation 4 respectively maintain the exact same marginals to
$p_{t}({\bm{x}}_{t})$ and $p_{t}({\bm{x}}_{t}|{\bm{y}})$.

###### Proof.

By borrowing the results from (Equation 4, (Karras et al., 2022)) and setting
$s_{t}=1,\sigma_{t}=t$, ${\bm{x}}_{t}$ determined by the following ODE
preserves the marginal $p_{t}({\bm{x}}_{t})$ for all $t\in[0,T]$:

$\mathrm{d}{\bm{x}}_{t}=-t\nabla_{{\bm{x}}_{t}}\log
p_{t}({\bm{x}}_{t})\mathrm{d}t,\ \ {\bm{x}}_{T}\sim p_{T}({\bm{x}}_{T})$ (38)

Using the posterior mean $\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]$ to represent
the score $\nabla_{{\bm{x}}_{t}}\log p_{t}({\bm{x}}_{t})$ using Lemma 5, we
recover equation 2:

$\mathrm{d}{\bm{x}}_{t}=\frac{{\bm{x}}_{t}-\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]}{t}\mathrm{d}t,\
\ {\bm{x}}_{T}\sim p_{T}({\bm{x}}_{T})$ (39)

Likewise, the following ODE preserves the marginal
$p_{t}({\bm{x}}_{t}|{\bm{y}})$ for all $t\in[0,T]$:

$\mathrm{d}{\bm{x}}_{t}=-t\nabla_{{\bm{x}}_{t}}\log
p_{t}({\bm{x}}_{t}|{\bm{y}})\mathrm{d}t,\ \ {\bm{x}}_{T}\sim
p_{T}({\bm{x}}_{T}|{\bm{y}})$ (40)

By Lemma 6, we recover equation 4:

$\mathrm{d}{\bm{x}}_{t}=\frac{{\bm{x}}_{t}-\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}}]}{t}\mathrm{d}t,\
\ {\bm{x}}_{T}\sim p_{T}({\bm{x}}_{T}|{\bm{y}})$ (41)

∎

###  A.2 Proof of Proposition 1

To relate $\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}}]$ to
$\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]$, we note that

$\nabla_{{\bm{x}}_{t}}\log
p_{t}({\bm{x}}_{t}|{\bm{y}})=\nabla_{{\bm{x}}_{t}}\log
p_{t}({\bm{x}}_{t})+\nabla_{{\bm{x}}_{t}}\log p_{t}({\bm{y}}|{\bm{x}}_{t})$
(42)

Using Lemma 5 and Lemma 6, we have

$\frac{1}{s_{t}^{2}\sigma_{t}^{2}}(s_{t}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}}]-{\bm{x}}_{t})=\frac{1}{s_{t}^{2}\sigma_{t}^{2}}(s_{t}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]-{\bm{x}}_{t})+\nabla_{{\bm{x}}_{t}}\log
p_{t}({\bm{y}}|{\bm{x}}_{t})$ (43)

and consequently,

$\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}}]=\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]+s_{t}\sigma_{t}^{2}\nabla_{{\bm{x}}_{t}}\log
p_{t}({\bm{y}}|{\bm{x}}_{t})$ (44)

###  A.3 Proof of Proposition 2

When $\sigma\rightarrow 0$, finding the minimizer of equation 3.2 is
equivalent to solving the following hard-constraint optimization problem:

$\min_{\hat{{\bm{x}}}_{0}}\lVert\hat{{\bm{x}}}_{0}-\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]\rVert^{2}_{2}\
\ \ \ \mathrm{s.t.}\ \ \ \ {\bm{y}}={\bm{A}}\hat{{\bm{x}}}_{0}$ (45)

We define the Lagrangian
$\mathcal{L}(\hat{{\bm{x}}}_{0},\lambda)=\frac{1}{2}\lVert\hat{{\bm{x}}}_{0}-\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]\rVert^{2}_{2}+\lambda^{T}({\bm{y}}-{\bm{A}}\hat{{\bm{x}}}_{0})$,
where $\lambda$ is the Lagrangian multiplier. By the optimality condition, we
have

$\displaystyle\nabla_{\hat{{\bm{x}}}_{0}}\mathcal{L}=\hat{{\bm{x}}}_{0}-\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]-{\bm{A}}^{T}\lambda=\bm{0}$
(46)
$\displaystyle\nabla_{\lambda}\mathcal{L}={\bm{y}}-{\bm{A}}\hat{{\bm{x}}}_{0}=\bm{0}$
(47)

Multiplying ${\bm{A}}$ to equation 46 and combining the condition of equation
47 gives:

$\displaystyle{\bm{A}}(\hat{{\bm{x}}}_{0}-\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]-{\bm{A}}^{T}\lambda)=\bm{0}$
(48)
$\displaystyle\Rightarrow{\bm{y}}-{\bm{A}}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]-{\bm{A}}{\bm{A}}^{T}\lambda=\bm{0}$
(49)
$\displaystyle\Rightarrow{\bm{A}}{\bm{A}}^{T}\lambda={\bm{y}}-{\bm{A}}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]$
(50)

Multiplying ${\bm{A}}^{\dagger}$ to equation 50 and leveraging the property
${\bm{A}}^{\dagger}{\bm{A}}{\bm{A}}^{T}={\bm{A}}^{T}$, we have

${\bm{A}}^{T}\lambda={\bm{A}}^{\dagger}{\bm{y}}-{\bm{A}}^{\dagger}{\bm{A}}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]$
(51)

and consequently,

$\displaystyle\hat{{\bm{x}}}_{0}$
$\displaystyle=\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]+{\bm{A}}^{\dagger}{\bm{y}}-{\bm{A}}^{\dagger}{\bm{A}}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]$
(52)
$\displaystyle={\bm{A}}^{\dagger}{\bm{y}}+({\bm{I}}-{\bm{A}}^{\dagger}{\bm{A}})\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]$
(53)

###  A.4 Proof of Theorem 4

We follow (Bao et al., 2022b) to derive the relationship between the optimal
reverse variances ${\bm{v}}_{t}^{2}({\bm{x}}_{t})$ and optimal posterior
variances ${\bm{r}}_{t}^{2}({\bm{x}}_{t})$ based on more general non-Markov
forward process introduced by (Song et al., 2021a). To find the optimal
solution to equation 19, we present a much simpler proof than (Bao et al.,
2022b) using functional derivatives motivated by (Rezende & Viola, 2018).
Given a noise schedule $\\{\beta_{t}\\}_{t=1}^{T}$ and
$\alpha_{t}=1-\beta_{t}$, the forward process $q({\bm{x}}_{1:T}|{\bm{x}}_{0})$
is defined as

$\displaystyle
q({\bm{x}}_{1:T}|{\bm{x}}_{0})=q({\bm{x}}_{T}|{\bm{x}}_{0})\prod_{t=2}^{T}q({\bm{x}}_{t-1}|{\bm{x}}_{t},{\bm{x}}_{0})$
(54) $\displaystyle
q({\bm{x}}_{t-1}|{\bm{x}}_{t},{\bm{x}}_{0})=\mathcal{N}({\bm{x}}_{t-1}|\tilde{{\bm{\mu}}}_{t}({\bm{x}}_{t},{\bm{x}}_{0}),\lambda_{t}^{2}{\bm{I}})$
(55)
$\displaystyle\tilde{{\bm{\mu}}}_{t}({\bm{x}}_{t},{\bm{x}}_{0})=\sqrt{\bar{\alpha}_{t-1}}{\bm{x}}_{0}+\sqrt{\bar{\beta}_{t-1}-\lambda_{t}^{2}}\cdot\frac{{\bm{x}}_{t}-\sqrt{\bar{\alpha}_{t}}{\bm{x}}_{0}}{\sqrt{\bar{\beta}_{t}}}$
(56)

where $\bar{\alpha}_{t}=\prod_{i=1}^{t}\alpha_{i}$ and
$\bar{\beta}_{t}=1-\bar{\alpha}_{t}$. (Song et al., 2021a) show that for
arbitrary choice of $\lambda_{t}$, the marginal distributions
$q({\bm{x}}_{t}|{\bm{x}}_{0})$ maintain
$q({\bm{x}}_{t}|{\bm{x}}_{0})=\mathcal{N}({\bm{x}}_{t}|\sqrt{\bar{\alpha}_{t}}{\bm{x}}_{0},\bar{\beta}_{t}{\bm{I}})$.
DDPM forward process is a special case when
$\lambda_{t}^{2}=\tilde{\beta}_{t}$ with
$\tilde{\beta}_{t}=\frac{\bar{\beta}_{t-1}}{\bar{\beta}_{t}}\beta_{t}$, which
is used in (Nichol & Dhariwal, 2021) for pre-training DDPM model.

To fit the data distribution $q({\bm{x}}_{0})$, we define the reverse process,
given by

$p({\bm{x}}_{0:T})=p({\bm{x}}_{T})\prod_{t=1}^{T}p({\bm{x}}_{t-1}|{\bm{x}}_{t}),\
\
p({\bm{x}}_{t-1}|{\bm{x}}_{t})=\mathcal{N}({\bm{x}}_{t-1}|{\bm{m}}_{t}({\bm{x}}_{t}),\mathrm{diag}[{\bm{v}}_{t}^{2}({\bm{x}}_{t})])$
(57)

To train $p({\bm{x}}_{0})$ we minimize the KL divergence between the forward
and the reverse process:

$\min_{p}D_{KL}(q({\bm{x}}_{0:T})||p({\bm{x}}_{0:T})),\ \
q({\bm{x}}_{0:T})=q({\bm{x}}_{0})q({\bm{x}}_{1:T}|{\bm{x}}_{0})$ (58)

which is equivalent to minimizing the variational bound
$\mathbb{E}_{q}[L_{\text{vb}}]$ on negative log-likelihood of data
distribution $q({\bm{x}}_{0})$ with $L_{\text{vb}}$ given as follows:

$\displaystyle L_{\text{vb}}=L_{0}+L_{1}+...+L_{T}$ (59) $\displaystyle
L_{0}=-\log p({\bm{x}}_{0}|{\bm{x}}_{1})$ (60) $\displaystyle
L_{t-1}=D_{KL}(q({\bm{x}}_{t-1}|{\bm{x}}_{0},{\bm{x}}_{t})||p({\bm{x}}_{t-1}|{\bm{x}}_{t}))$
(61) $\displaystyle
L_{T}=D_{KL}(q({\bm{x}}_{T}|{\bm{x}}_{0})||p({\bm{x}}_{T}))$ (62)

For $t\in[2,T]$, $L_{t-1}$ are KL divergences between two Gaussians, which
possess analytical forms:

$\displaystyle L_{t-1}$
$\displaystyle\equiv\log\frac{|\mathrm{diag}[{\bm{v}}_{t}^{2}({\bm{x}}_{t})]|}{|\lambda_{t}^{2}{\bm{I}}|}+\lVert\tilde{{\bm{\mu}}}_{t}({\bm{x}}_{t},{\bm{x}}_{0})-{\bm{m}}({\bm{x}}_{t})\rVert^{2}_{\mathrm{diag}[{\bm{v}}_{t}^{2}({\bm{x}}_{t})]^{-1}}+\mathrm{tr}[\lambda_{t}^{2}\mathrm{diag}[{\bm{v}}_{t}^{2}({\bm{x}}_{t})]^{-1}]$
(63)
$\displaystyle\equiv\sum_{i=1}^{d}\log{\bm{v}}_{t}^{2}({\bm{x}}_{t})_{i}+\frac{(\tilde{{\bm{\mu}}}_{t}({\bm{x}}_{t},{\bm{x}}_{0})_{i}-{\bm{m}}_{t}({\bm{x}}_{t})_{i})^{2}}{{\bm{v}}_{t}^{2}({\bm{x}}_{t})_{i}}+\frac{\lambda_{t}^{2}}{{\bm{v}}_{t}^{2}({\bm{x}}_{t})_{i}}$
(64)

where “$\equiv$” denotes “equals up to a constant and a scaling factor” and
$i$ indexes the elements of an vector.

Note that minimizing $\mathbb{E}_{q}[L_{\text{vb}}]$ can be decomposed into
$T$ independent optimization sub-problems:

$\min_{{\bm{m}}_{t},{\bm{v}}_{t}}\mathbb{E}_{q({\bm{x}}_{0},{\bm{x}}_{t})}[L_{t-1}],\
\ t\in[1,T]$ (65)

The optimal ${\bm{m}}_{t}$ and ${\bm{v}}_{t}$ can be found by taking the
functional derivatives of $\mathbb{E}_{q({\bm{x}}_{0},{\bm{x}}_{t})}[L_{t-1}]$
w.r.t ${\bm{m}}_{t}$ and ${\bm{v}}_{t}^{2}$ then set to zero:

$\displaystyle\frac{\delta\mathbb{E}_{q({\bm{x}}_{0},{\bm{x}}_{t})}[L_{t-1}]}{\delta{\bm{m}}_{t}({\bm{x}}_{t})_{i}}\equiv\mathbb{E}_{q({\bm{x}}_{0})}[q({\bm{x}}_{t}|{\bm{x}}_{0})\frac{{\bm{m}}_{t}({\bm{x}}_{t})_{i}-\tilde{{\bm{\mu}}}_{t}({\bm{x}}_{t},{\bm{x}}_{0})_{i}}{{\bm{v}}_{t}^{2}({\bm{x}}_{t})_{i}}]=0$
(66)
$\displaystyle\frac{\delta\mathbb{E}_{q({\bm{x}}_{0},{\bm{x}}_{t})}[L_{t-1}]}{\delta{\bm{v}}_{t}^{2}({\bm{x}}_{t})_{i}}\equiv\mathbb{E}_{q({\bm{x}}_{0})}[q({\bm{x}}_{t}|{\bm{x}}_{0})(\frac{1}{{\bm{v}}_{t}^{2}({\bm{x}}_{t})_{i}}-\frac{(\tilde{{\bm{\mu}}}_{t}({\bm{x}}_{t},{\bm{x}}_{0})_{i}-{\bm{m}}({\bm{x}}_{t})_{i})^{2}}{({\bm{v}}_{t}^{2}({\bm{x}}_{t})_{i})^{2}}-\frac{\lambda_{t}^{2}}{({\bm{v}}_{t}^{2}({\bm{x}}_{t})_{i})^{2}})]=0$
(67)

We can solve for optimal ${\bm{m}}_{t}({\bm{x}}_{t})_{i}$ by rearranging
equation 66:

$\displaystyle\mathbb{E}_{q({\bm{x}}_{0})}[q({\bm{x}}_{t}|{\bm{x}}_{0})]{\bm{m}}_{t}({\bm{x}}_{t})_{i}=\mathbb{E}_{q({\bm{x}}_{0})}[q({\bm{x}}_{t}|{\bm{x}}_{0})\tilde{{\bm{\mu}}}_{t}({\bm{x}}_{t},{\bm{x}}_{0})_{i}]$
(68) $\displaystyle\Rightarrow
q({\bm{x}}_{t}){\bm{m}}_{t}({\bm{x}}_{t})_{i}=\int
q({\bm{x}}_{0})q({\bm{x}}_{t}|{\bm{x}}_{0})\tilde{{\bm{\mu}}}_{t}({\bm{x}}_{t},{\bm{x}}_{0})_{i}\mathrm{d}{\bm{x}}_{0}$
(69) $\displaystyle\Rightarrow{\bm{m}}_{t}({\bm{x}}_{t})_{i}=\int
q({\bm{x}}_{0}|{\bm{x}}_{t})\tilde{{\bm{\mu}}}_{t}({\bm{x}}_{t},{\bm{x}}_{0})_{i}\mathrm{d}{\bm{x}}_{0}$
(70)
$\displaystyle\Rightarrow{\bm{m}}_{t}({\bm{x}}_{t})_{i}=\tilde{{\bm{\mu}}}_{t}({\bm{x}}_{t},\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}])_{i}$
(71)

where equation 71 is due to the linearity of $\tilde{{\bm{\mu}}}_{t}$ w.r.t
${\bm{x}}_{0}$.

Likewise, rearranging equation 67 gives

${\bm{v}}_{t}^{2}({\bm{x}}_{t})_{i}=\lambda_{t}^{2}+\mathbb{E}_{q({\bm{x}}_{0}|{\bm{x}}_{t})}[(\tilde{{\bm{\mu}}}_{t}({\bm{x}}_{t},{\bm{x}}_{0})_{i}-{\bm{m}}({\bm{x}}_{t})_{i})^{2}]$
(72)

By plugging the optimal ${\bm{m}}_{t}({\bm{x}}_{t})_{i}$ determined by
equation 71 into equation 72 and dropping the element index $i$, we conclude
the proof:

$\displaystyle{\bm{v}}_{t}^{2}({\bm{x}}_{t})$
$\displaystyle=\lambda_{t}^{2}+\mathbb{E}_{q({\bm{x}}_{0}|{\bm{x}}_{t})}[(\tilde{{\bm{\mu}}}_{t}({\bm{x}}_{t},{\bm{x}}_{0})-\tilde{{\bm{\mu}}}_{t}({\bm{x}}_{t},\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]))^{2}]$
(73)
$\displaystyle=\lambda_{t}^{2}+\mathbb{E}_{q({\bm{x}}_{0}|{\bm{x}}_{t})}[(\tilde{{\bm{\mu}}}_{t}(\bm{0},{\bm{x}}_{0}-\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]))^{2}]$
(74)
$\displaystyle=\lambda_{t}^{2}+(\sqrt{\bar{\alpha}_{t-1}}-\sqrt{\bar{\beta}_{t-1}-\lambda_{t}^{2}}\sqrt{\frac{\bar{\alpha}_{t}}{\bar{\beta}_{t}}})^{2}\cdot\mathbb{E}_{q({\bm{x}}_{0}|{\bm{x}}_{t})}[({\bm{x}}_{0}-\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}])^{2}]$
(75)
$\displaystyle=\lambda_{t}^{2}+(\sqrt{\bar{\alpha}_{t-1}}-\sqrt{\bar{\beta}_{t-1}-\lambda_{t}^{2}}\sqrt{\frac{\bar{\alpha}_{t}}{\bar{\beta}_{t}}})^{2}\cdot{\bm{r}}_{t}^{2}({\bm{x}}_{t})$
(76)

We are often given a pre-trained DDPM model with learned reverse variances.
The following Corollary of Theorem 4 gives a simplified relationship between
optimal posterior variances and optimal reverse variances under DDPM case:

######  Corollary 7.

For the DDPM forward process $\lambda_{t}^{2}=\tilde{\beta}_{t}$, the optimal
posterior variances ${\bm{r}}_{t}^{*2}({\bm{x}}_{t})$ and optimal reverse
variances ${\bm{v}}_{t}^{*2}({\bm{x}}_{t})$ are related by

${\bm{v}}_{t}^{*2}({\bm{x}}_{t})=\tilde{\beta}_{t}+(\frac{\sqrt{\bar{\alpha}_{t-1}}\beta_{t}}{1-\bar{\alpha}_{t}})^{2}\cdot{\bm{r}}_{t}^{*2}({\bm{x}}_{t})$
(77)

######  Remark 8.

From Theorem 4, we also know that to perform optimal ancestral sampling, we
only need to provide the MMSE estimator
$\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]$ to compute the reverse mean
${\bm{m}}^{*}_{t}({\bm{x}}_{t})$ and the posterior variances
${\bm{r}}^{*2}_{t}({\bm{x}}_{t})$ to compute the reverse variances
${\bm{v}}^{*2}_{t}({\bm{x}}_{t})$, which can be both obtained from the
proposed pre-training (equation 12). This may provide an alternative way for
pre-trainning DDPM model that differs from (Nichol & Dhariwal, 2021).

###  A.5 Proof of Proposition 3

Deriving the optimal solution to equation 12 under the diagonal posterior
covariance case, i.e.,
${\bm{\Sigma}}_{t}({\bm{x}}_{t})=\mathrm{diag}[{\bm{r}}_{t}^{2}({\bm{x}}_{t})]$,
is similar to Appendix A.4. Note that we seek for point-wise maximizer of
equation 12, i.e., find the optimum of

$\mathbb{E}_{p_{t}({\bm{x}}_{0},{\bm{x}}_{t})}[\log
q_{t}({\bm{x}}_{0}|{\bm{x}}_{t})]\equiv\mathbb{E}_{p_{t}({\bm{x}}_{0},{\bm{x}}_{t})}[\sum_{i=1}^{d}\frac{1}{{\bm{r}}_{t}^{2}({\bm{x}}_{t})_{i}}({\bm{x}}_{0i}-{\bm{\mu}}_{t}({\bm{x}}_{t})_{i})^{2}+\log{\bm{r}}_{t}^{2}({\bm{x}}_{t})_{i}]$
(78)

For any $t$, taking the functional derivatives of
$\mathbb{E}_{p_{t}({\bm{x}}_{0},{\bm{x}}_{t})}[\log
q_{t}({\bm{x}}_{0}|{\bm{x}}_{t})]$ w.r.t ${\bm{\mu}}_{t}({\bm{x}}_{t})_{i}$
and ${\bm{r}}_{t}^{2}({\bm{x}}_{t})_{i}$ and then set to zero, we obtain the
optimality conditions:

$\displaystyle\mathbb{E}_{p({\bm{x}}_{0})}[p({\bm{x}}_{t}|{\bm{x}}_{0})({\bm{\mu}}_{t}({\bm{x}}_{t})_{i}-{\bm{x}}_{0i})]=0$
(79)
$\displaystyle\mathbb{E}_{p({\bm{x}}_{0})}[p({\bm{x}}_{t}|{\bm{x}}_{0})(-\frac{1}{({\bm{r}}_{t}^{2}({\bm{x}}_{t})_{i})^{2}}({\bm{x}}_{0i}-{\bm{\mu}}_{t}({\bm{x}}_{t})_{i})^{2}+\frac{1}{{\bm{r}}_{t}^{2}({\bm{x}}_{t})_{i}}]=0$
(80)

Combining the optimality conditions given by equation 79 and equation 80, we
conclude the proof.

##  Appendix B Closed-form Solutions for Implementing Guidances

In this section, we provide important closed-form results for implementing
efficient Type I and Type II guidance under isotropic posterior covariance
case: ${\bm{\Sigma}}_{t}({\bm{x}}_{t})=r_{t}^{2}({\bm{x}}_{t}){\bm{I}}$.
Before we delve into the closed-form results, we first give some important
notations. We define the downsampling operator given sampling position
${\bm{m}}\in\\{0,1\\}^{d\times 1}$ as ${\bm{D}}_{{\bm{m}}}\in\\{0,1\\}^{\lVert
m\rVert_{0}\times d}$, which selects rows of a given matrix that corresponds
to one in ${\bm{m}}$ and when performing left multiplication. We use
${\bm{D}}_{\downarrow s}$ to denote the standard $s$-folds downsampling
operator, which is equivalent to ${\bm{D}}_{{\bm{m}}}$ when ones in ${\bm{m}}$
are spaced evenly. For image signal, it selects the upper-left pixel for each
distinct $s\times s$ patch (Zhang et al., 2020) when performing left
multiplication to the vectorized image. We use ${\bm{D}}_{\Downarrow s}$ to
denote the distinct block downsampler, i.e., averaging $s$ length $d/s$
distinct blocks of a vector. For image signal, it averaging distinct
$d/s\times d/s$ blocks (Zhang et al., 2020) when performing left
multiplication to the vectorized image. We denote the Fourier transform matrix
for $d$-dimensional signal as ${\bm{F}}$, Fourier transform matrix for
$d/s$-dimensional signal as ${\bm{F}}_{\downarrow s}$, the Fourier transform
of a vector ${\bm{v}}$ as $\hat{{\bm{v}}}$, and the complex conjugate of a
complex vector ${\bm{v}}$ as $\bar{{\bm{v}}}$. We use $\odot$ to denote
element-wise multiplication, and the divisions used below are also element-
wise.

######  Lemma 9.

Performing $s$-fold standard downsampling in spacial domain is equivalent to
performing $s$-fold block downsampling in frequency domain:
${\bm{D}}_{\Downarrow s}={\bm{F}}_{\downarrow s}{\bm{D}}_{\downarrow
s}{\bm{F}}^{-1}$.

###### Proof.

Considering an arbitrary $d$-dimensional signal in frequency domain
$\hat{{\bm{x}}}[k],\ k=0,2,..,d-1$. Multiplying ${\bm{F}}_{\downarrow
s}{\bm{D}}_{\downarrow s}{\bm{F}}^{-1}$ to $\hat{{\bm{x}}}$ is equivalent to
letting $\hat{{\bm{x}}}$ go through the following linear system and obtain the
output $\hat{{\bm{y}}}$:

$\displaystyle{\bm{x}}[n]=\frac{1}{d}\sum_{k=0}^{d-1}\hat{{\bm{x}}}[k]e^{j\frac{2\pi}{d}kn}$
(81) $\displaystyle{\bm{x}}_{\downarrow s}[n]={\bm{x}}[ns]$ (82)
$\displaystyle\hat{{\bm{y}}}[k]=\sum_{n=0}^{d/s-1}{\bm{x}}_{\downarrow
s}[n]e^{-j\frac{2\pi}{d/s}kn}$ (83)

We now use $\hat{{\bm{x}}}$ to represent $\hat{{\bm{y}}}$:

$\displaystyle\hat{{\bm{y}}}[k]$
$\displaystyle=\sum_{n=0}^{d/s-1}\frac{1}{d}\sum_{k^{\prime}=0}^{d-1}\hat{{\bm{x}}}[k^{\prime}]e^{j\frac{2\pi}{d}k^{\prime}ns}e^{-j\frac{2\pi}{d/s}kn}$
(84)
$\displaystyle=\sum_{n=0}^{d/s-1}\frac{1}{d}e^{-j\frac{2\pi}{d/s}kn}\sum_{k^{\prime}=0}^{d-1}\hat{{\bm{x}}}[k^{\prime}]e^{j\frac{2\pi}{d/s}k^{\prime}n}$
(85)
$\displaystyle=\sum_{n=0}^{d/s-1}\frac{1}{d}e^{-j\frac{2\pi}{d/s}kn}(\sum_{k^{\prime}=0}^{d/s-1}+\sum_{k^{\prime}=d/s}^{2d/s-1}+\sum_{k^{\prime}=2d/s}^{3d/s-1}+...+\sum_{k^{\prime}=(s-1)d/s}^{sd/s-1})\hat{{\bm{x}}}[k^{\prime}]e^{j\frac{2\pi}{d/s}k^{\prime}n}$
(86)
$\displaystyle=\sum_{n=0}^{d/s-1}\frac{1}{d}e^{-j\frac{2\pi}{d/s}kn}\sum_{k^{\prime}=0}^{d/s-1}(\hat{{\bm{x}}}[k^{\prime}]+\hat{{\bm{x}}}[k^{\prime}+d/s]+...+\hat{{\bm{x}}}[k^{\prime}+(s-1)d/s])e^{j\frac{2\pi}{d/s}k^{\prime}n}$
(87)
$\displaystyle=\sum_{n=0}^{d/s-1}e^{-j\frac{2\pi}{d/s}kn}\frac{1}{d/s}\sum_{k^{\prime}=0}^{d/s-1}\frac{\hat{{\bm{x}}}[k^{\prime}]+\hat{{\bm{x}}}[k^{\prime}+d/s]+...+\hat{{\bm{x}}}[k^{\prime}+(s-1)d/s]}{s}e^{j\frac{2\pi}{d/s}k^{\prime}n}$
(88)
$\displaystyle=\frac{\hat{{\bm{x}}}[k]+\hat{{\bm{x}}}[k+d/s]+...+\hat{{\bm{x}}}[k+(s-1)d/s]}{s}$
(89)

where equation 87 is because $e^{j\frac{2\pi}{d/s}k^{\prime}n}$ has period of
$d/s$ in $k^{\prime}$ and equation 89 is because $d/s$-dimensional inverse
Fourier transform and Fourier transform are canceled out. From equation 89 we
have $\hat{{\bm{y}}}={\bm{D}}_{\Downarrow s}\hat{{\bm{x}}}$. So
${\bm{D}}_{\Downarrow s}\hat{{\bm{x}}}={\bm{F}}_{\downarrow
s}{\bm{D}}_{\downarrow s}{\bm{F}}^{-1}\hat{{\bm{x}}}$ for any
$\hat{{\bm{x}}}$, and consequently, ${\bm{D}}_{\Downarrow
s}={\bm{F}}_{\downarrow s}{\bm{D}}_{\downarrow s}{\bm{F}}^{-1}$. ∎

###  B.1 Closed-form Solutions to Vector-Jacobian Products in Type I Guidance

When using Type I guidance, we are required to approximate the likelihood
score, which can be approximated via vector-Jacobian product (Song et al.,
2023):

$\nabla_{{\bm{x}}_{t}}\log
p({\bm{y}}|{\bm{x}}_{t})\approx({\bm{v}}^{T}\frac{\partial\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]}{\partial{\bm{x}}_{t}})^{T},$
(90)

where
${\bm{v}}={\bm{A}}^{T}(\sigma^{2}{\bm{I}}+{\bm{A}}{\bm{\Sigma}}_{t}({\bm{x}}_{t}){\bm{A}}^{T})^{-1}({\bm{y}}-{\bm{A}}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}])$.
Below we provide closed-form solution to ${\bm{v}}$ under three common
degradation operators: i) inpainting, ii) debluring, and iii) super
resolution, then we can approximate the conditional posterior mean
$\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}}]$ based on proposition 1.

Inpainting. The observation model for image inpainting can be expressed as:

${\bm{y}}=\underbrace{{\bm{D}}_{{\bm{m}}}}_{{\bm{A}}}{\bm{x}}_{0}+{\bm{n}},$
(91)

and the closed-form solution to ${\bm{v}}$ in image inpainting is given by the
following:

${\bm{v}}=\frac{\tilde{{\bm{y}}}-{\bm{m}}\odot\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]}{\sigma^{2}+r_{t}^{2}({\bm{x}}_{t})}$
(92)

where
$\tilde{{\bm{y}}}={\bm{D}}_{{\bm{m}}}^{T}{\bm{y}}={\bm{m}}\odot({\bm{x}}_{0}+\tilde{{\bm{n}}}),\,\tilde{{\bm{n}}}\sim\mathcal{N}(\bm{0},{\bm{I}})$
is the zero-filling measurements that fills the masked region with zeros and
processes the exact same size to ${\bm{x}}_{0}$. In practice, the measurements
in inpainting are usually stored in the form of $\tilde{{\bm{y}}}$, while
${\bm{y}}$ used here is for mathematical convenient.

###### Proof.

$\displaystyle{\bm{v}}$
$\displaystyle={\bm{D}}_{\bm{m}}^{T}(\sigma^{2}{\bm{I}}+{\bm{D}}_{\bm{m}}r_{t}^{2}({\bm{x}}_{t}){\bm{I}}{\bm{D}}_{\bm{m}}^{T})^{-1}({\bm{y}}-{\bm{D}}_{\bm{m}}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}])$
(93)
$\displaystyle={\bm{D}}_{\bm{m}}^{T}(\sigma^{2}{\bm{I}}+r_{t}^{2}({\bm{x}}_{t}){\bm{D}}_{\bm{m}}{\bm{D}}_{\bm{m}}^{T})^{-1}({\bm{y}}-{\bm{D}}_{\bm{m}}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}])$
(94)
$\displaystyle={\bm{D}}_{\bm{m}}^{T}((\sigma^{2}+r_{t}^{2}({\bm{x}}_{t})){\bm{I}})^{-1}({\bm{y}}-{\bm{D}}_{\bm{m}}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}])$
(95)
$\displaystyle=\frac{{\bm{D}}_{\bm{m}}^{T}({\bm{y}}-{\bm{D}}_{\bm{m}}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}])}{\sigma^{2}+r_{t}^{2}({\bm{x}}_{t})}$
(96)
$\displaystyle=\frac{\tilde{{\bm{y}}}-{\bm{m}}\odot\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]}{\sigma^{2}+r_{t}^{2}({\bm{x}}_{t})}$
(97)

where equation 95 and equation 97 are because
${\bm{D}}_{\bm{m}}{\bm{D}}_{\bm{m}}^{T}={\bm{I}}$ and
${\bm{D}}_{\bm{m}}^{T}{\bm{D}}_{\bm{m}}=\mathrm{diag}[{\bm{m}}]$. ∎

Debluring. The observation model for image debluring can be expressed as:

${\bm{y}}={\bm{x}}_{0}\ast{\bm{k}}+{\bm{n}},$ (98)

where ${\bm{k}}$ is the blurring kernel and $\ast$ is convolution operator. By
assuming $\ast$ is a circular convolution operator, we can convert equation 98
to the canonical form ${\bm{y}}={\bm{A}}{\bm{x}}_{0}+{\bm{n}}$ by leveraging
the convolution property of Fourier transform:

${\bm{y}}=\underbrace{{\bm{F}}^{-1}\mathrm{diag}[\hat{{\bm{k}}}]{\bm{F}}}_{{\bm{A}}}{\bm{x}}_{0}+{\bm{n}},$
(99)

and the closed-form solution to ${\bm{v}}$ in image debluring is given by the
following:

${\bm{v}}={\bm{F}}^{-1}(\bar{\hat{{\bm{k}}}}\odot\frac{{\bm{F}}({\bm{y}}-{\bm{A}}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}])}{\sigma^{2}+r_{t}^{2}({\bm{x}}_{t})\bar{\hat{{\bm{k}}}}\odot\hat{{\bm{k}}}}).$
(100)

###### Proof.

Since ${\bm{A}}$ is a real matrix, we have
${\bm{A}}^{T}={\bm{A}}^{H}={\bm{F}}^{-1}\mathrm{diag}[\bar{\hat{{\bm{k}}}}]{\bm{F}}$,
then

$\displaystyle{\bm{v}}$
$\displaystyle={\bm{F}}^{-1}\mathrm{diag}[\bar{\hat{{\bm{k}}}}]{\bm{F}}(\sigma^{2}{\bm{I}}+{\bm{F}}^{-1}\mathrm{diag}[\hat{{\bm{k}}}]{\bm{F}}r_{t}^{2}({\bm{x}}_{t}){\bm{I}}{\bm{F}}^{-1}\mathrm{diag}[\bar{\hat{{\bm{k}}}}]{\bm{F}})^{-1}({\bm{y}}-{\bm{A}}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}])$
(101)
$\displaystyle={\bm{F}}^{-1}\mathrm{diag}[\bar{\hat{{\bm{k}}}}]{\bm{F}}(\sigma^{2}{\bm{I}}+r_{t}^{2}({\bm{x}}_{t}){\bm{F}}^{-1}\mathrm{diag}[\hat{{\bm{k}}}]\mathrm{diag}[\bar{\hat{{\bm{k}}}}]{\bm{F}})^{-1}({\bm{y}}-{\bm{A}}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}])$
(102)
$\displaystyle={\bm{F}}^{-1}\mathrm{diag}[\bar{\hat{{\bm{k}}}}]{\bm{F}}(\sigma^{2}{\bm{I}}+r_{t}^{2}({\bm{x}}_{t}){\bm{F}}^{-1}\mathrm{diag}[\hat{{\bm{k}}}\odot\bar{\hat{{\bm{k}}}}]{\bm{F}})^{-1}({\bm{y}}-{\bm{A}}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}])$
(103)
$\displaystyle={\bm{F}}^{-1}\mathrm{diag}[\bar{\hat{{\bm{k}}}}]{\bm{F}}({\bm{F}}^{-1}(\sigma^{2}{\bm{I}}+r_{t}^{2}({\bm{x}}_{t})\mathrm{diag}[\hat{{\bm{k}}}\odot\bar{\hat{{\bm{k}}}}]){\bm{F}})^{-1}({\bm{y}}-{\bm{A}}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}])$
(104)
$\displaystyle={\bm{F}}^{-1}\mathrm{diag}[\bar{\hat{{\bm{k}}}}]{\bm{F}}({\bm{F}}^{-1}\mathrm{diag}[\sigma^{2}+r_{t}^{2}({\bm{x}}_{t})\hat{{\bm{k}}}\odot\bar{\hat{{\bm{k}}}}]{\bm{F}})^{-1}({\bm{y}}-{\bm{A}}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}])$
(105)
$\displaystyle={\bm{F}}^{-1}\mathrm{diag}[\bar{\hat{{\bm{k}}}}]{\bm{F}}{\bm{F}}^{-1}\mathrm{diag}[\sigma^{2}+r_{t}^{2}({\bm{x}}_{t})\hat{{\bm{k}}}\odot\bar{\hat{{\bm{k}}}}]^{-1}{\bm{F}}({\bm{y}}-{\bm{A}}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}])$
(106)
$\displaystyle={\bm{F}}^{-1}(\bar{\hat{{\bm{k}}}}\odot\frac{{\bm{F}}({\bm{y}}-{\bm{A}}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}])}{\sigma^{2}+r_{t}^{2}({\bm{x}}_{t})\hat{{\bm{k}}}\odot\bar{\hat{{\bm{k}}}}})$
(107)

∎

Super resolution. According to (Zhang et al., 2020), the observation model for
image super resolution can be approximately expressed as:

${\bm{y}}=({\bm{x}}_{0}\ast{\bm{k}})_{\downarrow s}+{\bm{n}}$ (108)

By leveraging the convolution property of Fourier transform, we can convert
equation 108 to the canonical form ${\bm{y}}={\bm{A}}{\bm{x}}_{0}+{\bm{n}}$:

${\bm{y}}=\underbrace{{\bm{D}}_{\downarrow
s}{\bm{F}}^{-1}\mathrm{diag}[\hat{{\bm{k}}}]{\bm{F}}}_{{\bm{A}}}{\bm{x}}_{0}+{\bm{n}}$
(109)

and the closed-form solution to ${\bm{v}}$ in image super resolution is given
by the following:

${\bm{v}}={\bm{F}}^{-1}(\bar{\hat{{\bm{k}}}}\odot_{s}\frac{{\bm{F}}_{\downarrow_{s}}({\bm{y}}-{\bm{A}}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}])}{\sigma^{2}+r_{t}^{2}({\bm{x}}_{t})(\bar{\hat{{\bm{k}}}}\odot\hat{{\bm{k}}})_{\Downarrow_{s}}})$
(110)

where $\odot_{s}$ denotes block processing operator with element-wise
multiplication (Zhang et al., 2020).

###### Proof.

Since ${\bm{A}}$ is a real matrix, we have
${\bm{A}}^{T}={\bm{A}}^{H}={\bm{F}}^{-1}\mathrm{diag}[\bar{\hat{{\bm{k}}}}]{\bm{F}}{\bm{D}}_{\downarrow
s}^{T}$, then

$\displaystyle{\bm{v}}$
$\displaystyle={\bm{F}}^{-1}\mathrm{diag}[\bar{\hat{{\bm{k}}}}]{\bm{F}}{\bm{D}}_{\downarrow
s}^{T}(\sigma^{2}{\bm{I}}+{\bm{D}}_{\downarrow
s}{\bm{F}}^{-1}\mathrm{diag}[\hat{{\bm{k}}}]{\bm{F}}r_{t}^{2}({\bm{x}}_{t}){\bm{I}}{\bm{F}}^{-1}\mathrm{diag}[\bar{\hat{{\bm{k}}}}]{\bm{F}}{\bm{D}}_{\downarrow
s}^{T})^{-1}({\bm{y}}-{\bm{A}}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}])$ (111)
$\displaystyle={\bm{F}}^{-1}\mathrm{diag}[\bar{\hat{{\bm{k}}}}]{\bm{F}}{\bm{D}}_{\downarrow
s}^{T}(\sigma^{2}{\bm{I}}+r_{t}^{2}({\bm{x}}_{t}){\bm{D}}_{\downarrow
s}{\bm{F}}^{-1}\mathrm{diag}[\hat{{\bm{k}}}\odot\bar{\hat{{\bm{k}}}}]{\bm{F}}{\bm{D}}_{\downarrow
s}^{T})^{-1}({\bm{y}}-{\bm{A}}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}])$ (112)
$\displaystyle={\bm{F}}^{-1}\mathrm{diag}[\bar{\hat{{\bm{k}}}}]{\bm{F}}{\bm{D}}_{\downarrow
s}^{T}({\bm{D}}_{\downarrow
s}{\bm{F}}^{-1}(\sigma^{2}{\bm{I}}+r_{t}^{2}({\bm{x}}_{t})\mathrm{diag}[\hat{{\bm{k}}}\odot\bar{\hat{{\bm{k}}}}]){\bm{F}}{\bm{D}}_{\downarrow
s}^{T})^{-1}({\bm{y}}-{\bm{A}}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}])$ (113)
$\displaystyle={\bm{F}}^{-1}\mathrm{diag}[\bar{\hat{{\bm{k}}}}]{\bm{F}}{\bm{D}}_{\downarrow
s}^{T}({\bm{D}}_{\downarrow
s}{\bm{F}}^{-1}\mathrm{diag}[\sigma^{2}+r_{t}^{2}({\bm{x}}_{t})\hat{{\bm{k}}}\odot\bar{\hat{{\bm{k}}}}]{\bm{F}}{\bm{D}}_{\downarrow
s}^{T})^{-1}({\bm{y}}-{\bm{A}}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}])$ (114)

By Lemma 9, we have ${\bm{D}}_{\downarrow
s}{\bm{F}}^{-1}={\bm{F}}^{-1}_{\downarrow s}{\bm{D}}_{\Downarrow s}$. Taking
the hermitian transpose to both side and leveraging
${\bm{F}}^{-1}=\frac{1}{d}{\bm{F}}^{H}$ and ${\bm{F}}_{\downarrow
s}^{-1}=\frac{1}{d/s}{\bm{F}}_{\downarrow s}^{H}$, we also have
${\bm{F}}{\bm{D}}_{\downarrow s}^{T}=s{\bm{D}}_{\Downarrow
s}^{T}{\bm{F}}_{\downarrow s}$. So

$\displaystyle{\bm{v}}$
$\displaystyle={\bm{F}}^{-1}\mathrm{diag}[\bar{\hat{{\bm{k}}}}]s{\bm{D}}_{\Downarrow
s}^{T}{\bm{F}}_{\downarrow s}({\bm{F}}^{-1}_{\downarrow s}{\bm{D}}_{\Downarrow
s}\mathrm{diag}[\sigma^{2}+r_{t}^{2}({\bm{x}}_{t})\hat{{\bm{k}}}\odot\bar{\hat{{\bm{k}}}}]s{\bm{D}}_{\Downarrow
s}^{T}{\bm{F}}_{\downarrow
s})^{-1}({\bm{y}}-{\bm{A}}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}])$ (115)
$\displaystyle={\bm{F}}^{-1}\mathrm{diag}[\bar{\hat{{\bm{k}}}}]s{\bm{D}}_{\Downarrow
s}^{T}{\bm{F}}_{\downarrow s}({\bm{F}}^{-1}_{\downarrow
s}\mathrm{diag}[\sigma^{2}+r_{t}^{2}({\bm{x}}_{t})(\hat{{\bm{k}}}\odot\bar{\hat{{\bm{k}}}})_{\Downarrow
s}]{\bm{F}}_{\downarrow
s})^{-1}({\bm{y}}-{\bm{A}}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}])$ (116)
$\displaystyle={\bm{F}}^{-1}\mathrm{diag}[\bar{\hat{{\bm{k}}}}]s{\bm{D}}_{\Downarrow
s}^{T}{\bm{F}}_{\downarrow s}{\bm{F}}^{-1}_{\downarrow
s}\mathrm{diag}[\sigma^{2}+r_{t}^{2}({\bm{x}}_{t})(\hat{{\bm{k}}}\odot\bar{\hat{{\bm{k}}}})_{\Downarrow
s}]^{-1}{\bm{F}}_{\downarrow
s}({\bm{y}}-{\bm{A}}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}])$ (117)
$\displaystyle={\bm{F}}^{-1}(\bar{\hat{{\bm{k}}}}\odot s{\bm{D}}_{\Downarrow
s}^{T}\frac{{\bm{F}}_{\downarrow
s}({\bm{y}}-{\bm{A}}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}])}{\sigma^{2}+r_{t}^{2}({\bm{x}}_{t})(\hat{{\bm{k}}}\odot\bar{\hat{{\bm{k}}}})_{\Downarrow
s}})$ (118)
$\displaystyle={\bm{F}}^{-1}(\bar{\hat{{\bm{k}}}}\odot_{s}\frac{{\bm{F}}_{\downarrow_{s}}({\bm{y}}-{\bm{A}}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}])}{\sigma^{2}+r_{t}^{2}({\bm{x}}_{t})(\bar{\hat{{\bm{k}}}}\odot\hat{{\bm{k}}})_{\Downarrow_{s}}})$
(119)

∎

###  B.2 Closed-form Solutions to Proximal Problems in Type II Guidance

When using Type II guidance, we are required to solve the following auto-
weighted proximal problem:

$\mathbb{E}_{q}[{\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}}]=\arg\min_{{\bm{x}}_{0}}\lVert{\bm{y}}-{\bm{A}}{\bm{x}}_{0}\rVert^{2}+\sigma^{2}\lVert{\bm{x}}_{0}-\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]\rVert^{2}_{{\bm{\Sigma}}_{t}^{-1}({\bm{x}}_{t})}$
(120)

which has general closed-form solution given by

$\mathbb{E}_{q}[{\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}}]=({\bm{\Sigma}}_{t}({\bm{x}}_{t})^{-1}+\frac{1}{\sigma^{2}}{\bm{A}}^{T}{\bm{A}})^{-1}({\bm{\Sigma}}_{t}({\bm{x}}_{t})^{-1}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]+\frac{1}{\sigma^{2}}{\bm{A}}^{T}{\bm{y}})$
(121)

The derivation of equation 121 can be directly obtained via computing the mean
of $q_{t}({\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}})$ by using the Bayes’ theorem for
Gaussian variables (Bishop, 2006).

Below we provide closed-form solution to
$\mathbb{E}_{q}[{\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}}]$ under three common
degradation operators same as for Type I guidance, then we can directly
replace the denoising result in unconditional diffusion step with the closed-
form solution of $\mathbb{E}_{q}[{\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}}]$ to solve
target inverse problem. The closed-form results are obtained by borrowing the
results from (Zhu et al., 2023) and setting $\rho_{t}$ to
$\frac{\sigma^{2}}{r_{t}^{2}({\bm{x}}_{t})}$.

Inpainting. The observation model for image inpainting see equation 91, and
the closed-form solution to
$\mathbb{E}_{q}[{\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}}]$ is given as follows:

$\mathbb{E}_{q}[{\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}}]=\frac{\tilde{{\bm{y}}}+\rho_{t}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]}{\rho_{t}+{\bm{m}}},\
\ \rho_{t}=\frac{\sigma^{2}}{r_{t}^{2}({\bm{x}}_{t})}$ (122)

Debluring. The observation model for image debluring see equation 98 and
equation 99, and the closed-form solution to
$\mathbb{E}_{q}[{\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}}]$ is given as follows:

$\mathbb{E}_{q}[{\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}}]={\bm{F}}^{-1}\frac{\bar{\hat{{\bm{k}}}}\odot{\bm{F}}{\bm{y}}+\rho_{t}{\bm{F}}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]}{\bar{\hat{{\bm{k}}}}\odot\hat{{\bm{k}}}+\rho_{t}},\
\ \rho_{t}=\frac{\sigma^{2}}{r_{t}^{2}({\bm{x}}_{t})}$ (123)

Super resolution. The observation model for image super resolution see
equation 108 and equation 109, and the closed-form solution to
$\mathbb{E}_{q}[{\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}}]$ is given as follows:

$\mathbb{E}_{q}[{\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}}]={\bm{F}}^{-1}\frac{1}{\rho_{t}}({\bm{d}}-\bar{\hat{{\bm{k}}}}\odot_{s}\frac{(\hat{{\bm{k}}}\odot{\bm{d}})_{\Downarrow
s}}{(\bar{\hat{{\bm{k}}}}\odot\hat{{\bm{k}}})_{\Downarrow s}+\rho_{t}}),\ \
\rho_{t}=\frac{\sigma^{2}}{r_{t}^{2}({\bm{x}}_{t})}$ (124)

where
${\bm{d}}=\hat{{\bm{k}}}\odot({\bm{F}}{\bm{y}}_{\uparrow_{s}})+\rho_{t}{\bm{F}}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]$
and $\uparrow_{s}$ denotes the standard $s$-fold upsampling, i.e., upsampling
the spatial size by filling the new entries with zeros (Zhang et al., 2020).

##  Appendix C Using Conjugate Gradient Method for General Posterior
Covariance

For general posterior covariance, the closed-form solution for ${\bm{v}}$ or
$\mathbb{E}_{q}[{\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}}]$ is usually unavailable.
However, we can still approximate their solutions by using numerical solutions
of linear equations. Specifically, we can see that for Type II guidance,
$\mathbb{E}_{q}[{\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}}]$ satisfies the following
linear equation:

$({\bm{\Sigma}}_{t}({\bm{x}}_{t})^{-1}+\frac{1}{\sigma^{2}}{\bm{A}}^{T}{\bm{A}})\mathbb{E}_{q}[{\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}}]={\bm{\Sigma}}_{t}({\bm{x}}_{t})^{-1}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]+\frac{1}{\sigma^{2}}{\bm{A}}^{T}{\bm{y}}$
(125)

Since
${\bm{\Sigma}}_{t}({\bm{x}}_{t})^{-1}+\frac{1}{\sigma^{2}}{\bm{A}}^{T}{\bm{A}}$
is in fact the precision matrix of
$q_{t}({\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}})$, we know that it is symmetric and
positive-definite. Therefore, we can solve equation 125 using conjugate
gradient (CG) method (Hestenes et al., 1952).

Likewise, for Type I guidance, it easy to see that
$\sigma^{2}{\bm{I}}+{\bm{A}}{\bm{\Sigma}}_{t}({\bm{x}}_{t}){\bm{A}}^{T}$ is
symmetric and positive-definite. To approximate ${\bm{v}}$, we define a
temporal variable ${\bm{u}}$, and approximate it using CG to solve the
following linear equation:

$(\sigma^{2}{\bm{I}}+{\bm{A}}{\bm{\Sigma}}_{t}({\bm{x}}_{t}){\bm{A}}^{T}){\bm{u}}={\bm{y}}-{\bm{A}}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]$
(126)

Then ${\bm{v}}$ can be obtained as follows:

${\bm{v}}={\bm{A}}^{T}{\bm{u}}$ (127)

In our experiments, only the Convert-type posterior covariance (Section 4.1)
requires to use the CG method, as the posterior covariances used in prior
works and Analytic-type posterior covariance (Section 4.2) are isotropic that
possess closed-form solutions. We use the black box CG method implemented in
scipy.sparse.linalg.cg with tol=1e-4.

##  Appendix D Additional Experimental Details and Results

###  D.1 Additional Quantitative Results

In this section, we report the FID (Figure LABEL:fig:typeIIfid) and SSIM
(Figure LABEL:fig:typeIIssim) performance for supplementary results of Type II
guidance.

###  D.2 Converting Optimal Solutions Between Different Perturbation Kernels

Suppose we are given a family of optimal solutions $q_{t}$ defined by the
following MLE objectives for all $t\in[0,T]$:

$\max_{q_{t}}\mathbb{E}_{{\bm{x}}_{0}\sim
p({\bm{x}}_{0}),\epsilon\sim\mathcal{N}(\bm{0},{\bm{I}})}\log
q_{t}({\bm{x}}_{0}|{\bm{x}}_{t}=s_{t}({\bm{x}}_{0}+\sigma_{t}\epsilon))$ (128)

As can be seen, $q_{t}$ equals the optimal solutions to equation 12 when the
perturbation kernels $p_{t}({\bm{x}}_{t}|{\bm{x}}_{0})$ are set to
$\mathcal{N}({\bm{x}}_{t}|s_{t}{\bm{x}}_{0},s_{t}^{2}\sigma_{t}^{2}{\bm{I}})$.
Now, suppose we want to perform sampling based on the diffusion ODE (SDE)
under the perturbation kernels
$\mathcal{N}({\bm{x}}_{t}|\tilde{s}_{t}{\bm{x}}_{0},\tilde{s}_{t}^{2}\tilde{\sigma}_{t}^{2}{\bm{I}})$.
This means that we are required to provide optimal solutions $\tilde{q}_{t}$
defined by the following objectives for all $t\in[0,T]$:

$\max_{\tilde{q}_{t}}\mathbb{E}_{{\bm{x}}_{0}\sim
p({\bm{x}}_{0}),\epsilon\sim\mathcal{N}(\bm{0},{\bm{I}})}\log\tilde{q}_{t}({\bm{x}}_{0}|{\bm{x}}_{t}=\tilde{s}_{t}({\bm{x}}_{0}+\tilde{\sigma}_{t}\epsilon))$
(129)

The idea is that, $\tilde{q}_{t}$ can be directly represented by $q_{t}$, so
we do not require to perform re-training:

$\tilde{q}_{t}({\bm{x}}_{0}|{\bm{x}}_{t}={\bm{x}})=q_{t^{\prime}}({\bm{x}}_{0}|{\bm{x}}_{t^{\prime}}=\frac{s_{t^{\prime}}}{\tilde{s}_{t}}{\bm{x}}),\
\ \sigma_{t^{\prime}}=\tilde{\sigma}_{t}$ (130)

###### Proof.

We only need to show that $\tilde{q}_{t}$ defined by equation 128 and equation
130 are equivalent to $\tilde{q}_{t}$ defined by equation 129. From equation
130, we actually know that
$q_{t^{\prime}}({\bm{x}}_{0}|{\bm{x}}_{t^{\prime}}={\bm{x}})=\tilde{q}_{t}({\bm{x}}_{0}|{\bm{x}}_{t}=\frac{\tilde{s}_{t}}{s_{t^{\prime}}}{\bm{x}})$,
plug it in equation 128 at $t^{\prime}$ we have

$\max_{\tilde{q}_{t}}\mathbb{E}_{{\bm{x}}_{0}\sim
p({\bm{x}}_{0}),\epsilon\sim\mathcal{N}(\bm{0},{\bm{I}})}\log\tilde{q}_{t}({\bm{x}}_{0}|{\bm{x}}_{t}=\tilde{s}_{t}({\bm{x}}_{0}+\sigma_{t^{\prime}}\epsilon)),\
\ \sigma_{t^{\prime}}=\tilde{\sigma}_{t}$ (131)

which is equivalent to equation 129. ∎

For example, suppose we are given optimal solutions $q_{t}$ under DDPM
perturbation kernels, i.e.,
$p_{t}({\bm{x}}_{t}|{\bm{x}}_{0})=\mathcal{N}({\bm{x}}_{t}|\sqrt{\bar{\alpha}_{t}}{\bm{x}}_{0},\bar{\beta}_{t}{\bm{I}})$.
We aim to convert these solutions to optimal solutions $\tilde{q}_{t}$ under
perturbation kernels used in Section 2.2, i.e.,
$p_{t}({\bm{x}}_{t}|{\bm{x}}_{0})=\mathcal{N}({\bm{x}}_{t}|{\bm{x}}_{0},t^{2}{\bm{I}})$.
We can realize it using following two steps: (1) finding $t^{\prime}$ such
that $\sqrt{\frac{\bar{\beta}_{t^{\prime}}}{\bar{\alpha}_{t^{\prime}}}}=t$,
and (2) scaling the input of $q_{t^{\prime}}$ by the factor of
$\sqrt{\bar{\alpha}_{t^{\prime}}}$. Formally,

$\tilde{q}_{t}({\bm{x}}_{0}|{\bm{x}}_{t}={\bm{x}})=q_{t^{\prime}}({\bm{x}}_{0}|{\bm{x}}_{t^{\prime}}=\sqrt{\bar{\alpha}_{t^{\prime}}}{\bm{x}}),\
\ \sqrt{\frac{\bar{\beta}_{t^{\prime}}}{\bar{\alpha}_{t^{\prime}}}}=t$ (132)

###  D.3 Implementation of DDNM Solutions

Using DDNM solution (equation 11) for Type II guidance in noiseless inverse
problems requires to compute the pseudo-inverse ${\bm{A}}^{\dagger}$ for a
given linear operator ${\bm{A}}$. Below we present the implementations of the
pseudo-inverses used in our experiments.

Inpainting. We construct the DDNM solution
$\mathbb{E}_{q}[{\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}}]$ for the inpainting case
${\bm{A}}={\bm{D}}_{\bm{m}}$ as follows

$\mathbb{E}_{q}[{\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}}]=\tilde{{\bm{y}}}+(1-{\bm{m}})\odot\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]$
(133)

###### Proof.

Since ${\bm{D}}_{\bm{m}}{\bm{D}}_{\bm{m}}^{T}={\bm{I}}$ is non-singular, we
can directly obtain ${\bm{A}}^{\dagger}$ by
${\bm{A}}^{\dagger}={\bm{D}}_{{\bm{m}}}^{T}({\bm{D}}_{\bm{m}}{\bm{D}}_{\bm{m}}^{T})^{-1}={\bm{D}}_{{\bm{m}}}^{T}$.
Plug in equation 11, we have

$\displaystyle\mathbb{E}_{q}[{\bm{x}}_{0}|{\bm{x}}_{t},{\bm{y}}]$
$\displaystyle={\bm{D}}_{{\bm{m}}}^{T}{\bm{y}}+({\bm{I}}-{\bm{D}}_{{\bm{m}}}^{T}{\bm{D}}_{{\bm{m}}})\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]$
(134)
$\displaystyle=\tilde{{\bm{y}}}+({\bm{I}}-\mathrm{diag}[{\bm{m}}])\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]$
(135)
$\displaystyle=\tilde{{\bm{y}}}+\mathrm{diag}[1-{\bm{m}}]\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]$
(136)
$\displaystyle=\tilde{{\bm{y}}}+(1-{\bm{m}})\odot\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]$
(137)

∎

Debluring. We construct the pseudo inverse ${\bm{A}}^{\dagger}$ for the linear
operator in the debluring case
${\bm{A}}={\bm{F}}^{-1}\mathrm{diag}[\hat{{\bm{k}}}]{\bm{F}}$ as follows

${\bm{A}}^{\dagger}={\bm{F}}^{-1}\mathrm{diag}[\hat{{\bm{k}}}]^{\dagger}{\bm{F}}$
(138)

where $\mathrm{diag}[\hat{{\bm{k}}}]^{\dagger}$ is defined as
$\mathrm{diag}[\hat{{\bm{k}}}]^{\dagger}=\mathrm{diag}[[l_{1},l_{2},...]^{T}]$
with $l_{i}=0$ if $\hat{{\bm{k}}}_{i}=0$888Numerically, $|\hat{{\bm{k}}}_{i}|$
is always larger than zero. We threshold $\hat{{\bm{k}}}_{i}$ to zero when
$|\hat{{\bm{k}}}_{i}|<3\times 10^{-2}$ similar to (Wang et al., 2023)
otherwise $l_{i}=1/\hat{{\bm{k}}}_{i}$.

###### Proof.

It is easy to see that with the above construction, ${\bm{A}}^{\dagger}$
satisfies ${\bm{A}}{\bm{A}}^{\dagger}{\bm{A}}={\bm{A}}$. ∎

Super resolution. We directly use torch.nn.functional.interpolate in place of
${\bm{A}}^{\dagger}$ for the super resolution case.

###  D.4 Inference Time Analysis

The speed of inference of our method depends on whether Analytic or Convert is
chosen.

In the case of Analytic, our method does not impose any additional
computational cost for online inference compared to previous works. The
computational expense arises from the offline pre-calculation of
$\hat{r}_{t}^{2}$ using equation 24. Similar to (Bao et al., 2022b), a small
number of Monte Carlo samples is typically adequate for obtaining an accurate
$\hat{r}_{t}$. With a pre-trained unconditional model and training dataset, we
compute $\hat{r}_{t}$ offline and deploy it alongside the pre-trained
unconditional model. As a result, the online inference cost of our method is
the same as previous methods. In our paper, we initially calculate
$\hat{r}_{t}$ using a randomly selected $0.5\%$ of the training set and reuse
it for all our experiments. This process takes a couple of hours on a single
1080Ti GPU.

For Convert, the CG method can introduce additional computational cost.
However, we’ve observed that CG converges rapidly ($\sim 5$ steps when
initialize at
$D_{t}({\bm{x}}_{t})\approx\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]$) and is only
employed during the last few sampling steps. The additional computational cost
of CG is negligible compared to the overall inference time.

###  D.5 Correlation Reduction Using Orthonormal Transform

The main drawback of current plug-and-play posterior covariance optimization
methods is the oversimplified diagonal covariance structure, which assumes
pixels are uncorrelated. If we are allowed to train unconditional diffusion
models from scratch, we may wish to model the pixel correlation for potential
performance enhancement. In this section, we present a scalable way for
modeling correlation using orthonormal transform (OT), which is motivated by
transform coding (TC) in image compression literature (Ballé et al., 2020).
Specifically, we assume ${\bm{x}}_{0}$ is represented by some orthonormal
basis $\mathbf{\Psi}$, such that

${\bm{x}}_{0}=\mathbf{\Psi}{\bm{\theta}}_{0}$ (139)

By the property of covariance, we have

$\mathrm{Cov}[{\bm{x}}_{0}|{\bm{x}}_{t}]=\mathbf{\Psi}\mathrm{Cov}[{\bm{\theta}}_{0}|{\bm{x}}_{t}]\mathbf{\Psi}^{T}$
(140)

Assume the elements of ${\bm{\theta}}_{0}$ are mutually independent. Since
given condition ${\bm{x}}_{t}$ is equivalent to given ${\bm{\theta}}_{0}$ plus
isotrophic Gaussian noise, the elements of ${\bm{\theta}}_{0}|{\bm{x}}_{t}$
are also mutually independent. This motivates that, if the elements of
${\bm{\theta}}_{0}$ are more “closed to” mutually independent than
${\bm{x}}_{0}$ using proper basis,
$\mathrm{Cov}[{\bm{\theta}}_{0}|{\bm{x}}_{t}]$ could be better approximated by
a diagonal matrix than $\mathrm{Cov}[{\bm{x}}_{0}|{\bm{x}}_{t}]$. Thereby, we
parameterize the covariance of the variational Gaussian posterior
$q_{t}({\bm{x}}_{0}|{\bm{x}}_{t})$ as

${\bm{\Sigma}}_{t}({\bm{x}}_{t})=\mathbf{\Psi}\mathrm{diag}[{\bm{r}}_{t}^{2}({\bm{x}}_{t})]\mathbf{\Psi}^{T}$
(141)

The main benefit of using orthonormal basis for parameterizing
${\bm{\Sigma}}_{t}({\bm{x}}_{t})$ is that the number of values to predict in
${\bm{\Sigma}}_{t}({\bm{x}}_{t})$ is reduced from $d^{2}$ to $d$, which can be
a significant reduction for image data ($d>10^{5}$), and the training
procedure can be simply implemented as learning diagnoal Gaussian posterior in
the transform space (latent space).

Figure 8: Ground truths and model predictions. $\sigma_{t}$ is set to 1; Spatial Square Error refers to $({\bm{x}}_{0}-\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}])^{2}$, and DCT Square Error refers to $(\mathbf{\Psi}^{T}{\bm{x}}_{0}-\mathbf{\Psi}^{T}\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}])^{2}$, where $\mathbf{\Psi}$ is DCT basis; As can be seen, Spatial Square Error and DCT Square Error can be predicted by Spatial Variance and DCT Variance, respectively. Guidance | Method | Inpaint (Random) | Deblur (Gaussian) | Deblur (Motion) | Super resolution ($4\times$)
---|---|---|---|---|---
SSIM $\uparrow$ |  LPIPS $\downarrow$ |  FID $\downarrow$ |  SSIM $\uparrow$ |  LPIPS $\downarrow$ |  FID $\downarrow$ |  SSIM $\uparrow$ |  LPIPS $\downarrow$ |  FID $\downarrow$ |  SSIM $\uparrow$ |  LPIPS $\downarrow$ |  FID $\downarrow$
I |  DCTVar (Ours) | 0.9067 | 0.1128 | 40.42 | 0.7824 | 0.2077 | 66.80 | 0.7499 | 0.2407 | 78.23 | 0.7827 | 0.2171 | 70.24
SpatialVar (Ours) | 0.9271 | 0.0931 | 34.43 | 0.7872 | 0.2026 | 64.78 | 0.7513 | 0.2392 | 78.49 | 0.7835 | 0.2171 | 70.68
$\Pi\text{GDM}$ | 0.8759 | 0.1651 | 59.84 | 0.7844 | 0.2107 | 68.62 | 0.7497 | 0.2436 | 76.28 | 0.7793 | 0.2204 | 71.43
II |  DCTVar (Ours) | 0.8952 | 0.1160 | 37.93 | 0.7428 | 0.2290 | 79.05 | 0.6786 | 0.2876 | 108.74 | 0.7345 | 0.2473 | 92.89
SpatialVar (Ours) | 0.8446 | 0.1658 | 60.91 | 0.7299 | 0.2352 | 78.94 | 0.6405 | 0.3132 | 121.20 | 0.7165 | 0.2579 | 91.44
$\Pi\text{GDM}$ | 0.8597 | 0.1603 | 58.94 | 0.7340 | 0.2383 | 82.25 | 0.6377 | 0.3245 | 130.83 | 0.7124 | 0.2662 | 100.94
Table 5: Results for covariance with OT on FFHQ dataset. We use red when there
is no improvement over baseline.

We pretrain unconditional diffusion models based on k-diffusion code base for
predicting three values: 1) posterior mean
$\mathbb{E}[{\bm{x}}_{0}|{\bm{x}}_{t}]$, 2) ${\bm{r}}^{2}_{t}({\bm{x}}_{t})$
when $\mathbf{\Psi}$ is set to identity matrix, referred to as spatial
variance, and 3) ${\bm{r}}^{2}_{t}({\bm{x}}_{t})$ when $\mathbf{\Psi}$ is set
to DCT basis, referred to as DCT variance. Example model predictions see
Figure 8. To compare different posterior covariance for solving inverse
problems, similar to Section 5.2, the optimal posterior covariance is used
only at the low noise level (For spatial variance, $\sigma_{t}<0.2$; For DCT
variance, $\sigma_{t}<1$), while for high noise level, we use $\Pi\text{GDM}$
covariance.

Table 5 presents the evaluation results for different covariances.
$\Pi\text{GDM}$ is $\Pi\text{GDM}$ posterior covariance discussed in Section
3.1. DCTVar and SpatialVar are posterior covariance parameterized using
$\mathbf{\Psi}$ set to DCT basis and identity matrix respectively (equation
140). We observed that, a notable improvement with DCTVar is its performance
under high noise levels. Unlike SpatialVar, which performs poorly under high
noise levels and is only used when $\sigma_{t}<0.2$, DCTVar can be used for
significantly higher noise levels, and even for all sampling steps without
hurting the performance. This is because a Gaussian with DCTVar fits the
posterior more closely under high noise levels, while a Gaussian with
SpatialVar can only be effective under lower noise levels.

###  D.6 Qualitative Results

In this section, we present additional visual examples using FFHQ datasets to
showcase the capabilities of our method. Figure 9, 10, 11, and 12 serve as
complementary visual examples to Table 3. These illustrations consistently
demonstrate that our optimal covariances generate superior images compared to
previous heuristic covariances in Type I guidance. Furthermore, Figure
LABEL:fig:guidance_II_qual demonstrates how our methods can enhance the
robustness of DiffPIR against the hyper-parameter $\lambda$, taking random
inpainting as an example. We also display the notable performance improvement
of $\Pi\text{GDM}$ with the adaptive weight, i.e., complete version of
$\Pi\text{GDM}$, in Figure 15, 16, 17, and 18. Lastly, Figure
LABEL:fig:complete_dps illustrates that our methods can enhance the DPS
robustness to the hyper-parameter $\zeta$, using random inpainting as an
example.

Figure 9: Qualitative results for Table 3 on random inpainting. Figure 10:
Qualitative results for Table 3 on Gaussian debluring. Figure 11: Qualitative
results for Table 3 on motion debluring. Figure 12: Qualitative results for
Table 3 on super resolution ($4\times$). Figure 15: Qualitative results for
complete $\Pi\text{GDM}$ on random inpainting. Figure 16: Qualitative results
for complete $\Pi\text{GDM}$ on Gaussian debluring. Figure 17: Qualitative
results for complete $\Pi\text{GDM}$ on motion debluring. Figure 18:
Qualitative results for complete $\Pi\text{GDM}$ on super resolution
($4\times$).