Corresponding authors: ,

Monte Carlo guided Diffusion for Bayesian linear inverse problems

Ill-posed linear inverse problems arise frequently in various applications, from computational photography to medical imaging.
A recent line of research exploits Bayesian inference with informative priors to handle the ill-posedness of such problems.
Amongst such priors, score-based generative models (SGM) have recently been successfully applied to several different inverse problems.
In this study, we exploit the particular structure of the prior defined by the SGM to define a sequence of intermediate linear inverse problems. As the noise level decreases, the posteriors of these inverse problems get closer to the target posterior of the original inverse problem.
To sample from this sequence of posteriors, we propose the use of Sequential Monte Carlo (SMC) methods.
The proposed algorithm, , is shown to be theoretically grounded and we provide numerical simulations showing that it outperforms competing baselines when dealing with ill-posed inverse problems in a Bayesian setting.

§ INTRODUCTION

This paper is concerned with  linear inverse problems $\lmeas =\m{A} \lstate+ \noisestd \noise$, where $\lmeas \in \rset^\dimorig$ is a vector of indirect observations, $\lstate \in \rset^{\dimx}$ is the vector of unknowns,  $\m{A} \in \rset^{\dimorig \times \dimx}$ is the linear forward operator and $\noise \in \rset^{\dimorig}$ is an unknown noise vector. This general model is used throughout computational imaging, including various tomographic imaging applications such as common types of magnetic resonance imaging [44], X-ray computed tomography [14], radar imaging [6], and basic image restoration tasks such as deblurring, superresolution, and image inpainting [17].
The classical approach to solving linear inverse problems relies on prior knowledge about $\lstate$, such as its smoothness, sparseness in a dictionary, or its geometric properties.
These approaches attempt to estimate a $\widehat{\lstate}$ by minimizing a regularized inverse problem, $\hat{\lstate}= \operatorname{argmin}_{\lstate} \{ \| \lmeas - \m{A} \lstate \|^2 + \operatorname{Reg}(\lstate) \}$, where $\operatorname{Reg}$ is a regularization term that balances data fidelity and noise while enabling efficient computations.
However, a common difficulty in the regularized inverse problem is the selection of an appropriate regularizer, which has a decisive influence on the quality of the reconstruction.

Whereas regularized inverse problems continue to dominate the field, many alternative statistical formulations have been proposed; see [4, 22, 31] and the references therein - see [41] for a mathematical perspective. A main advantage of statistical approaches is that they allow for uncertainty quantification in the reconstructed solution; see [10].
The Bayes' formulation of the regularized inverse problem is based on considering the indirect measurement $\ormeas$, the state $\state$ and the noise $\noise$ as random variables, and to specify $p(\lmeas |\lstate)$ the likelihood (the conditional distribution of $\ormeas$ at $\state$)  and the prior $p(\lstate)$ (the distribution of the state). One can use Bayes' theorem to obtain the posterior distribution $p(\lstate|\lmeas) \propto p(\lmeas|\lstate) p(\lstate)$, where "$\propto$" means that the two sides are equal to each other up to a multiplicative constant that does not depend on $\lstate$.
Moreover, the use of an appropriate method for Bayesian inference allows the quantification of the uncertainty in the reconstructed solution $\lstate$. A variety of priors are available, including but not limited to Laplace [16], total variation (TV) [24] and mixture-of-Gaussians [15].
In the last decade, a variety of techniques have been proposed to design and train generative models capable of producing perceptually realistic samples from the original data, even in challenging high-dimensional data such as images or audio [28, 29, 18].
Denoising diffusion models have been shown to be particularly effective generative models in this context [36, 40, 37, 38, 3].
These models convert noise into the original data domain through a series of denoising steps.
A popular approach is to use a generic diffusion model that has been pre-trained, eliminating the need for re-training and making the process more efficient and versatile [43, 52].
Although this was not the main motivation for developing these models, they can of course be used as prior distributions in Bayesian inverse problems.
This simple observation has led to a new, fast-growing line of research on how linear inverse problems can benefit from the flexibility and expressive power of the recently introduced deep generative models; see [2, 46, 42, 25, 35, 54, 34].

§.§ Contributions

* We propose , a novel algorithm for sampling from the Bayesian posterior of Gaussian linear inverse problems with denoising diffusion model priors.  specifically exploits the structure of both the linear inverse problem and the denoising diffusion generative model to design an efficient SMC sampler.

* We establish under sensible assumptions that the empirical distribution of the samples produced by  converges to the target posterior when the number of particles goes to infinity. To the best of our knowledge,  is the first provably consistent algorithm for conditional sampling from the denoising diffusion posteriors.

* To evaluate the performance of , we perform numerical simulations on several examples for which the target posterior distribution is known. Simulation results support our theoretical results, i.e. the empirical distribution of samples from  converges to the target posterior distribution. This is not the case for the competing methods (using the same denoising diffusion generative priors) which are shown, when run with random initialization of the denoising diffusion, to generate a significant number of samples outside the support of the target posterior.
We also illustrate samples from  in imaging inverse problems.

Background and notations.
This section provides a concise overview of the diffusion model framework and notations used in this paper. We cover the elements that are important for understanding  our approach, and we recommend that readers refer to the original papers for complete details and derivations [36, 19, 40, 37].
A denoising diffusion model is a generative model consisting of a forward and a backward process.
The forward process involves sampling $\state_0 \sim \fwmarg{\mathrm{data}}{}$ from the data distribution, which is then converted to a sequence $\chunk{\state{}}{1}{n}$ of recursively corrupted versions of $\state_0$.
The backward process involves sampling $\state_n$ according to an easy-to-sample reference distribution on $\rset^{\dimx}$ and generating $\state_0 \in \rset^{\dimx}$ by a sequence of denoising steps.

Following [36, 37], the forward process can be chosen as a Markov chain with joint distribution

\begin{equation}
\label{eq:forward_def}
\fwmarg{0:n}{\chunk{\lstate}{0}{n}} = \fwmarg{\mathrm{data}}{\lstate_0} \textstyle\prod_{t=1}^{n} \fwtrans{t}{\lstate_{t-1}}{\lstate_{t}}, \quad \fwtrans{t}{\lstate_{t-1}}{\lstate_{t}}= \gauss( \lstate_t; (1 - \beta_t)^{1/2} \lstate_{t-1}, \beta_t \idm_{\dimx})\eqsp,
\end{equation}

where $\idm_{\dimx}$ is the identity matrix of size $\dimx$,  $\sequence{\beta}[t][\nset] \subset (0,1)$ is a non-increasing sequence and $\normpdf(\statevar{}; \mathbf{\mu}, \Sigma)$ is the p.d.f. of the Gaussian distribution with mean $\mathbf{\mu}$ and covariance matrix $\Sigma$ (assumed to be non-singular) evaluated at $\statevar{}$.
For all $t > 0$, set $\alphacumprod{t}{} = \prod_{\ell=1}^t (1 - \beta_\ell)$ with the convention $\alpha_0 = 1$. We have for all $0 \leq s < t \leq n$,

\begin{equation}
\fwtrans{t|s}{\lstate_{s}}{\lstate_t} \eqdef \textstyle \int \prod_{\ell=s+1}^{t} \fwtrans{\ell}{\lstate_{\ell-1}}{\lstate_\ell} \rmd \chunk{\lstate}{s+1}{t-1} = \normpdf(\lstate_{t}; (\alphacumprod{t}{} / \alphacumprod{s}{})^{1/2} \lstate_{s}, (1 - \alphacumprod{t}{} / \alphacumprod{s}{}) \idm_{\dimx}) \eqsp.
\end{equation}

For the standard choices of  $ \alphacumprod{t}{}$,  the sequence of distributions $(\fwmarg{t}{})_t$ converges weakly to the standard normal distribution as $t \to \infty$, which we chose as the reference distribution.   For the reverse process, [37, 38] introduce an inference distribution $\revbridge{1:n|0}{\lstate_0}{{\chunk{\lstate}{1}{n}}}$, depending on a sequence $\sequence{\sigma}[t][\nset]$ of hyperparameters satisfying
$\sigma_{t}^2 \in [0,1 - \alphacumprod{t-1}{}]$ for all $t \in \nset^*$, and defined as
\revbridge{n|0}{\lstate_0}{\lstate_n} \textstyle\prod_{t=n}^{2}
\revbridge{t-1|t, 0}{\lstate_t, \lstate_0}{\lstate_{t-1}} \eqsp,$ where $\revbridge{n|0}{\lstate_0}{\lstate_n}=\gauss\left(\lstate_n ; \alphacumprod{n}{}^{1/2} \lstate_0,\left(1-\alphacumprod{n}{}  \right) \Id \right)$ and
$\revbridge{t-1|t, 0}{\lstate_t,\lstate_0}{\lstate_{t-1}} =\gauss\left(\lstate_{t-1} ; \boldsymbol{\mu}_{t}(\lstate_0,\lstate_t), \sigma_{t}^2 \idm_{d} \right) \eqsp,$
$\boldsymbol{\mu}_{t}(\lstate_0,\lstate_t)= \alphacumprod{t-1}{} ^{1/2} \lstate_0+ (1- \alphacumprod{t-1}{} -\sigma^2 _t)^{1/2} (\lstate_t-\alphacumprod{t}{}^{1/2} \lstate_0) \big/ (1-\alphacumprod{t}{})^{1/2} \eqsp.$
For $t=n-1,\dots,1$, we define by backward induction the sequence $
\revbridge{t|0}{\lstate_0}{\lstate_t} %=  \int \revbridge{1:n|0}{\lstate_0}{\lstate_{1:n}} \rmd \chunk{\lstate}{1}{t-1} \rmd \chunk{\lstate}{t+1}{n}
=  \int \revbridge{t|t+1,0}{\lstate_{t+1}, \lstate_0}{\lstate_{t}} \revbridge{t+1|0}{\lstate_0}{\lstate_{t+1}} \rmd \lstate_{t+1}.$
It is shown in <cit.> that for all $t \in [1:n]$, the  distributions of the forward and inference process conditioned on the initial state coincide, i.e. that
$\revbridge{t|0}{\lstate_0}{\lstate_t} =  \fwtrans{t|0}{\lstate_0}{\lstate_t}.$
The backward process is derived from the inference distribution by replacing, for each $t \in [2:n]$, $\lstate_0$ in the definition $\revbridge{t-1|t, 0}{\lstate_t, \lstate_0}{\lstate_{t-1}}$ with a prediction where $\xpred{0|t}[\theta](\lstate_t) \eqdef \alphacumprod{t}{}^{-1/2} \left(\lstate_t - (1 - \alphacumprod{t}{})^{1/2}\epsprop{\lstate_t, t}{\theta}\right)$

where $\epsprop{\lstate,t}{\theta}$ is typically a neural network parameterized by $\theta$.
More formally, the backward distribution is defined as
\oldbwmargparam{0:n}{\lstate_{0:n}}{\param} =\oldbwmargparam{n}{\lstate_n}{} \textstyle\prod_{t=0}^{n+1} \bwtrans{\theta}{t}{\lstate_{t+1}}{\lstate_{t}}\eqsp,
$ where
$\oldbwmargparam{n}{\lstate_n}{}  =  \normpdf{\left(\lstate_n; 0_\dimx, \Id_{\dimx} \right)}$ and
for all $t \in [1:n-1]$,

\begin{equation}
\label{eq:bwker:defn}
\bwtrans{\theta}{t}{\lstate_{t+1}}{\lstate_{t}} \eqdef \revbridge{t|t+1, 0}{\lstate_{t+1}, \xpred{0|t+1}[\theta](\lstate_{t+1})}{\lstate_{t}} = \normpdf(\lstate_t, \bwmean^\param _{t+1}(\lstate_{t+1}), \sigma^2 _{t+1} \Id_{\dimx})\eqsp,
\end{equation}

where $\bwmean_{t+1}(\lstate_{t+1}) \eqdef \bridgemean(\xpred{0|t+1}[\theta](\lstate_{t+1}), \lstate_{t+1})$ and $0_\dimx$ is the null vector of size $\dimx$. At step $0$, we set $\bwtrans{}{0}{\lstate_1}{\lstate_0} \eqdef \normpdf(\lstate_0; \xpred{0|1}[\theta](\lstate_{1}), \sigma^2 _1 \Id_{\dimx})$.
The parameter  $\theta$ is obtained <cit.> by solving the following optimization problem:

\begin{equation}
\label{eq:optimal-denoising}
\theta_{*} \in \textstyle \argmin_{\theta} \sum_{t=1}^{n} (2 \dimx \sigma^2_t \alpha_t)^{-1} \int \|\epsilon - \epsprop{\sqrt{\alpha_t}x_0 + \sqrt{1 - \alpha_t}\epsilon, t}{\theta}\|_2^2 \normpdf(\epsilon; 0_{\dimx}, \idm_{\dimx}) \fwmarg{\operatorname{data}}{\rmd \lstate_0} \rmd \epsilon \eqsp.
\end{equation}

Thus, $\epsprop{\state_t,t}{\theta_*}$ might be seen as the predictor of the noise added to $\state_0$ to obtain $\state_t$ (in the forward pass) and justifies the ”prediction” terminology.
The time $0$ marginal $\smash{\oldbwmargparam{0}{\lstate_0}{\param_*} = \int \oldbwmargparam{0:n}{\lstate_{0:n}}{\param_*} \rmd \lstate_{1:n}}$
which we will refer to as the prior  is used as an approximation of $\fwmarg{\mathrm{data}}{}$ and the time $s$ marginal is $\smash{\oldbwmargparam{s}{\lstate_s}{\param_*} = \int \oldbwmargparam{0:n}{\lstate_{0:n}}{\param_*} \rmd \lstate_{1:s-1} \rmd \lstate_{s+1:n}}$.
In the rest of the paper we drop the dependence on the parameter $\param_*$. We define for all $v \in \R^{\ell},w \in \R^{k}$, the concatenation operator $\cat{v}{w} = [v^T, w^T]^T \in \R^{\ell + k}$. For $i \in [1:\ell]$, we let $\vecidx{v}{}{i}$ the $i$-th coordinate of $v$.

Related works.
The subject of Bayesian problems is very vast, and it is impossible to discuss here all the results obtained in this very rich literature.
One of such domains is image restoration problems, such as deblurring, denoising inpainting, which are challenging problems in computer vision that involves restoring a partially observed degraded image.
Deep learning techniques are widely used for this task [2, 49, 48, 46]  with many of them relying on auto-encoders, VAEs [23, 32, 53], GANs [49, 51], or autoregressive transformers [50, 45].
In what follows, we focus on methods based on denoising diffusion that has recently emerged as a way to produce high-quality realistic samples from the original data distribution on par with the best GANs in terms of image and audio generation, without the intricacies of adversarial training; see [36, 40, 39].
Diffusion-based approaches do not require specific training for degradation types, making them much more versatile and computationally efficient. In [39], noisy linear inverse problems are proposed to be solved by diffusing the degraded observation forward, leading to intermediate observations $\{ \lmeas_s \}_{s = 0} ^n$, and then running a modified backward process that promotes consistency with $\lmeas_s$ at each step $s$.
The Denoising-Diffusion-Restoration model () [26] also modifies the backward process so that the unobserved part of the state follows the backward process while the observed part is obtained as a noisy weighted sum between the noisy observation and the prediction of the state. As observed by [30],  is very efficient, but the simple blending used occasionally causes inconsistency in the restoration process.
[8] considers a backward process targeting the posterior.
approximates the score of the posterior using the Tweedie formula, which incorporates the learned score of the prior.
The approximation error is quantified and shown to decrease when the noise level is large, i.e., when the posterior is close to the prior distribution. As shown in <Ref> with a very simple example, neither  nor  can be used to sample the target posterior and therefore do not solve the Bayesian recovery problem (even if we run  and  several time with independent initializations). Indeed, we show that  and  produce samples under the "prior" distribution (which is generally captured very well by the denoising diffusion model), but which are not consistent with the observations (many samples land in areas with very low likelihood).
In [43] the authors introduce , a Sequential Monte Carlo-based denoising diffusion model that aims at solving specifically the inpainting problem.  produces a particle approximation of the conditional distribution of the non observed part of the state conditionally on a forward-diffused trajectory of the observation. The resulting particle approximation is shown to converge to the true posterior of the SGM under the assumption that the joint laws of the forward and backward processes coincide, which fails to be true in realistic setting.

In comparison with ,  is a versatile approach that solves any Bayesian linear inverse problem while being consistent under practically no assumption. In parallel to our work, [47] also developed a similar SMC based methodology but with a different proposal kernel.

§ THE  ALGORITHM

In this section we present our methodology for the inpainting problem (<ref>), both with noise and without noise.
The more general case is treated in <Ref>.
Let $\dimorig \in [1:\dimx - 1]$.
In what follows we denote the $\dimorig$ top coordinates of a vector $\lstate \in \rset^{\dimx}$ by $\ltopstate$ and the remaining coordinates by $\lbottomstate$, so that $\lstate = \cat{\ltopstate}{\lbottomstate}$. The inpainting problem is defined as

\begin{equation}
\label{eq:inpainting}
\ormeas =  \topstate + \noisestd \noise \eqsp, \quad \varepsilon \sim \ndist{0}{\idm_{\dimorig}} \eqsp, \quad \sigma \geq 0 \eqsp,
\end{equation}

where $\topstate$ are the first $\dimorig$ coordinates of a random variable $\state \sim \bwmarg{0}{}$. The goal is then to recover the law of the complete state $\state$ given a realisation $\ltrmeas$ of the incomplete observation $\trmeas$ and the model (<ref>).

Noiseless case.

We begin by the case $\noisestd=0$. As the first $\dimorig$ coordinates are observed exactly, we aim at infering the remaining coordinates of $\state$,
which correspond to $\bottomstate$.
As such, given an observation $\lmeas$, we aim at sampling from the posterior $\smash{\filter{0}{\lmeas}{\lbottomstate_0}} \propto \bwmarg{0}{\cat{\lmeas}{\lbottomstate_0}}$ with integral form

\begin{equation}
\label{eq:posterior:integral_form}
\filter{0}{\lmeas}{\lbottomstate_0} \propto \textstyle \int \bwmarg{n}{\lstate_n} \left\{ \prod_{s = 1}^{n-1} \bwtrans{}{s}{\lstate_{s+1}}{\lstate_s} \right\} \bwtrans{}{0}{\lstate_1}{\cat{\lmeas}{\lbottomstate_0}} \rmd \lstate_{1:n} \eqsp.
\end{equation}

To solve this problem, we propose to use SMC algorithms [13, 5, 7], where a set of $N$ random samples, referred to as particles, is iteratively updated to approximate the posterior distribution. The updates involve, at iteration $s$, selecting promising particles from the pool of particles $\pchunk{\particle_{s+1}}{1}{N} = (\particle^1 _{s+1}, \dotsc, \particle^N _{s+1})$ based on a weight function $\uweight{}{s}$, and then apply a Markov transition $\bwopt{s}{\lmeas}{}{}$ to obtain the samples $\pchunk{\particle_{s}}{1}{N}$.
The transition $\bwopt{s}{\lmeas}{\lstate_{s+1}}{\lstate_s}$ is designed to follow the backward process while guiding the $\dimorig$ top coordinates of the pool of particles $\pchunk{\particle_s}{1}{N}$ towards the measurement $\lmeas$.

Note that under the backward dynamics (<ref>), $\topstate_t$ and $\bottomstate_t$ are independent conditionally on $\state_{t+1}$ with transition kernels respectively
$\tbwtrans{}{t}{\lstate_{t+1}}{\ltopstate_{t}}  \eqdef \normpdf(\ltopstate_{t}; \tbwmean_{t+1}(\lstate_{t+1}), \sigma^2 _{t+1} \Id_{\dimorig})$ and $
\bbwtrans{}{t}{\lstate_{t+1}}{\lbottomstate_{t}}  \eqdef \normpdf(\lbottomstate_{t}; \bbwmean_{t+1}(\lstate_{t+1}), \sigma^2 _{t+1} \Id_{\dimx - \dimorig})$
where $\tbwmean_{t+1}(\lstate_{t+1}) \in \rset^{\dimorig}$ and $\bbwmean_{t+1}(\lstate_{t+1}) \in \rset^{\dimx - \dimorig}$ are such that $\bwmean_{t+1}(\lstate_{t+1}) = \cat{\tbwmean_{t+1}(\lstate_{t+1})}{\bbwmean_{t+1}(\lstate_{t+1})}$ and the above kernels satisfy $\bwtrans{}{t}{\lstate_{s+1}}{\lstate_s} = \tbwtrans{}{t}{\lstate_{t+1}}{\ltopstate_{t}} \bbwtrans{}{t}{\lstate_{t+1}}{\lbottomstate_{t}}$. We consider the following proposal kernels for $t \in [1:n]$,

\begin{equation}
\label{eq:optkernel}
\bwopt{s}{\lmeas}{\lstate_{t+1}}{\lstate_t} \propto \bwtrans{}{t}{\lstate_{t+1}}{\lstate_t} \tfwtrans{t|0}{\lmeas}{\ltopstate_t} \eqsp, \quad \mathrm{where} \quad \tfwtrans{t|0}{\lmeas}{\ltopstate_t} \eqdef \normpdf(\ltopstate_t; \alphacumprod{t}{}^{1/2} \lmeas, (1 - \alphacumprod{t}{}) \Id_{\dimorig}) \eqsp,
\end{equation}

and $\bwopt{n}{\lmeas}{}{}(\lstate_n) \propto \bwmarg{n}{\lstate_n} \tfwtrans{n|0}{\lmeas}{\ltopstate_n}$. For the final step, we define $\bbwopt{0}{\lmeas}{\lstate_1}{\lbottomstate_0} = \bbwtrans{}{0}{\lstate_1}{\lbottomstate_0}$. Using standard Gaussian conjugation formulas, we obtain

\begin{align*}
\bwopt{t}{\lmeas}{\lstate_{t+1}}{\lstate_t} & = \bbwtrans{}{t}{\lstate_{t+1}}{\lbottomstate_t} \cdot \normpdf\left(\ltopstate_t; \mathsf{K}_t \alpha^{1/2} _t \lmeas + (1 - \mathsf{K}_t) \tbwmean_{t+1}(\lstate_{t+1}), (1 - \alphacumprod{t}{}) \mathsf{K}_t \cdot \Id_{\dimorig} \right) \eqsp,\\
\bwmarg{n}{\lstate_n} & = \normpdf(\lbottomstate_n; 0_{\dimx - \dimorig}, \Id_{\dimx - \dimorig})  \cdot \normpdf(\ltopstate_n; \mathsf{K}_n \alphacumprod{n}{}^{1/2} \lmeas, (1 - \alphacumprod{n}{}) \mathsf{K}_n \cdot \Id_{\dimorig})
\end{align*}

where $\mathsf{K}_{t} \eqdef \sigma^2 _{t+1} \big/ (\sigma^2 _{t+1} + 1 - \alpha_t)$. For this procedure to target the posterior $\smash{\filter{0}{\lmeas}{}}$, the weight function $\uweight{}{s}$ is chosen as follows; we set $
\uweight{}{n-1}(\lstate_n) \eqdef \int \bwtrans{}{n-1}{\lstate_{n}}{\lstate_{n-1}} \tfwtrans{n-1|0}{\lmeas}{\ltopstate_{n-1}} \rmd \lstate_{n-1}
= \normpdf\left(\alpha_{n-1}^{1/2} \lmeas; \tbwmean_{n}(\lstate_n), \sigma_{n}^2 + 1 - \alpha_s\right)$
and for $\smash{t \in [1:n-2]}$,

\begin{equation}
\label{eq:weightfunc}
\uweight{}{t}(\lstate_{t+1}) \eqdef  \frac{\int \tbwtrans{}{t}{\lstate_{t+1}}{\ltopstate_{t}} \tfwtrans{t|0}{\lmeas}{\ltopstate_{t}} \rmd \lstate_{t}}{\tfwtrans{t+1|0}{\lmeas}{\ltopstate_{t+1}}}
= \frac{\normpdf\left(\alpha_{t}^{1/2} \lmeas; \tbwmean_{t+1}(\lstate_{s+1}), (\sigma_{t+1}^2 + 1 - \alpha_t) \Id_{\dimorig}\right)}{\normpdf\left(\alpha_{t+1}^{1/2}\lmeas; \ltopstate_{t+1}, (1 - \alpha_{t+1}) \Id_{\dimorig} \right)} \eqsp.
\end{equation}

For the final step, we set $\uweight{}{0}(\lstate_1) \eqdef \tbwtrans{}{0}{\ltopstate_1}{\lmeas} \big/ \tfwtrans{1|0}{\lmeas}{\ltopstate_1}$.
The overall SMC algorithm targeting $\filter{0}{\lmeas}{}$ using the instrumental kernel (<ref>) and weight function (<ref>) is summarized in <Ref>.

($\sigma = 0$)

Number of particles $N$
$\pchunk{\particle _0}{1}{N}$
Operations involving index $i$ are repeated for each $i \in [1:N]$
$\overline{z}^i _{n} \sim \normpdf(\mathbf{0}_{\dimorig}, \idm_{\dimorig}), \quad \underline{z}^i _{n} \sim \normpdf(\mathbf{0}_{\dimx - \dimorig}, \idm_{\dimx - \dimorig}), \quad \tparticle^i _n = \mathsf{K}_n \alphacumprod{n}{}^{1/2} \lmeas + (1 - \alphacumprod{n}{}) \mathsf{K} _n \overline{z}^i _n , \quad \particle^i _n = \cat{\tparticle^i _n}{\underline{z}^i _n}$;

$s \gets n-1:0$
$s = n-1$
$\uweight{}{n-1}(\particle^i _n) = \normpdf\big(\alphacumprod{n}{}^{1/2} \lmeas; \topvec{\bwmean}_{n}(\particle^i _n), 2 - \alphacumprod{n}{} \big)$;

$\uweight{}{s}(\particle^i _{s+1}) = \normpdf\left(\alphacumprod{s}{}^{1/2} \lmeas; \topvec{\bwmean}_{s+1}(\particle^i _{s+1}), \sigma_{s+1}^2 + 1 - \alphacumprod{s}{}\right) \big/ \normpdf\left(\alphacumprod{s+1}{}^{1/2}\lmeas; \tparticle^i _{s+1}, 1 - \alphacumprod{s+1}{} \right)$;

$\anc{i}{s+1} \sim \categorical\big(\{ \uweight{}{s}(\particle^{j} _{s+1}) / \sum_{k = 1}^N \uweight{}{s}(\particle^{k} _{s+1})\} _{j = 1} ^N \big), \quad \overline{z}^i _{s} \sim \normpdf(\mathbf{0}_{\dimorig}, \idm_{\dimorig}), \quad \underline{z}^i _{s} \sim \normpdf(\mathbf{0}_{\dimx - \dimorig}, \idm_{\dimx - \dimorig})$;

$\tparticle^i _s = \mathsf{K}_s \alphacumprod{s}{}^{1/2} \lmeas + (1 - \mathsf{K}_s) \topvec{\bwmean}_{s+1}(\particle^i _{s+1}) + (1 - \alpha_s)^{1/2} \mathsf{K}^{1/2} _s \overline{z}^i _{s}, \quad \bparticle^{i}_{s} = \botvec{\bwmean}_{s+1}(\particle^i _{s+1}) + \sigma _{s+1} \underline{z}^i _{s}$;

Set $\particle^i _s = \cat{\tparticle^i _s}{\bparticle^i _s}$;

We now provide a justification to <Ref>. Let $\{ \pot{s}{\lmeas}{} \}_{s = 1} ^{n}$ be a sequence of positive functions. Consider the sequence of distributions $\left\{\smash{\filter{s}{\lmeas}{}}\right\}_{s=1}^{n}$ defined as follows; $\filter{n}{\lmeas}{\lstate_n} \propto \bwmarg{n}{\lstate_n} \pot{n}{\lmeas}{\lstate_n}$ and for $t \in [1:n-1]$

\begin{equation}
\filter{t}{\lmeas}{\lstate_{t}}
\propto \textstyle \int {\pot{t+1}{\lmeas}{\lstate_{t+1}}}^{-1}{\pot{t}{\lmeas}{\lstate_t}} \bwtrans{}{t}{\lstate_{t+1}}{\lstate_{t}}  \filter{t+1}{\lmeas}{\rmd \lstate_{t+1}} \eqsp,  \label{eq:FKrecursion}
\end{equation}

By construction, the time $t$ marginal (<ref>) is $\filter{t}{\lmeas}{\lstate_t} \propto \bwmarg{t}{\lstate_t} \pot{t}{\lmeas}{\lstate_t}$ for all $t \in [1:n]$.
Then, using $\smash{\filter{1}{\lmeas}{}}$ and (<ref>), we have that $\filter{0}{\lmeas}{\lbottomstate_0} \propto  \int  {\pot{1}{\lmeas}{\lstate_1}}^{-1}{\tbwtrans{}{0}{\ltopstate_1}{\lmeas}} \bbwtrans{}{0}{\lstate_1}{\lbottomstate_0}  \filter{1}{\lmeas}{\rmd \lstate_1}$.

The recursion (<ref>) suggests a way of obtaining a particle approximation of $\filter{0}{\lmeas}{}$; by sequentially approximating each $\filter{t}{\lmeas}{}$ we can effectively derive a particle approximation of the posterior.
To construct the intermediate particle approximations we use the framework of auxiliary particle filters (APF) [33].
We focus on the case $\pot{t}{\lmeas}{\lstate_t} = \tfwtrans{t|0}{\lmeas}{\ltopstate_t}$ which corresponds to <Ref>. The initial particle approximation $\filter{n}{\lmeas}{}$ is obtained by drawing $N$ i.i.d. samples $\pchunk{\particle_{n}}{1}{N}$ from $\bwopt{n}{\lmeas}{}{}$ and setting $\filter{n}{N}{} = N^{-1} \sum_{i = 1}^N \delta_{\particle^{i} _{n}}$ where $\delta_\particle$ is the Dirac mass at $\xi$.
Assume that the empirical approximation of $\filter{t+1}{\lmeas}{}$ is $\filter{t+1}{N}{} = N^{-1} \sum_{i = 1}^N \delta_{\particle^{i} _{t+1}}\eqsp,$ where $\pchunk{\particle_{t+1}}{1}{N}$ are $N$ random variables.
Substituting $\smash{\filter{t+1}{N}{}}$ into the recursion (<ref>) and introducing the instrumental kernel (<ref>), we obtain the mixture

\begin{equation}
\label{eq:particleapprox:intermediate}
\widehat\phi^N _t (\lstate_t) = \textstyle  \sum_{i = 1}^N {\uweight{}{t}(\particle^{i} _{t+1})}  \bwopt{t}{\lmeas}{\particle^i _{t+1}}{\lstate_t} \big / \sum_{j = 1}^N  \uweight{}{t}(\particle^{j} _{t+1}) \eqsp.
\end{equation}

Then, a particle approximation of (<ref>) is obtained by
sampling $N$ conditionally i.i.d. ancestor indices  $\anc{1:N}{t+1} \iid \categorical(\{ \uweight{}{t}(\particle^{i} _{t+1}) / \sum_{j = 1}^N \uweight{}{t}(\particle^{j} _{t+1})\} _{i = 1} ^N)$,
and then propagating each ancestor particle $\particle^{I^i _{t+1}} _{t+1}$ according to the instrumental kernel (<ref>).
The final particle approximation is given by
$\filter{0}{N}{} = N^{-1} \sum_{i = 1}^N \delta_{\bparticle^{i} _{0}},$ where $\smash{\bparticle^{i} _{0} \sim \bbwtrans{}{0}{\particle^{\anc{i}{0}} _1}{\cdot}}$, $\smash{\anc{i}{0} \sim \categorical(\{ \uweight{}{0}(\particle^{k} _{1}) \big/ \sum_{j = 1}^N \uweight{}{0}(\particle^j _1) \} _{k = 1} ^N)}$.
The sequence of distributions $\{ \bwmarg{t}{} \}_{t=0}^{n}$ approximating the marginals of the forward process initialized at $\bwmarg{0}{}$ defines a path that bridges between $\bwmarg{n}{}$ and the prior $\bwmarg{0}{}$ such that the discrepancy between $\bwmarg{t}{}$ and $\bwmarg{t+1}{}$ is small. SMC samplers based on this path are robust to multi-modality and offer an interesting alternative to the geometric and tempering paths traditionally used in the SMC literature, see [9]. Our proposals $\filter{t}{\lmeas}{\lstate_t} \propto \bwmarg{t}{\lstate_t} \tfwtrans{t|0}{\lmeas}{\ltopstate_t}$ inherit the behavior of $\{ \bwmarg{t}{} \}_{t \in \nset}$ and bridge the initial distribution $\filter{n}{\lmeas}{}$ and posterior $\filter{0}{\lmeas}{}$. Indeed, as $\lmeas$ is a noiseless observation of $\state_0 \sim \bwmarg{0}{}$, we may consider $\smash{\alphacumprod{t}{}^{1/2} \lmeas + (1 - \alphacumprod{t}{})^{1/2} \noise_t}$, with $\noise_t \sim \normpdf(\mathbf{0}_{\dimorig}, \Id_{\dimorig})$, as a noisy observation of $X_t \sim \bwmarg{t}{}$ and thus, $\filter{t}{\lmeas}{}$ is the associated posterior.  We illustrate this intuition by considering the following Gaussian mixture (GM) example. We assume that $\bwmarg{0}{\lstate_0} = \sum_{i = 1}^M w_i \cdot \normpdf(\lstate_0; \mu_i, \Id_{\dimx})$ where $M > 1$ and $\{ w_i \}_{i = 1}^M$ are drawn uniformly on the simplex. The marginals of the forward process are available in closed form and are given by $\bwmarg{t}{\lstate_t} = \sum_{i = 1}^M w_i \cdot \normpdf(\lstate_t; \alphacumprod{t}{}^{1/2} \mu_i, \Id_{\dimx})$, which shows that the discrepancy between $\bwmarg{t}{}$ and $\bwmarg{t+1}{}$ is small as long as $\alphacumprod{t}{}^{1/2} - \alphacumprod{t+1}{}^{1/2} \approx 0$. The posteriors $\{ \filter{t}{\lmeas}{} \}_{t \in [0:n]}$ are also available in closed form and displayed in <Ref>, which illustrates that our choice of potentials ensures that the discrepancy between consecutive posteriors is small.

$t = 450$      $t = 100$      $t = 80$     $t = 70$     $t=50$     $t=20$     $t=15$     $t=5$

Display of samples from $\protect\filter{t}{\lmeas}{\lstate_t} \propto \protect\bwmarg{t}{\lstate_t} \protect\tfwtrans{t|0}{\lmeas}{\ltopstate_t}$ for the GM prior. In yellow the samples from $\protect\filter{t}{\lmeas}{}$, in purple those from the prior and in light blue those from the posterior $\protect\filter{0}{\lmeas}{}.$ $n$ is set to $500$.

The idea of using the forward diffused observation to guide the observed part of the state, as we do here through $\tfwtrans{t}{\lmeas}{\ltopstate_t}$, has been exploited in prior works but in a different way. For instance, in [40, 39] the observed part of the state is directly replaced by the forward noisy observation and, as it has been noted [43], this introduces an irreducible bias. Instead,  weights the backward process by the density of the forward one conditioned on $\lmeas$, resulting in a natural and consistent algorithm.

We now establish the convergence of  with a general sequence of potentials $\{ \pot{s}{\lmeas}{} \}_{s = 1} ^n$. We consider the following assumption on the sequence of potentials $\{ \pot{t}{\lmeas}{} \}_{t = 1} ^n$.

$\displaystyle \sup_{\lstate \in \R^{\dimx}} \tbwtrans{}{0}{\lstate}{\lmeas} / \pot{1}{\lmeas}{\lstate} < \infty$ and $\displaystyle\sup_{\lstate \in \R^{\dimx}} \int \pot{t}{\lmeas}{\lstate_{t}} \bwtrans{}{t}{\lstate}{\lstate_{t}} \rmd \lstate_{t} \big/ \pot{t+1}{\lmeas}{\lstate} < \infty$ for all $t \in [1:n-1]$.

The following exponential deviation inequality is standard and is a direct application of <cit.>. In particular, it implies a $\bigo(1 / \sqrt{N})$ bound on the mean squared error $\| \filter{0}{N}{h} - \filter{0}{\lmeas}{h} \|_2$.
Assume <ref>. There exist constants $c_{1,n}, c_{2,n} \in (0, \infty)$ such that, for all $N \in \nset$, $\varepsilon > 0$ and bounded function $h : \rset^{\dimx} \mapsto \rset$,
\pP \left[ \big| \filter{0}{N}{h} - \filter{0}{\ltrmeas}{h} \big| \geq \varepsilon \right] \leq c_{1,n} \exp(- c_{2,n} N \varepsilon^2 \big/ |h|^2 _\infty )
where $|h|_\infty \eqdef \sup_{x \in \rset^{\dimx}} |h(x)|$.

We also furnish our estimator with an explicit non-asymptotic bound on its bias. Define $\smash{\bfilter{0}{N}{} = \pE \big[ \filter{0}{N}{} ]}$ where $\filter{0}{N}{} = N^{-1} \sum_{i = 1}^N \delta_{\bparticle^i _0}$ is the particle approximation produced by <Ref> and the expectation is with respect to the law of $(\particle^{1:N}_{1:n}, \anc{1:N}{0:n-1})$.
Define for all $t \in [1:n]$, $\filter{t}{\star}{\lstate_t} \propto \bwmarg{t}{\lstate_t} \int \delta_\lmeas(\rmd \ltopstate_0) \bwtrans{}{0|t}{\lstate_t}{\lstate_0} \rmd \lbottomstate_0 \eqsp,$ where $\bwtrans{}{0|t}{\lstate_t}{\lstate_0} \eqdef \int \left\{ \prod_{s = 0} ^{t-1} \bwtrans{}{s}{\lstate_{s+1}}{\lstate_s} \right\} \rmd \lstate_{1:t-1}$.

It holds that

\begin{equation}
\kldivergence{\filter{0}\lmeas{}}{\bfilter{0}{N}{}} \leq \mathsf{C}^{\lmeas} _{0:n}(N-1)^{-1} + \mathsf{D}^{\lmeas} _{0:n} N^{-2} \eqsp,
\end{equation}

where  $\mathsf{D}^\lmeas _{0:n} > 0$,
$\mathsf{C}^\lmeas _{0:n} \eqdef
\textstyle\sum_{t = 1} ^n \int \frac{\normconst_t / \normconst_0}{\pot{t}{\lmeas}{z_t}} \left\{ \int \delta_{\lmeas}(\rmd \ltopstate_0)  \bwtrans{}{0|t}{z_t}{x_0} \rmd \lbottomstate_0 \right\} \filter{t}{\star}{\rmd z_t} $
and $\normconst_t \eqdef \int \pot{t}{\lmeas}{\lstate_t} \bwmarg{t}{\rmd \lstate_t}$ for all $t \in [1:n]$ and $\normconst_0 \eqdef \int \delta_\lmeas(\rmd \ltopstate_0) \bwmarg{0}{\lstate_0} \rmd \lbottomstate_0$. If furthermore <ref> holds then both $\mathsf{C}^{\lmeas} _{0:n}$ and $\mathsf{D}^{\lmeas} _{0:n}$ are finite.

The proof of <Ref> is postponed to <Ref>. <ref> is an assumption on the equivalent of the weights $\{\uweight{}{t}\}_{t = 0}^n$ with a general sequence of potentials $\{ \pot{t}{\lmeas}{}\}_{t = 1} ^n$ and is not restrictive as it can be satisfied by setting for example $\pot{s}{\lmeas}{\lstate_s} = \tfwtrans{s|0}{\lmeas}{\ltopstate_s} + \delta$ where $\delta > 0$. The resulting algorithm is then only a slight modification of the one described above, see <Ref> for more details.
It is also worth noting that <Ref> combined with Pinsker's inequality implies that the bias of  goes to $0$ with the number of particle samples $N$ for fixed $n$. We have chosen to present a bound in Kullback–Leibler (KL) divergence, inspired by [1, 21], as it allows an explicit dependence on the modeling choice $\{ \pot{s}{\lmeas}{} \}_{s = 1} ^n$, see <Ref>.
Finally, unlike the theoretical guarantees established for  in [43], proving the asymptotic exactness of our methodology w.r.t. to the generative model posterior does not require having $\bwmarg{s+1}{\lstate_{s+1}} \bwtrans{}{s}{\lstate_{s+1}}{\lstate_{s}} = \bwmarg{s}{\lstate_s} \fwtrans{s+1}{\lstate_s}{\lstate_{s+1}}$ for all $\smash{s \in [0:n-1]}$,
which does not in hold in practice.

Noisy case.

We consider the case $\noisestd > 0$. The posterior density is given by
\noisyfilter{0}{\lmeas}{\lstate_0} \propto \pot{0}{\lmeas}{\ltopstate_0} \bwmarg{0}{\lstate_0}$, where $\pot{0}{\lmeas}{\lstate_0} \eqdef \normpdf(\lmeas; \ltopstate_0, \noisestd^2 \Id_{\dimorig})$.
In what follows, assume that there exists $\tau \in [1:n]$ such that $\sigma^2 = (1 - \alphacumprod{\tau}{}) \big/ \alphacumprod{\tau}{}$. We denote $\tilde\lmeas_\tau = \alphacumprod{\tau}{}^{1/2} \lmeas$. We can then write that

\begin{equation}
\label{eq:pot_as_foward}
\pot{0}{\lmeas}{\ltopstate_0} = \alphacumprod{\tau}{}^{1/2} \cdot \normpdf(\tilde\lmeas_\tau; \alphacumprod{\tau}{}^{1/2} \lstate_0, (1 - \alphacumprod{\tau}{}) \cdot \Id_{\dimorig})
= \alphacumprod{\tau}{}^{1/2} \cdot \tfwtrans{\tau|0}{\ltopstate_0}{\tilde\lmeas_\tau} \eqsp,
\end{equation}

which hints that the likelihood function $\pot{0}{\lmeas}{}$ is closely related to the forward process (<ref>). We may then write the posterior $\noisyfilter{0}{\lmeas}{\lstate_0}$ as
$\filter{0}{\lmeas}{\lstate_0}  \propto \tfwtrans{\tau|0}{\ltopstate_0}{\tilde\lmeas_\tau} \bwmarg{0}{\lstate_0} \propto \int \delta_{\tilde\lmeas_\tau}(\rmd \ltopstate_\tau) \fwtrans{\tau|0}{\lstate_0}{\lstate_\tau} \bwmarg{0}{\lstate_0} \rmd \lbottomstate_\tau$.

Next, assume that the forward process (<ref>) is the reverse of the backward one (<ref>), i.e. that

\begin{equation}
\label{eq:bweqfw}
\bwmarg{t}{\lstate_t} \fwtrans{t+1}{\lstate_t}{\lstate_{t+1}} = \bwmarg{t+1}{\lstate_{t+1}} \bwtrans{}{t}{\lstate_{t+1}}{\lstate_t}\eqsp,  \quad \forall t \in [0:n-1] \eqsp.
\end{equation}

This is similar to the assumption made in   [43]. Then, it is easily seen that it implies $\bwmarg{0}{\lstate_0} \fwtrans{\tau|0}{\lstate_0}{\lstate_{\tau}} = \bwmarg{\tau}{\lstate_{\tau}} \bwtrans{}{0|\tau}{\lstate_{\tau}}{\lstate_0}$ and thus

\begin{equation}
\label{eq:noiseless_at_tau}
\filter{0}{\lmeas}{\lstate_0} = {\int \bwtrans{}{0|\tau}{\lstate_\tau}{\lstate_0} \delta_{\tilde\lmeas_\tau}(\rmd \ltopstate_\tau) \bwmarg{\tau}{\lstate_\tau} \rmd \lbottomstate_\tau}\big/{\int \delta_{\tilde\lmeas_\tau}(\rmd \overline{z}_\tau) \bwmarg{\tau}{z_\tau} \rmd \underline{z}_\tau}
= \int \bwtrans{}{0|\tau}{\cat{{\tilde\lmeas_\tau}}{\lbottomstate_\tau}}{\lstate_0} \filter{\tau}{\tilde\lmeas_\tau}{\rmd \lbottomstate_\tau} \eqsp,
\end{equation}

where $\filter{\tau}{\tilde\lmeas_\tau}{\lbottomstate_\tau} \propto \bwmarg{\tau}{\cat{\tilde\lmeas_\tau}{\lbottomstate_\tau}}$.
(<ref>) highlights that solving the inverse problem (<ref>) with $\noisestd > 0$ is equivalent to solving an inverse problem on the intermediate state $\state_{\tau} \sim \bwmarg{\tau}{}$ with noiseless observation $\tilde\lmeas_\tau$ of the $\dimorig$ top coordinates and then propagating the resulting posterior back to time $0$ with the backward kernel $\bwtrans{}{0|\tau}{}{}$.

The assumption (<ref>) does not always holds in realistic settings.
Therefore, while (<ref>) also holds only approximately in practice, we can still use it as inspiration for designing potentials when the assumption is not valid.
Consider then $\{ \pot{t}{\lmeas}{} \}_{t = \tau} ^n$ and sequence of probability measures $\{ \filter{t}{\lmeas}{} \}_{t = \tau} ^n$ defined for all $t \in [\tau: n]$ as
\filter{t}{\lmeas}{x_t} \propto \pot{t}{\lmeas}{\lstate_t} \bwmarg{t}{\lstate_t}$, where $\pot{t}{\lmeas}{\lstate_t }\eqdef\normpdf\left(\lstate_t; \tilde\lmeas_\tau, \big(1 - (1 - \inflationstd) \alphacumprod{t}{} / \alphacumprod{\tau}{} \big) \Id_{\dimorig}\right)$, $\inflationstd \geq 0$.
In the case of $\inflationstd = 0$, we have $\pot{t}{\lmeas}{\lstate_t} = \tfwtrans{t|\tau}{\tilde\lmeas_\tau}{\ltopstate_t}$ for $t \in [\tau + 1:n]$ and $\filter{\tau}{\lmeas}{} = \filter{\tau}{\tilde\lmeas_\tau}{}$. The recursion (<ref>) holds for $t \in [\tau:n]$ and assuming $\inflationstd > 0$, we find that
\filter{0}{\lmeas}{\lstate_0} \propto {\pot{0}{\lmeas}{\lstate_0}} \int {\pot{\tau}{\lmeas}{\lstate_\tau}}^{-1} \bwtrans{}{0|\tau}{\lstate_\tau}{\lstate_0} \filter{\tau}{\lmeas}{\rmd \lstate_\tau} \eqsp,
which resembles the recursion (<ref>). In practice we take $\kappa$ to be small in order to mimick the Dirac delta mass at $\ltopstate_\tau$ in (<ref>). Having a particle approximation $\filter{\tau}{N}{} = N^{-1} \sum_{i = 1}^N \delta_{\particle^i _\tau}$ of $\filter{\tau}{\lmeas}{}$ by adapting <Ref>, we estimate $\filter{0}{\lmeas}{}$ with $
\filter{0}{N}{} = \textstyle \sum_{i = 1}^N \weight{i}{0} \delta_{\particle^i _0}$ where $\particle^i _0 \sim \bwtrans{}{0|\tau}{\particle^i _\tau}{\cdot }$ and $\weight{i}{0} \propto \pot{0}{\lmeas}{\particle^i _0}$.
In the next section we extend this methodology to general linear Gaussian observation models.
Finally, note that (<ref>) allows us to extend   to handle noisy inverse problems in a principled manner. This is detailed in <Ref>.

§.§ Extension to general linear inverse problems

§ NUMERICS

The focus of this work is on providing an algorithm that consistently approximates the posterior distribution of a linear inverse problem with Gaussian noise.
A prerequisite for quantitative evaluation in ill-posed inverse problems in a Bayesian setting is to have access to samples of the posterior distribution.
This generally requires having at least an unnormalized proxy of the posterior density, so that one can run MCMC samplers such as the No U-turn sampler (NUTS) [20].
Therefore, this section focus on mixture models of two types of basis distribution, the Gaussian and the Funnel distributions.
We then present a brief illustration of  on image data. However, in this setting, the actual posterior distribution is unknown and the main goal is to explore the potentially multimodal posterior distribution, which makes a comparison with a "real image" meaningless.
Therefore, metrics such as Fréchet Inception Distance (FID) and LPIPS score, which require comparison to a ground truth, are not useful for evaluating Bayesian reconstruction methods in such settings.[The code for the experiments is available at <https://anonymous.4open.science/r/mcgdiff/README.md>.]

Mixture Models:

§.§.§ Image datasets

We now present samples from  in different image dataset and different kinds of inverse problems.

*Super Resolution

We start by super resolution.
We set $\sigma_y = 0.05$ for all the datasets
and $\zeta_{\mathrm{coeff}} = 0.1$ for  . We use $100$ steps of  with $\eta=1$. The results are shown in <Ref>.
We use a downsampling ratio of $4$ for the CIFAR-10 dataset, $8$ for both Flowers and Cats datasets and $16$ for the others.
The dimension of the datasets are recalled in <ref>.
We display in <ref> samples from , and over several different image datasets (<ref>).
For each algorithm, we
generate $1000$ samples and we show the pair of samples that are the furthest apart in $L^2$ norm from each other in the pool of samples.
For  we ran several parallel particle filters with $N=64$ to generate $1000$ samples.

CIFAR-10     Flowers     Cats     Bedroom     Church     CelebaHQ

Ratio 4 for CIFAR, 8 for flowers and Cats and 16 for CELEB

CIFAR-10     Flowers     Cats     Bedroom     Church     CelebaHQ
$(W, H, C)$     $(32, 32, 3)$     $(64, 64, 3)$     $(128, 128, 3)$     $(256, 256, 3)$     $(256, 256, 3)$     $(256, 256, 3)$

The datasets used for the inverse problems over image datasets.

*Gaussian 2D debluring
We consider a Gaussian 2D square kernel with sizes $(w / 6, h / 6)$ and standard deviation $w / 30$ where $(w, h)$ are the width and height of the image.
We set $\sigma_y = 0.1$ for all the datasets
and $\zeta_{\mathrm{coeff}} = 0.1$ for  . We use $100$ steps of  with $\eta=1$.
We display in <ref> samples from , and over several different image datasets (<ref>).
For each algorithm, we
generate $1000$ samples and we show the pair of samples that are the furthest appart in $L^2$ norm from each other in the pool of samples.
For  we ran several parallel particle filters with $N=64$ to generate $1000$ samples.

CIFAR 10     Flowers     Cats     Bedroom     Church     CelebaHQ

Inpainting on CelebA

Inpainting with different masks on the CelebA test set.

We consider the inpainting problem on the CelebA dataset with several different masks in <ref>.
We show in <ref> the evolution of the particle cloud with $s$.

$1000$     $900$     $800$     $700$     $600$     $500$

$400$     $300$     $200$     $100$     $50$     $4$

Evolution of the particle cloud for one of the masks. The numbers on top and bottom indicate the step $s$ of the approximation.

§ CONCLUSION

In this paper, we present   a novel method for solving Bayesian linear Gaussian inverse problems with SGM priors.
We show that  is theoretically grounded and provided numerical experiments that reflect the adequacy of  in a Bayesian framework, as opposed to recent works.
This difference is of the uttermost importance when the relevance of the generated samples is hard to verify, as in safety critical applications.
is a first step towards robust approaches for addressing the challenges of Bayesian linear inverse problems with SGM priors.

Acknowledgments GC, YJEL and EM would like to thank the Isaac Newton Institute for Mathematical Sciences for support and hospitality during the program The mathematical and statistical foundation of future data-driven engineering when work on this paper was undertaken

[1]
Andrieu, C., Lee, A., and Vihola, M. (2018).
Uniform ergodicity of the iterated conditional SMC and geometric
ergodicity of particle Gibbs samplers.
Bernoulli, 24(2):842 – 872.

[2]
Arjomand Bigdeli, S., Zwicker, M., Favaro, P., and Jin, M. (2017).
Deep mean-shift priors for image restoration.
Advances in Neural Information Processing Systems, 30.

[3]
Benton, J., Shi, Y., De Bortoli, V., Deligiannidis, G., and Doucet, A. (2022).
From denoising diffusions to denoising markov models.
arXiv preprint arXiv:2211.03595.

[4]
Besag, J., York, J., and Mollié, A. (1991).
Bayesian image restoration, with two applications in spatial
Annals of the institute of statistical mathematics, 43:1–20.

[5]
Cappé, O., Moulines, E., and Ryden, T. (2005).
Inference in Hidden Markov Models (Springer Series in
Springer-Verlag, Berlin, Heidelberg.

[6]
Cheney, M. and Borden, B. (2009).
Fundamentals of radar imaging.

[7]
Chopin, N., Papaspiliopoulos, O., et al. (2020).
An introduction to sequential Monte Carlo, volume 4.

[8]
Chung, H., Kim, J., Mccann, M. T., Klasky, M. L., and Ye, J. C. (2023).
Diffusion posterior sampling for general noisy inverse problems.
In The Eleventh International Conference on Learning

[9]
Dai, C., Heng, J., Jacob, P. E., and Whiteley, N. (2022).
An invitation to sequential monte carlo samplers.
Journal of the American Statistical Association,

[10]
Dashti, M. and Stuart, A. M. (2017).
The bayesian approach to inverse problems.
In Handbook of uncertainty quantification, pages 311–428.

[11]
Dinh, L., Sohl-Dickstein, J., and Bengio, S. (2017).
Density estimation using real NVP.
In 5th International Conference on Learning Representations,
ICLR 2017, Toulon, France, April 24-26, 2017, Conference Track
Proceedings. OpenReview.net.

[12]
Douc, R., Moulines, E., and Stoffer, D. (2014).
Nonlinear time series: Theory, methods and applications with R
CRC press.

[13]
Doucet, A., De Freitas, N., Gordon, N. J., et al. (2001).
Sequential Monte Carlo methods in practice, volume 1.

[14]
Elbakri, I. A. and Fessler, J. A. (2002).
Statistical image reconstruction for polyenergetic x-ray computed
IEEE transactions on medical imaging, 21(2):89–99.

[15]
Fergus, R., Singh, B., Hertzmann, A., Roweis, S. T., and Freeman, W. T. (2006).
Removing camera shake from a single photograph.
In Acm Siggraph 2006 Papers, pages 787–794.

[16]
Figueiredo, M. A., Bioucas-Dias, J. M., and Nowak, R. D. (2007).
Majorization–minimization algorithms for wavelet-based image
IEEE Transactions on Image processing, 16(12):2980–2991.

[17]
González, R. C., Woods, R. E., and Masters, B. R. (2009).
Digital image processing, third edition.
Journal of Biomedical Optics, 14:029901.

[18]
Gui, J., Sun, Z., Wen, Y., Tao, D., and Ye, J. (2021).
A review on generative adversarial networks: Algorithms, theory, and
IEEE transactions on knowledge and data engineering.

[19]
Ho, J., Jain, A., and Abbeel, P. (2020).
Denoising diffusion probabilistic models.
Advances in Neural Information Processing Systems,

[20]
Hoffman, M. and Gelman, A. (2011).
The no-u-turn sampler: Adaptively setting path lengths in hamiltonian
monte carlo.
Journal of Machine Learning Research, 15.

[21]
Huggins, J. H. and Roy, D. M. (2019).
Sequential Monte Carlo as approximate sampling: bounds, adaptive
resampling via $\infty$-ESS, and an application to particle Gibbs.
Bernoulli, 25(1):584 – 622.

[22]
Idier, J. (2013).
Bayesian approach to inverse problems.
John Wiley & Sons.

[23]
Ivanov, O., Figurnov, M., and Vetrov, D. (2018).
Variational autoencoder with arbitrary conditioning.
arXiv preprint arXiv:1806.02382.

[24]
Kaipio, J. P., Kolehmainen, V., Somersalo, E., and Vauhkonen, M. (2000).
Statistical inversion and monte carlo sampling methods in electrical
impedance tomography.
Inverse problems, 16(5):1487.

[25]
Kaltenbach, S., Perdikaris, P., and Koutsourelakis, P.-S. (2023).
Semi-supervised invertible neural operators for bayesian inverse
Computational Mechanics, pages 1–20.

[26]
Kawar, B., Elad, M., Ermon, S., and Song, J.
Denoising diffusion restoration models.
In Advances in Neural Information Processing Systems.

[27]
Kingma, D. P. and Ba, J. (2015).
Adam: A method for stochastic optimization.
In Bengio, Y. and LeCun, Y., editors, 3rd International
Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May
7-9, 2015, Conference Track Proceedings.

[28]
Kingma, D. P., Welling, M., et al. (2019).
An introduction to variational autoencoders.
Foundations and Trends® in Machine Learning,

[29]
Kobyzev, I., Prince, S. J., and Brubaker, M. A. (2020).
Normalizing flows: An introduction and review of current methods.
IEEE transactions on pattern analysis and machine intelligence,

[30]
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, pages 11461–11471.

[31]
Marnissi, Y., Zheng, Y., Chouzenoux, E., and Pesquet, J.-C. (2017).
A variational bayesian approach for image restoration—application
to image deblurring with poisson–gaussian noise.
IEEE Transactions on Computational Imaging, 3(4):722–737.

[32]
Peng, J., Liu, D., Xu, S., and Li, H. (2021).
Generating diverse structure for image inpainting with hierarchical
In Proceedings of the IEEE/CVF Conference on Computer Vision and
Pattern Recognition, pages 10775–10784.

[33]
Pitt, M. K. and Shephard, N. (1999).
Filtering via simulation: Auxiliary particle filters.
Journal of the American statistical association,

[34]
Sahlström, T. and Tarvainen, T. (2023).
Utilizing variational autoencoders in the bayesian inverse problem of
photoacoustic tomography.
SIAM Journal on Imaging Sciences, 16(1):89–110.

[35]
Shin, H. and Choi, M. (2023).
Physics-informed variational inference for uncertainty quantification
of stochastic differential equations.
Journal of Computational Physics, page 112183.

[36]
Sohl-Dickstein, J., Weiss, E., Maheswaranathan, N., and Ganguli, S. (2015).
Deep unsupervised learning using nonequilibrium thermodynamics.
In International Conference on Machine Learning, pages
2256–2265. PMLR.

[37]
Song, J., Meng, C., and Ermon, S. (2021a).
Denoising diffusion implicit models.
In International Conference on Learning Representations.

[38]
Song, Y., Durkan, C., Murray, I., and Ermon, S. (2021b).
Maximum likelihood training of score-based diffusion models.
Advances in Neural Information Processing Systems,

[39]
Song, Y., Shen, L., Xing, L., and Ermon, S. (2022).
Solving inverse problems in medical imaging with score-based
generative models.
In International Conference on Learning Representations.

[40]
Song, Y., Sohl-Dickstein, J., Kingma, D. P., Kumar, A., Ermon, S., and Poole,
B. (2021c).
Score-based generative modeling through stochastic differential
In International Conference on Learning Representations.

[41]
Stuart, A. M. (2010).
Inverse problems: a bayesian perspective.
Acta numerica, 19:451–559.

[42]
Su, J., Xu, B., and Yin, H. (2022).
A survey of deep learning approaches to image restoration.
Neurocomputing, 487:46–65.

[43]
Trippe, B. L., Yim, J., Tischer, D., Baker, D., Broderick, T., Barzilay, R.,
and Jaakkola, T. S. (2023).
Diffusion probabilistic modeling of protein backbones in 3d for the
motif-scaffolding problem.
In The Eleventh International Conference on Learning

[44]
Vlaardingerbroek, M. T. and Boer, J. A. (2013).
Magnetic resonance imaging: theory and practice.
Springer Science & Business Media.

[45]
Wan, Z., Zhang, J., Chen, D., and Liao, J. (2021).
High-fidelity pluralistic image completion with transformers.
In Proceedings of the IEEE/CVF International Conference on
Computer Vision, pages 4692–4701.

[46]
Wei, X., van Gorp, H., Gonzalez-Carabarin, L., Freedman, D., Eldar, Y. C., and
van Sloun, R. J. (2022).
Deep unfolding with normalizing flow priors for inverse problems.
IEEE Transactions on Signal Processing, 70:2962–2971.

[47]
Wu, L., Trippe, B. L., Naesseth, C. A., Blei, D. M., and Cunningham, J. P.
Practical and asymptotically exact conditional sampling in diffusion
arXiv preprint arXiv:2306.17775.

[48]
Xiang, H., Zou, Q., Nawaz, M. A., Huang, X., Zhang, F., and Yu, H. (2023).
Deep learning for image inpainting: A survey.
Pattern Recognition, 134:109046.

[49]
Yeh, R. A., Lim, T. Y., Chen, C., Schwing, A. G., Hasegawa-Johnson, M., and Do,
M. N. (2018).
Image restoration with deep generative models.
In 2018 IEEE International Conference on Acoustics, Speech and
Signal Processing (ICASSP), pages 6772–6776. IEEE.

[50]
Yu, J., Lin, Z., Yang, J., Shen, X., Lu, X., and Huang, T. S. (2018).
Generative image inpainting with contextual attention.
In Proceedings of the IEEE conference on computer vision and
pattern recognition, pages 5505–5514.

[51]
Zeng, Y., Fu, J., Chao, H., and Guo, B. (2022).
Aggregated contextual transformations for high-resolution image
IEEE Transactions on Visualization and Computer Graphics.

[52]
Zhang, G., Ji, J., Zhang, Y., Yu, M., Jaakkola, T., and Chang, S. (2023).
Towards coherent image inpainting using denoising diffusion implicit
arXiv preprint arXiv:2304.03322.

[53]
Zheng, C., Cham, T.-J., and Cai, J. (2019).
Pluralistic image completion.
In Proceedings of the IEEE/CVF Conference on Computer Vision and
Pattern Recognition, pages 1438–1447.

[54]
Zhihang, X., Yingzhi, X., and Qifeng, L. (2023).
A domain-decomposed vae method for bayesian inverse problems.
arXiv preprint arXiv:2301.05708.

§  EXTENSION

§ PROOFS

§.§ Proof of <Ref>

§.§ Preliminary definitions.

We preface the proof with notations and definitions of a few quantities that will be used throughout.

For a probability measure $\mu$ and $f$ a bounded measurable function, we write $\mu(f) \eqdef \int f(x) \mu(\rmd x)$  the expectation of $f$ under $\mu$ and if $K(\rmd x|z)$ is a transition kernel we write $K(f)(z) \eqdef \int f(x) K(\rmd x|z)$.

Define the smoothing distribution

\begin{equation}
\filter{0:n}{y}{\rmd \lstate_{0:n}} \propto \delta_{\lmeas}(\rmd \ltopstate_0) \bwmarg{0:n}{\lstate_{0:n}}  \rmd \lbottomstate_0 \rmd \lstate_{1:n} \eqsp,
\end{equation}

which admits the posterior $\filter{0}{\lmeas}{}$ as time $0$ marginal. Its particle estimate known as the poor man smoother is given by

\begin{equation}
\filter{0:n}{N}{\rmd x_{0:n}} = N^{-1} \sum_{k_{0:n} \in [1:N]^{n+1}}  \delta_{\cat{\lmeas}{\bparticle^{k_0}_{0}}}(\rmd x_0) \prod_{s = 1}^n \1 \big\{ k_s = \anc{k_{s-1}}{s} \big\} \delta_{\particle^{k_s} _{s}}(\rmd x_s) \eqsp.
\end{equation}

We also let $\bfilter{0:n}{N}{}$ be the probability measure defined for any $B \in \sigalgebra(\R^{\dimx})^{\otimes n+1}$ by

\[
\bfilter{0:n}{N}{B} = \pE \big[ \filter{0:n}{N}{B} \big] \eqsp,
\]

where the expectation is with respect to the probability measure

\begin{multline}
\lawSMC{0:n}\big(\rmd (\lstate^{1:N} _{0:n}, \ianc{1:N}{1:n}) \big) = \prod_{i = 1}^N \bwopt{n}{\lmeas}{}{}(\rmd \lstate^i _n) \prod_{\ell = 2}^n \left\{ \prod_{j = 1}^{N} \sum_{k = 1}^N \weight{k}{\ell-1} \delta_k(\rmd \ianc{j}{\ell}) \bwopt{\ell-1}{y}{\iparticle{\ianc{j}{\ell}}{\ell}}{\rmd \iparticle{j}{\ell-1}} \right\} \\
\times \prod_{j = 1}^{N} \sum_{k = 1}^N \weight{k}{0} \delta_k(\rmd \ianc{j}{1}) \bbwopt{0}{y}{\iparticle{\ianc{j}{1}}{1}}{\rmd \lbottomstate^j _0} \delta_{\lmeas}(\rmd \ltopstate^j _0) \eqsp,
\end{multline}

where $\weight{i}{t} \eqdef \uweight{}{t}(\particle^i _{t+1}) / \sum_{j = 1}^N \uweight{}{t}(\particle^j _{t+1})$ and
which corresponds to the joint law of all the random variables generated by <Ref>.
It then follows by definition that for any $C \in \sigalgebra(\R^{\dimx})$,

\[
\int \bfilter{0:n}{N}{\rmd z_{0:n}} \1_C(z_0) = \pE \left[ \int \filter{0:n}{N}{\rmd z_{0:n}} \1_C(z_0) \right] = \pE \big[ \filter{0}{N}{C} \big] = \bfilter{0}{N}{C} \eqsp.
\]

Define also the law of the conditional particle cloud

\begin{equation}
\begin{aligned}
\label{eq:cSMCkernel}
& \lawcSMC{}\big(\rmd (\iparticle{1:N}{0:n}, \ianc{1:N}{1:n})\big| z_{0:n} \big) =  \delta_{z_n}(\rmd \iparticle{N}{n}) \prod_{i = 1}^{N-1} \bwopt{n}{y}{}{}(\rmd \iparticle{i}{n}) \\
& \hspace{2cm} \times\prod_{\ell = 2}^n \delta_{z_{\ell -1}}(\rmd \iparticle{N}{\ell-1}) \delta_N(\rmd \ianc{N}{\ell-1}) \prod_{j = 1}^{N-1} \sum_{k = 1}^N \weight{k}{\ell-1} \delta_k(\rmd \ianc{j}{\ell}) \bwopt{\ell-1}{y}{\iparticle{\ianc{j}{\ell}}{\ell}}{\rmd \iparticle{j}{\ell-1}} \\
& \hspace{2cm} \times \delta_{z_0}(\rmd \lstate^{N} _{0}) \delta_N(\rmd \ianc{N}{1}) \prod_{j = 1}^{N-1} \sum_{k = 1}^N \weight{k}{0} \delta_k(\rmd \ianc{j}{1}) \bbwopt{0}{y}{\iparticle{\ianc{j}{1}}{1}}{\rmd \lbottomstate^j _0} \delta_\lmeas(\rmd \ltopstate^j _0)\eqsp.
\end{aligned}
\end{equation}

In what follows $\pE _{z_{0:n}}$ refers to expectation with respect to $\lawcSMC{}(\cdot | z_{0:n})$.
Finally, for $s \in [0:n-1]$ we let $\sumweight{s}$ denote the sum of the filtering weights at step $s$, i.e. $\sumweight{s} = \sum_{i = 1}^N \uweight{}{s}(\particle^{i}_{s+1})$. We also write $\normconst_0 = \int \bwmarg{0}{\lstate_0} \delta_\lmeas(\rmd \ltopstate_0) \rmd \lbottomstate_0$ and for all $\ell \in [1:n]$, $\normconst_\ell = \int \tfwtrans{\ell|0}{\lmeas}{\ltopstate_\ell} \bwmarg{\ell}{\rmd \lstate_\ell}.$

The proof of <Ref> relies on two Lemmata stated below and proved in <Ref>; in <Ref> we provide an expression for the Radon-Nikodym derivative $\rmd \filter{0:n}{y}{} / \rmd \bfilter{0:n}{y}{}$ and in <Ref> we explicit its leading term.
$\filter{0:n}{y}{}$ and $\bfilter{0:n}{N}{}$ are equivalent and we have that

\begin{equation}
\bfilter{0:n}{N}{\rmd z_{0:n}} = \pE _{z_{0:n}} \left[ \frac{N^{n} \normconst_0 / \normconst_n}{\prod_{s = 0}^{n-1} \sumweight{s}} \right] \filter{0:n}{y}{\rmd z_{0:n}} \eqsp.
\end{equation}

It holds that

\begin{multline}
\label{eq:normconst_closedform}
\frac{\normconst_n}{\normconst_0} \pE_{z_{0:n}} \left[ \prod_{s = 0} ^{n-1} N^{-1} \sumweight{s} \right]
= \left( \frac{N-1}{N} \right)^{n}  \\
+ \frac{(N-1)^{n-1}}{N^n} \sum_{s = 1} ^n \frac{\normconst_s / \normconst_0}{\tfwtrans{s|0}{\lmeas}{\overline{z}_s} } \int \bwtrans{}{0|s}{z_s}{x_0} \delta_{\lmeas}(\rmd \ltopstate_0) \rmd \lbottomstate_0 + \frac{\mathsf{D}^\lmeas _{0:n}}{N^2} \eqsp.
\end{multline}

where $\mathsf{D}^\lmeas _{0:n}$ is a positive constant.

Before proceeding with the proof of <Ref>, let us note that having $z \mapsto \uweight{}{\ell}(z)$ bounded on $\R^{\dimx}$ for all $\ell \in [0:n-1]$ is sufficient to guarantee that $\mathsf{C}^\lmeas _{0:n}$ and $\mathsf{D}^\lmeas _{0:n}$ are finite since in this case it follows immediately that $\pE_{z_{0:n}} \left[ \prod_{s = 0} ^{n-1} N^{-1} \sumweight{s} \right]$ is bounded and so is the right hand side of (<ref>). This can be achieved with a slight modification of (<ref>). Indeed, consider instead the following recursion for $s \in [0:n]$ where $\delta > 0$,

\begin{align*}
\filter{n}{\lmeas}{\lstate_n} & \propto \big( \tfwtrans{n|0}{\lmeas}{\ltopstate_n} + \delta \big)  \bwmarg{n}{\lstate_n} \eqsp, \\
\filter{s}{\lmeas}{\lstate_s} & \propto \int \filter{s+1}{\lmeas}{\lstate_{s+1}} \bwtrans{}{s}{\lstate_{s+1}}{\rmd \lstate_s} \frac{\tfwtrans{s}{\lmeas}{\ltopstate_s} + \delta}{\tfwtrans{s+1}{\lmeas}{\ltopstate_{s+1}} + \delta} \rmd \lstate_{s+1} \eqsp.
\end{align*}

Then we have that

\begin{equation*}
\filter{0}{\lmeas}{\lbottomstate_0} \propto \int \filter{1}{\lmeas}{\lstate_1} \bbwtrans{}{0}{\lstate_1}{\lbottomstate_0} \frac{\tbwtrans{}{0}{\lstate_1}{\lmeas}}{\tfwtrans{1|0}{\lmeas}{\ltopstate_1} + \delta} \rmd \lstate_1 \eqsp.
\end{equation*}

We can then use <Ref> to produce a particle approximation of  $\filter{0}{\lmeas}{}$ using the following transition and weight function,

\begin{align*}
\bwopt{s}{\lmeas, \delta}{\lstate_{s+1}}{\lstate_s} & = \frac{\gamma_s(\lmeas | \lstate_{s+1})}{\gamma_s(\lmeas | \lstate_{s+1}) + \delta} \bwopt{s}{\lmeas}{\lstate_{s+1}}{\lstate_s} + \frac{\delta}{\gamma_s(\lmeas | \lstate_{s+1}) + \delta} \bwtrans{}{s}{\lstate_{s+1}}{\lstate_s} \eqsp,  \\
\uweight{}{s}(\lstate_{s+1}) & = \big( \gamma_s(\lmeas | \lstate_{s+1}) + \delta \big) \big/ \big( \tfwtrans{s+1|0}{\lmeas}{\ltopstate_{s+1}} + \delta \big) \eqsp,
\end{align*}

where $\gamma_s(\lmeas | \lstate_{s+1}) = \int \tfwtrans{s|0}{\lmeas}{\ltopstate_s} \bwtrans{}{s}{\lstate_{s+1}}{\lstate_s} \rmd \lstate_s$ is available in closed form and $\bwopt{s}{\lmeas}{}{}$ is defined in (<ref>). $\uweight{}{s}$ is thus clearly bounded for all $s \in [0:n-1]$ and it is still possible to sample from $\bwopt{s}{\lmeas, \delta}{}{}$ since it is simply a mixture between the transition (<ref>) and the “prior” transition.

Consider the forward Markov kernel

\begin{equation}
\fwdker{1:n}{z_0, }{\rmd z_{1:n}} = \frac{\bwmarg{1:n}{\rmd z_{1:n}} \bwtrans{}{0}{z_1}{z_0}}{\int \bwmarg{1:n}{\rmd \tilde{z}_{1:n}} \bwtrans{}{0}{\tilde{z}_1}{\tilde{z}_0}} \eqsp,
\end{equation}

which satisfies

\[
\filter{0:n}{y}{\rmd z_{0:n}} = \filter{0}{y}{\rmd z_0} \fwdker{1:n}{z_0,}{\rmd z_{1:n}} \eqsp.
\]

By <Ref> we have for any $C \in \sigalgebra(\R^{\dimx})$ that

\begin{align*}
\bfilter{0}{N}{C} & = \int \bfilter{0:n}{N}{\rmd z_{0:n}} \1_C(z_0) \\
& = \int \1_C(z_0) \pE _{z_{0:n}} \left[ \frac{N^n \normconst_0 / \normconst_n}{\prod_{s = 0}^{n-1} \sumweight{s}} \right] \filter{0:n}{y}{\rmd z_{0:n}} \\
& = \int \1_C(z_0) \int \fwdker{1:n}{z_0,}{\rmd z_{1:n}} \pE _{z_{0:n}} \left[ \frac{N^n \normconst_0 / \normconst_n}{\prod_{s = 0}^{n-1} \sumweight{s}} \right] \filter{0}{y}{\rmd z_{0}} \eqsp,
\end{align*}

which shows that the Radon-Nikodym derivative $\rmd \bfilter{0}{N}{} / \rmd \filter{0}{y}{}$ is,

\[
\frac{\rmd \bfilter{0}{N}{}}{\rmd \filter{0}{y}{}}(z_0) = \int \fwdker{1:n}{z_0,}{\rmd z_{1:n}} \pE _{z_{0:n}} \left[ \frac{N^n \normconst_0 / \normconst_n}{\prod_{s = 0}^{n-1} \sumweight{s}} \right] \eqsp.
\]

Applying Jensen's inequality twice yields

\[
\frac{\rmd \bfilter{0}{N}{}}{\rmd \filter{0}{y}{}}(z_{0}) \geq \frac{N^n \normconst_0 / \normconst_n}{\int \fwdker{1:n}{z_0,}{\rmd z_{1:n}} \pE _{z_{0:n}} \left[ \prod_{s = 0}^{n-1} \sumweight{s} \right]} \eqsp,
\]

and it then follows that

\[
\kldivergence{\filter{0}{y}{}}{\bfilter{0}{N}{}} \leq \int \log \left( \frac{\normconst_n}{\normconst_0}  \int \fwdker{1:n}{z_0,}{\rmd z_{1:n}} \pE _{z_{0:n}} \left[ \prod_{s = 0}^{n-1} N^{-1} \sumweight{s} \right] \right) \filter{0}{y}{\rmd z_0} \eqsp.
\]

Finally, using <Ref> and the fact that $\log(1 + x) < x$ for $x > 0$ we get

\begin{equation*}
\kldivergence{\filter{0}{y}{}}{\bfilter{0}{N}{}} \leq \frac{\mathsf{C}^\lmeas _{0:n}}{N-1} + \frac{\mathsf{D}^\lmeas _{0:n}}{N^2}
\end{equation*}

\[
\mathsf{C}^\lmeas _{0:n} \eqdef \sum_{s = 1} ^n \int \frac{\normconst_s / \normconst_0}{\tfwtrans{s|0}{\lmeas}{\overline{z}_s} } \left(  \bwtrans{}{0|s}{z_s}{x_0} \delta_{\lmeas}(\rmd \ltopstate_0) \rmd \lbottomstate_0 \right) \filter{s}{\lmeas}{\rmd z_s} \eqsp,
\]

and $\filter{s}{\lmeas}{z_s} \propto \bwmarg{s}{z_s} \int \bwtrans{}{0|s}{z_s}{z_0} \delta_{\lmeas}(\rmd \underline{z}_0) \rmd \overline{z}_0.$

§.§ Proof of <Ref> and <Ref>

We have that

\begin{align*}
& \bfilter{0:n}{N}{\rmd z_{0:n}} \\
& = N^{-1} \int \lawSMC{0:n}(\rmd \iparticle{1:N}{0:n}, \rmd \ianc{1:N}{1:n}) \sum_{k_{0:n} \in [1:N]^{n+1}}  \delta_{\cat{\lmeas}{\lbottomstate^{k_0}_{0}}}(\rmd z_0) \prod_{s = 1}^n \1 \big\{ k_s = \ianc{k_{s-1}}{s} \big\} \delta_{\lstate^{k_s} _{s}}(\rmd z_s) \\
& = N^{-1} \int \sum_{k_{0:n}} \sum_{\ianc{1:N}{1:n}}  \delta_{\cat{\lmeas}{\lbottomstate^{k_0}_{0}}}(\rmd z_0) \prod_{s = 1}^n \1 \big\{k_s = \ianc{k_{s-1}}{s}\big\} \delta_{\iparticle{k_s}{s}}(\rmd z_s) \\
& \hspace{.5cm} \times \prod_{j = 1}^N \bwopt{n}{y}{}{}(\rmd \iparticle{j}{n}) \left\{ \prod_{\ell = 2}^n  \prod_{i = 1}^N \weight{\ianc{i}{\ell}}{\ell-1} \bwopt{\ell-1}{y}{\iparticle{\ianc{i}{\ell}}{\ell}}{\rmd \iparticle{i}{\ell-1}} \right\} \prod_{r = 1}^N  \weight{\ianc{r}{1}}{0} \bbwopt{\ell-1}{y}{\iparticle{\ianc{r}{1}}{1}}{\rmd \lbottomstate^{r} _{0}} \delta_\lmeas(\ltopstate^r _0)\\
& = N^{-1}\int \sum_{k_{0:n}} \sum_{\ianc{1:N}{1:n}} \bwopt{n}{y}{}{}(\rmd \iparticle{k_n}{n}) \delta_{\iparticle{k_n}{n}}(\rmd z_n) \prod_{j \neq k_n} \bwopt{n}{y}{}{}(\rmd \iparticle{j}{n}) \prod_{\ell = 2}^n \bigg\{ \prod_{i \neq k_{\ell-1} } \weight{\ianc{i}{\ell}}{\ell-1} \bwopt{\ell-1}{y}{\iparticle{\ianc{i}{\ell}}{\ell}}{\rmd \iparticle{i}{\ell-1}}\\
& \hspace{.5cm} \times \1 \big\{ \ianc{k_{\ell-1}}{\ell} = k_\ell\} \wratio{}{\ell-1}{\iparticle{\ianc{k_{\ell-1}}{\ell}}{\ell}} \bwopt{\ell-1}{y}{\iparticle{\ianc{k_{\ell-1}}{\ell}}{\ell}}{\rmd \iparticle{k_{\ell-1}}{\ell}} \delta_{\iparticle{k_{\ell-1}}{\ell-1}}(\rmd z_{\ell-1}) \bigg\} \\
& \hspace{.5cm} \times \bigg\{ \prod_{r \neq k_0} \weight{\ianc{r}{1}}{0} \bbwopt{0}{y}{\iparticle{\ianc{r}{1}}{1}}{\rmd \lbottomstate^{r} _{0}} \delta_\lmeas(\rmd \ltopstate^r _0) \bigg\} \1 \big\{ \ianc{k_0}{1} = k_1 \big\} \wratio{}{0}{\iparticle{\ianc{k_{0}}{1}}{1}} \bwopt{0}{y}{\iparticle{\ianc{k_{0}}{1}}{0}}{\rmd \lbottomstate^{k_0} _0} \delta_{\cat{\lmeas}{\lbottomstate^{k_0} _0}}(\rmd z_{0}) \eqsp.
\end{align*}

Then, using that for all $s \in [2:n]$

\[
\uweight{}{s-1}(\iparticle{k_s}{s}) \bwopt{s-1}{y}{\iparticle{k_s}{s}}{\rmd \iparticle{k_{s-1}}{s-1}} = \frac{\tfwtrans{s-1|0}{\lmeas}{\ltopstate^{k_{s-1}} _{s-1}}}{\tfwtrans{s|0}{\lmeas}{\ltopstate^{k_s} _{s}}} \bwtrans{}{s}{\iparticle{k_s}{s}}{\rmd \iparticle{k_{s-1}}{s-1}} \eqsp,
\]

we recursively get that

\begin{align*}
& \bwopt{n}{y}{}{}(\rmd \iparticle{k_n}{n}) \delta_{\iparticle{k_n}{n}}(\rmd z_n) \prod_{s = 2}^{n}\1 \big\{ \ianc{k_{s-1}}{s} = k_s\} \wratio{}{s-1}{\iparticle{\ianc{k_{s-1}}{s}}{s}} \bwopt{s-1}{y}{\iparticle{\ianc{k_{s-1}}{s}}{s}}{\rmd \iparticle{k_{s-1}}{s-1}} \delta_{\iparticle{k_{s-1}}{s-1}}(\rmd z_{s-1})\\
& \hspace{2cm}  \times \1 \big\{ \ianc{k_0}{1} = k_1 \big\} \wratio{}{0}{\iparticle{\ianc{k_{0}}{1}}{1}} \bwopt{0}{y}{\iparticle{\ianc{k_{0}}{1}}{1}}{\rmd \lbottomstate^{k_0} _0} \delta_{\cat{\lmeas}{\lbottomstate^{k_0} _0}}(\rmd z_{0}) \\
& = \frac{\tfwtrans{n|0}{\lmeas}{z_n} \bwmarg{n}{\rmd z_n}}{\normconst_n} \delta_{z_n}(\rmd \iparticle{k_n}{n}) \prod_{s=2}^{n} \1 \big\{ \ianc{k_{s-1}}{s} = k_s\} \frac{\tfwtrans{s-1|0}{\lmeas}{\overline{z}_{s-1}}}{\sumweight{s-1} \tfwtrans{s|0}{\lmeas}{\overline{z}_s}} \bwtrans{}{s-1}{z_s}{\rmd z_{s-1}} \delta_{z_{s-1}}(\rmd \iparticle{k_{s-1}}{s-1}) \\
& \hspace{2cm} \times \1 \big\{ \ianc{k_0}{1} = k_1 \big\} \frac{\tbwtrans{}{0}{z_1}{\lmeas}}{\sumweight{0} \tfwtrans{1|0}{\lmeas}{\overline{z}_1}} \bwtrans{}{0}{z_1}{\rmd \underline{z}_0} \delta_\lmeas (\rmd \overline{z}_0) \delta_{z_0}(\rmd \lstate^{k_0} _{0}) \\
& = \frac{\normconst_0}{\normconst_n} \filter{0:n}{y}{\rmd z_{0:n}}  \delta_{z_n}(\rmd \iparticle{k_n}{n}) \prod_{s=1}^{n} \1 \big\{ \ianc{k_{s-1}}{s} = k_s\} \frac{1}{\sumweight{s-1}} \delta_{z_{s-1}}(\rmd \iparticle{k_{s-1}}{s-1}) \eqsp.
\end{align*}

Thus, we obtain

\begin{align*}
\bfilter{0:n}{N}{\rmd z_{0:n}} & = N^{-1} \int \sum_{k_{0:n}} \sum_{\ianc{1:N}{1:n}} \filter{0:n}{y}{\rmd z_{0:n}} \frac{\normconst_0 / \normconst_n}{\prod_{s = 0}^{n-1} \sumweight{s}} \delta_{z_n}(\rmd \iparticle{k_n}{n}) \prod_{j \neq k_n} \bwopt{n}{y}{}{}(\rmd \iparticle{j}{n})  \\
& \hspace{.5cm} \times \prod_{\ell = 2}^n \1\big\{ \ianc{k_{\ell-1}}{\ell} = k_\ell \big\} \delta_{z_{\ell -1}}(\rmd \iparticle{k_{\ell-1}}{\ell-1}) \prod_{i \neq k_{\ell-1}} \weight{\ianc{i}{\ell}}{\ell-1} \bwopt{\ell-1}{y}{\iparticle{\ianc{i}{\ell}}{\ell}}{\rmd \iparticle{i}{\ell-1}}  \\
& \hspace{.5cm} \times \1 \big\{ \ianc{k_0}{1} = k_1 \big\} \delta_{z_0}(\rmd \lstate^{k_0} _0) \prod_{i \neq k_0} \weight{\ianc{i}{1}}{0} \bbwtrans{}{0}{\lstate^{\ianc{i}{1}} _1}{\lbottomstate^i _0} \delta_\lmeas(\rmd \ltopstate^i _0) \\
& = N^{-1} \sum_{k_{0:n}} \filter{0:n}{y}{\rmd z_{0:n}} \pE^{k_{0:n}} _{z_{0:n}} \left[ \frac{\normconst_0 / \normconst_n}{\prod_{s = 0}^{n-1} \sumweight{s}} \right] \eqsp,
\end{align*}

where for all $k_{0:n} \in [1:N]^{n+1}$ $\pE^{k_{0:n}} _{z_{0:n}}$ denotes the expectation under the Markov kernel

\begin{align*}
& \lawcSMC{k_{0:n}}\big(\rmd (\iparticle{1:N}{0:n}, \ianc{1:N}{1:n})\big| z_{0:n} \big) =  \delta_{z_n}(\rmd \iparticle{k_n}{n}) \prod_{i \neq k_n} \bwopt{n}{y}{}{}(\rmd \iparticle{i}{n}) \\
& \hspace{2cm} \times\prod_{\ell = 2}^n \delta_{z_{\ell -1}}(\rmd \iparticle{k_{\ell-1}}{\ell-1}) \delta_{k_\ell}(\rmd \ianc{k _{\ell-1}}{\ell}) \prod_{j \neq k_{\ell-1}} \sum_{k = 1}^N \weight{k}{\ell-1} \delta_k(\rmd \ianc{j}{\ell}) \bwopt{\ell-1}{y}{\iparticle{\ianc{j}{\ell}}{\ell}}{\rmd \iparticle{j}{\ell-1}} \\
& \hspace{2cm} \times \delta_{z_0}(\rmd \lstate^{k_0} _{0}) \delta_{k_1}(\rmd \ianc{k_0}{1}) \prod_{j \neq k_0} \sum_{k = 1}^N \weight{k}{0} \delta_k(\rmd \ianc{j}{1}) \bbwopt{0}{y}{\iparticle{\ianc{j}{1}}{1}}{\rmd \lbottomstate^j _0} \delta_\lmeas(\rmd \ltopstate_0)\eqsp.
\end{align*}

Note however that for all $(k_{0:n}, \ell_{0:n}) \in ([1:N]^{n+1})^2$,
\[
\pE^{k_{0:n}} _{z_{0:n}} \left[ \frac{1}{\prod_{s = 0}^{n-1} \sumweight{s}} \right] = \pE^{\ell_{0:n}} _{z_{0:n}} \left[ \frac{1}{\prod_{s = 0}^{n-1} \sumweight{s}} \right]
\]

and thus it follows that

\begin{equation}
\bfilter{0:n}{N}{\rmd z_{0:n}} = \pE _{z_{0:n}} \left[ \frac{N^{n} \normconst_0 / \normconst_n}{\prod_{s = 0}^{n-1} \sumweight{s}} \right] \filter{0:n}{y}{\rmd z_{0:n}} \eqsp.
\end{equation}

Denote by $\{ \filtration{s} \}_{s = 0} ^n$ the filtration generated by a conditional particle cloud sampled from the kernel $\lawcSMC{}$ (<ref>), i.e. for all $\ell \in [0:n-1]$

\[
\filtration{s} = \sigma\big( \particle^{1:N} _{s:n}, \anc{1:N}{s+1:n} \big) \eqsp.
\]

and $\filtration{n} = \sigma\big( \particle^{1:N} _{n} \big)$. Define for all bounded $f$ and $\ell \in [0:n-1]$

\begin{equation}
\normconstest{\ell:n}(f) = \left\{ \prod_{s = \ell + 1}^{n-1} N^{-1} \sumweight{s} \right\} N^{-1} \sum_{k = 1}^N \uweight{}{\ell}(\particle^{k} _{\ell+1}) f(\particle^{k} _{\ell+1}) \eqsp,
\end{equation}

with the convention $\normconstest{\ell:n}(f) = 1$ if $\ell \geq n$.
Define also the transition Kernel

\begin{equation}
\Qker{\ell-1|\ell+1}{y}: \R^{\dimx} \times \sigalgebra(\R^{\dimx}) \ni (x_{\ell+1}, A) \mapsto \int \1_A(x_\ell) \uweight{}{\ell-1}(\lstate_\ell) \bwopt{\ell}{y}{\lstate_{\ell+1}}{\rmd \lstate_{\ell}}  \eqsp.
\end{equation}

Using <ref>, it is easily seen that for all $\ell \in [0:n-1]$,

\begin{equation}
\label{eq:Qidentity}
\uweight{}{\ell}(\lstate_{\ell+1}) \Qker{\ell-1|\ell+1}{y}(f)(\lstate_{\ell+1}) = \frac{1}{\tfwtrans{\ell+1|0}{\lmeas}{\ltopstate_{\ell+1}}} \int \tfwtrans{\ell|0}{\lmeas}{\ltopstate_s} \uweight{}{\ell-1}(x_\ell) f(\lstate_\ell) \bwtrans{}{\ell}{x_{\ell+1}}{\rmd x_\ell} \eqsp.
\end{equation}

Define $\mathbf{1}: \lstate \in \R^{\dimx} \mapsto 1$. We may thus write that $\normconstest{\ell:n}(f) = N^{-1} \normconstest{\ell+1:n}(\mathbf{1}) \sum_{k = 1}^N \uweight{}{\ell}(\particle^{k} _{\ell+1}) f(\particle^{k} _{\ell+1}).$

For all $\ell \in [0:n-1]$ it holds that

\begin{equation*}
\pE_{z_{0:n}} \big[ \normconstest{\ell-1:n}(f) \big] = \frac{N-1}{N} \pE_{z_{0:n}} \left[ \normconstest{\ell:n}\left( \Qker{\ell-1|\ell+1}{y}(f)\right)\right] + \frac{1}{N} \pE_{z_{0:n}} \left[ \normconstest{\ell:n}(\mathbf{1}) \right] \uweight{}{\ell-1}( z_{\ell}) f(z_{\ell}) \eqsp.
\end{equation*}

By the tower property and the fact that $\normconstest{\ell:n}(f)$ is $\filtration{\ell+1}$-measurable, we have that

\begin{equation*}
\pE_{z_{0:n}} \big[ \normconstest{\ell-1:n}(f) \big] = \pE_{z_{0:n}} \left[ N^{-1} \normconstest{\ell+1:n}(\mathbf{1}) \sumweight{\ell} \pE_{z_{0:n}} \left[ N^{-1} \sum_{k = 1}^N \uweight{}{\ell-1}(\particle^{k} _{\ell}) f(\particle^{k} _{\ell})\bigg| \filtration{\ell+1} \right]\right] \eqsp.
\end{equation*}

Note that for all $\ell \in [0:n-1]$, $(\particle^{1} _{\ell}, \dotsc, \particle^{N-1} _\ell)$ are identically distributed conditionally on $\filtration{\ell+1}$ and

\[
\pE_{z_{0:n}} \left[\uweight{}{\ell-1}(\particle^{j} _{\ell}) f(\particle^{j} _{\ell}) \bigg| \filtration{\ell+1} \right] = \frac{1}{\sumweight{\ell}} \sum_{k = 1}^N \uweight{}{\ell}(\particle^{k} _{\ell+1}) \int \uweight{}{\ell-1}(x_\ell) f(x_{\ell }) \bwopt{\ell}{y}{\particle^{k} _{\ell+1}}{\rmd x_{\ell}} \eqsp,
\]

leading to

\begin{multline*}
\pE_{z_{0:n}} \left[ N^{-1} \sum_{k = 1}^N \uweight{}{\ell-1}(\particle^{k} _{\ell}) f(\particle^{k} _{\ell}) \bigg| \filtration{\ell+1} \right] \\
= \frac{N-1}{N \sumweight{\ell}} \sum_{k = 1}^N \uweight{}{\ell}(\particle^{k} _{\ell+1}) \int \uweight{}{\ell-1}(x_\ell) f(x_{\ell }) \bwopt{\ell}{y}{\particle^{k} _{\ell+1}}{\rmd x_{\ell}} + \frac{1}{N} \uweight{}{\ell-1}(z_\ell) f(z_{\ell}) \eqsp,
\end{multline*}

and the desired recursion follows.

We proceed by induction and show for all $\ell \in [0:n-2]$

\begin{equation}
\begin{aligned}
\label{eq:normconstest_exp}
& \pE_{z_{0:n}}\big[ \normconstest{\ell:n}(f) ] \\
& = \left( \frac{N-1}{N} \right)^{n-\ell} \frac{\int \bwmarg{\ell+1}{\rmd x_{\ell+1}} \tfwtrans{\ell+1|0}{\lmeas}{\ltopstate_{\ell+1}} \uweight{}{\ell}(x_{\ell+1}) f(x_{\ell+1})}{\normconst_{n}}  \\
& + \frac{(N-1)^{n-\ell-1}}{N^{n-\ell}} \bigg[ (\normconst_{\ell+1} / \normconst_n)f(z_{\ell+1}) \uweight{}{\ell}(z_{\ell+1}) \\
& + \sum_{s=\ell+2} ^n \frac{\normconst_s / \normconst_n}{\tfwtrans{s|0}{\lmeas}{\overline{z}_s}} \int \uweight{}{\ell}(x_{\ell+1}) \tfwtrans{\ell+1|0}{\lmeas}{\ltopstate_{\ell+1}} f(x_{\ell+1}) \bwtrans{}{\ell+1|s}{z_s}{\rmd x_{\ell+1}}\bigg] + \frac{\mathsf{D}^\lmeas _{\ell:n}}{N^2} \eqsp.
\end{aligned}
\end{equation}

where $f$ is a bounded function and  $\mathsf{D}^\lmeas _{\ell:n}$ is a a positive constant.
The desired result in <Ref> then follows by taking $\ell = 0$ and $f = \mathbf{1}$.

Assume that (<ref>) holds at step $\ell$. To show that it holds at step $\ell-1$ we use <Ref> and we compute $\pE_{z_{0:n}} \left[ \normconstest{\ell:n}\left( \Qker{\ell-1|\ell+1}{y}(f)\right)\right]$ and $\pE_{z_{0:n}} \left[ \normconstest{\ell:n}(\mathbf{1}) \right] \uweight{}{\ell-1}( z_{\ell}) f(z_{\ell})$.

Using the following identities which follow from (<ref>)

\begin{multline*}
\int  \tfwtrans{\ell+1|0}{\lmeas}{\ltopstate_{\ell+1}} \uweight{}{\ell}(x_{\ell+1}) \Qker{\ell-1|\ell+1}{y}(f)(x_{\ell+1}) \bwmarg{\ell+1}{\rmd x_{\ell+1}}\\
= \int  \tfwtrans{\ell|0}{\lmeas}{\ltopstate_\ell} \uweight{}{\ell-1}(x_\ell) f(x_\ell) \bwmarg{\ell}{\rmd x_\ell}\eqsp,
\end{multline*}

\begin{multline*}
\int \uweight{}{\ell}(x_{\ell+1}) \tfwtrans{\ell+1|0}{\lmeas}{\ltopstate_{\ell+1}} \Qker{\ell-1|\ell+1}{y}(f)(x_{\ell+1}) \bwtrans{}{\ell+1|s}{x_s}{\rmd x_{\ell+1}} \\
= \int \uweight{}{\ell-1}(x_\ell)  \tfwtrans{\ell|0}{\lmeas}{\ltopstate_\ell} f(x_\ell) \bwtrans{}{\ell|s}{x_s}{\rmd x_\ell}\eqsp,
\end{multline*}

we get by (<ref>) that

\begin{equation}
\begin{aligned}
& \frac{N-1}{N} \pE_{z_{0:n}} \left[ \normconstest{\ell:n}\left( \Qker{\ell-1|\ell+1}{y}(f)\right)\right]  \\
& \hspace{.7cm} = \left( \frac{N-1}{N} \right)^{n-\ell + 1} \frac{\int \tfwtrans{\ell|0}{\lmeas}{\ltopstate_\ell} \uweight{}{\ell-1}(x_\ell) f(x_\ell) \bwmarg{\ell}{\rmd x_\ell}}{\normconst_n} \\
& \hspace{1.2cm} + \frac{(N-1)^{n-\ell}}{N^{n-\ell+1}} \bigg[ \frac{\normconst_{\ell+1} / \normconst_{n}}{\tfwtrans{\ell+1|0}{\lmeas}{\overline{z}_{\ell+1}}} \int  \tfwtrans{\ell|0}{\lmeas}{\ltopstate_\ell} \uweight{}{\ell-1}(x_\ell) f(x_\ell) \bwtrans{}{\ell}{z_{\ell+1}}{\rmd x_\ell}\\
& \hspace{1.2cm} + \sum_{s=\ell+2} ^n \frac{\normconst_s / \normconst_n}{\tfwtrans{s|0}{\lmeas}{\overline{z}_s}} \int \uweight{}{\ell-1}(x_\ell)  \tfwtrans{\ell|0}{\lmeas}{\ltopstate_\ell} f(x_\ell) \bwtrans{}{\ell|s}{z_s}{\rmd x_\ell} \bigg] + \frac{\mathsf{D}^\lmeas _{\ell:n}}{N^2}\\
& \hspace{.7cm} = \left( \frac{N-1}{N} \right)^{n-\ell + 1} \frac{\int  \tfwtrans{\ell|0}{\lmeas}{\ltopstate_\ell} \uweight{}{\ell-1}(x_\ell) f(x_\ell) \bwmarg{\ell}{\rmd x_\ell}}{\normconst_n} \\
&  \hspace{1.2cm}  + \frac{(N-1)^{n-\ell}}{N^{n-\ell+1}} \sum_{s=\ell+1} ^n \frac{\normconst_s / \normconst_n}{\tfwtrans{s|0}{\lmeas}{\overline{z}_s}} \int \uweight{}{\ell-1}(x_\ell)  \tfwtrans{\ell|0}{\lmeas}{\ltopstate_s} f(x_\ell) \bwtrans{}{\ell|s}{z_s}{\rmd x_\ell} + \frac{\mathsf{D}^\lmeas _{\ell:n}}{N^2} \eqsp.
\end{aligned}
\end{equation}

The induction step is finished by using again (<ref>) and noting that

\begin{align*}
\frac{1}{N} \pE_{z_{0:n}} \left[ \normconstest{\ell:n}(\mathbf{1}) \right] \uweight{}{\ell-1}( z_{\ell}) f(z_{\ell}) = \frac{(N-1)^{n-\ell}}{N^{n - \ell +1}} \big(\normconst_\ell / \normconst_n\big) \uweight{}{\ell-1}( z_{\ell}) f(z_{\ell}) + \frac{\widetilde{\mathsf{D}}^\lmeas _{\ell:n}}{N^2} \eqsp.
\end{align*}

and then setting $\mathsf{D}^\lmeas _{\ell-1:n} = \mathsf{D}^\lmeas _{\ell:n} + \widetilde{\mathsf{D}}^\lmeas _{\ell:n}$.

It remains to compute the initial value at $\ell = n-2$. Note that

\begin{equation}
\pE_{z_{0:n}} \left[ \normconstest{n-1:n}(f) \right] = \frac{N-1}{N} \int \bwopt{n}{y}{}{}(\rmd x_n) \uweight{}{n-1}(x_n) f(x_n) + \frac{1}{N} \uweight{}{n-1}(z_n) f(z_n)
\end{equation}

and thus by <Ref> and similarly to the previous computations

\begin{align*}
& \pE_{z_{0:n}} \left[ \normconstest{n-2:n}(f) \right] \\
&  = \left( \frac{N-1}{N} \right)^2 \int \bwopt{n}{y}{}{}(\rmd x_n) \uweight{}{n-1}(x_n) \Qker{n-2|n}{y}(f)(x_n) + \frac{N-1}{N^2} \bigg[ \uweight{}{n-1}(z_n)\Qker{n-2|n}{y}(f)(z_n) \\
& +  \uweight{}{n-2}(z_{n-1}) f(z_{n-1}) \int \bwopt{n}{y}{}{}(\rmd x_n) \uweight{}{n-1|n}(x_n) \bigg] + \frac{\mathsf{D}^\lmeas _{n-2fa:n}}{N^2} \\
& = \left( \frac{N-1}{N} \right)^2  \frac{\int \tfwtrans{n-1|0}{\lmeas}{\lstate_{n-1}}  \uweight{}{n-2}(\lstate_{n-1}) \bwmarg{n-1}{\rmd \lstate_{n-1}}}{\normconst_n} \\
& + \frac{N-1}{N^2} \bigg[ \big( \normconst_{n-1} / \normconst_n \big)\uweight{}{n-2}(z_{n-1}) f(z_{n-1}) \\
& + \frac{1}{\tfwtrans{n|0}{\lmeas}{\ltopstate_n}} \int \tfwtrans{n-1|0}{\lmeas}{\ltopstate_{n-1}} \uweight{}{n-2}(\lstate_{n-1}) f(\lstate_{n-1}) \bwtrans{}{n-1}{z_n}{\rmd \lstate_{n-1}}\bigg] + \frac{\mathsf{D}^\lmeas _{n-2:n}}{N^2} \eqsp.
\end{align*}

§.§ Proof of <Ref>

In this section and only in this section we make the following assumption

For all $s \in [0:n-1]$,
\bwmargV{s}{\lstate_s} \fwtrans{s+1}{\lstate_s}{\lstate_{s+1}} = \bwmargV{s+1}{\lstate_{s+1}} \bwtransV{}{s}{\lstate_{s+1}}{\lstate_s} \eqsp.

We also consider $\sigma_{\delta} = 0$. In what follows we let $\tau_{\dimtr+1} = n$ and we write $\tau_{1:\dimtr} = \{ \tau_1, \dotsc, \tau_\dimtr \}$ and $\overline{\tau_{1:\dimtr}} = [1:n] \setminus \tau_{1:t}$. Define the measure

\begin{equation}
\label{eq:joint_noiseless}
\Gamma^{\diffltrmeas} _{0:n}(\rmd \ltrstate_{0:n}) = \bwmargV{n}{\rmd \ltrstate_n} \prod_{s \in \overline{\tau_{1:\dimtr}}} \bwtransV{}{s}{\ltrstate_{s+1}}{\rmd \ltrstate_{s}}  \prod_{i = 1}^{\dimtr} \bwtransV{}{\tau_i}{\ltrstate_{\tau_i + 1}}{\ltrstate_{\tau_i}} \rmd \projidx{\ltrstate}{\tau_i}{i} \delta_{\vecidx{\diffltrmeas}{}{i}} (\rmd \vecidx{\ltrstate}{\tau_i}{i}) \eqsp.
\end{equation}

Under <ref> it has the following alternative forward expression,

\begin{equation}
\label{eq:joint_noiseless:forward}
\Gamma^\diffltrmeas _{0:n}(\rmd \chunk{\ltrstate}{0}{n})  = \bwmargV{0}{\rmd \ltrstate_0} \prod_{s \in \overline{\tau_{1:\dimtr}}} \fwtrans{s+1}{\ltrstate_{s}}{\rmd \ltrstate_{s+1}} \prod_{i = 1}^\dimtr \fwtrans{\tau_i}{\ltrstate_{\tau_i - 1}}{\ltrstate_{\tau_i}} \rmd \projidx{\ltrstate}{\tau_i}{i} \delta_{\vecidx{\diffltrmeas}{}{i}} (\rmd \vecidx{\ltrstate}{\tau_i}{i}) \eqsp.
\end{equation}

Since the forward kernels decompose over the dimensions of the states, i.e.
$$\fwtrans{s+1}{\ltrstate_s}{\ltrstate_{s+1}} = \prod_{\ell = 1}^{\dimx} \idxfwtrans{s+1}{\vecidx{\ltrstate}{s}{\ell}}{\vecidx{\ltrstate}{s+1}{\ell}}{\ell}$$
where $\idxfwtrans{s+1}{\vecidx{\ltrstate}{s}{\ell}}{\vecidx{\ltrstate}{s+1}{\ell}}{\ell} = \normpdf(\vecidx{\ltrstate}{s+1}{\ell}; (\alpha_{s+1} / \alpha_s)^{1/2} \vecidx{\ltrstate}{s}{\ell}, 1 - (\alpha _{s+1} / \alpha_s) )$, we can write

\begin{equation}
\label{eq:joint_noiseless:forward:dim_decomp}
\Gamma^{\diffltrmeas} _{0:n}(\ltrstate_{0:n}) = \bwmargV{0}{\ltrstate_0} \prod_{\ell = 1}^{\dimx} \Gamma^{\diffltrmeas} _{1:n|0, \ell}\big(\vecidx{\ltrstate}{1}{\ell}, \dotsc, \vecidx{\ltrstate}{n}{\ell} \big| \vecidx{\ltrstate}{0}{\ell} \big) \eqsp,
\end{equation}

where for $\ell \in [1:\dimtr]$

\begin{equation}
\label{eq:joint_noiseless:forward_observed}
\Gamma^{\diffltrmeas} _{1:n|0, \ell}\big(\vecidx{\ltrstate}{1}{\ell}, \dotsc, \vecidx{\ltrstate}{n}{\ell} | \vecidx{\ltrstate}{0}{\ell} \big) = \idxfwtrans{\tau_\ell}{\vecidx{\ltrstate}{\tau_{\ell} - 1}{\ell}}{\vecidx{\diffltrmeas}{}{\ell}}{\ell}  \prod_{s \neq \tau_\ell}  \idxfwtrans{s}{\vecidx{\ltrstate}{s-1}{\ell}}{\rmd \vecidx{\ltrstate}{s}{\ell}}{\ell}  \eqsp,
\end{equation}

and for $\ell \in [\dimtr+1:\dimx]$,

\begin{equation}
\label{eq:joint_noiseless:forward_unobserved}
\Gamma^{\diffltrmeas} _{1:n|0, \ell}(\vecidx{\ltrstate}{1}{\ell}, \dotsc, \vecidx{\ltrstate}{n}{\ell} |\vecidx{\ltrstate}{0}{\ell}) = \prod_{s = 0}^{n-1} \idxfwtrans{s+1}{\vecidx{\ltrstate}{s}{\ell}}{\vecidx{\ltrstate}{s+1}{\ell}}{\ell} \eqsp.
\end{equation}

With these quantities in hand we can now prove <Ref>.

Note that for $\ell \in [1:\dimtr]$,

\begin{multline*}
\normpdf(\vecidx{\diffltrmeas}{}{\ell}; \alpha_{\tau_\ell} \vecidx{\ltrstate}{0}{\ell}, 1 - \alpha_{\tau_\ell}) = \idxfwtrans{\tau_\ell | 0}{\vecidx{\ltrstate}{0}{\ell}}{\vecidx{\diffltrmeas}{}{\ell}}{\ell}
=  \int \idxfwtrans{\tau_\ell}{\vecidx{\ltrstate}{\tau_{\ell} - 1}{\ell}}{\vecidx{\diffltrmeas}{}{\ell}}{\ell}  \prod_{s \neq \tau_\ell}  \idxfwtrans{s}{\vecidx{\ltrstate}{s-1}{\ell}}{\rmd \vecidx{\ltrstate}{s}{\ell}}{\ell}
\\ = \int \Gamma^{\diffltrmeas} _{1:n|0, \ell}\big(\rmd (\vecidx{\ltrstate}{1}{\ell}, \dotsc, \vecidx{\ltrstate}{n}{\ell}) |\vecidx{\ltrstate}{0}{\ell} \big)
\end{multline*}

and thus by (<ref>) we have that

\begin{align*}
\bwmargV{0}{\ltrstate_0} \pot{0}{\lmeas}{\ltrstate_0} & \propto \bwmargV{0}{\ltrstate_0} \prod_{\ell = 1}^\dimtr \normpdf(\vecidx{\diffltrmeas}{}{\ell}; \alpha_{\tau_\ell} \vecidx{\ltrstate}{0}{\ell}, 1 - \alpha_{\tau_\ell}) \\
& = \bwmargV{0}{\ltrstate_0} \prod_{\ell = 1}^\dimtr \int \Gamma^{\diffltrmeas} _{1:n|0, \ell}\big(\rmd (\vecidx{\ltrstate}{1}{\ell}, \dotsc, \vecidx{\ltrstate}{n}{\ell}) |\vecidx{\ltrstate}{0}{\ell} \big)  \\
& = \bwmargV{0}{\ltrstate_0} \prod_{\ell = 1}^\dimx \int \Gamma^{\diffltrmeas} _{1:n|0, \ell}\big(\rmd (\vecidx{\ltrstate}{1}{\ell}, \dotsc, \vecidx{\ltrstate}{n}{\ell}) |\vecidx{\ltrstate}{0}{\ell} \big) \eqsp.
\end{align*}

By (<ref>) it follows that

\[
\filter{0}{\ltrmeas}{\ltrstate_0} = \frac{1}{\int \Gamma^{\diffltrmeas} _{0:n}(\chunk{\tilde\ltrstate}{0}{n}) \rmd \chunk{\tilde\ltrstate}{0}{n} }\int \Gamma^{\diffltrmeas} _{0:n}(\chunk{\ltrstate}{0}{n}) \rmd \chunk{\ltrstate}{1}{n} \eqsp,
\]

and hence by (<ref>) and (<ref>) we get

\[
\noisyfilter{0}{\ltrmeas}{\ltrstate_{0}} \propto
\int \bwmargV{\tau_\dimtr}{\ltrstate_{\tau_\dimtr}} \delta_{\vecidx{\diffltrmeas}{}{\dimtr}}(\rmd \vecidx{\ltrstate}{\tau_\dimtr}{\dimtr}) \rmd \projidx{\ltrstate}{\tau_\dimtr}{\dimtr} \left\{ \prod_{i = 1}^{\dimtr-1} \bwtransV{}{\tau_i | \tau_{i + 1}}{\ltrstate_{\tau_{i+1}}}{\ltrstate_{\tau_{i}}} \delta_{\vecidx{\diffltrmeas}{}{i}}(\rmd \vecidx{\ltrstate}{\tau_i}{i}) \rmd \projidx{\ltrstate}{\tau_i}{i} \right\} \bwtransV{}{0|\tau_1}{\ltrstate_{\tau_1}}{\ltrstate_0} \eqsp. \]

This completes the proof.

§.§ Algorithmic details and numerics