# Discrete Markov Probabilistic Models

Le-Tuyet-Nhi Pham, Dario Shariatian, Antonio Ocello, Giovanni Conforti, Alain Durmus

► To cite this version:

Le-Tuyet-Nhi Pham, Dario Shariatian, Antonio Ocello, Giovanni Conforti, Alain Durmus. Discrete Markov Probabilistic Models. 2025. hal-04951334

HAL Id: hal-04951334

https://hal.science/hal-04951334v1

Preprint submitted on 17 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Discrete Markov Probabilistic Models

Le-Tuyet-Nhi Pham¹, Dario Shariatian², Antonio Ocello¹, Giovanni Conforti³, and
Alain Durmus¹

¹Centre de Mathématiques Appliquées (CMAP), École Polytechnique, France

²Sierra lab, Inria, Paris, France

³Department of Mathematics, University of Padova, Italy

## Abstract

This paper introduces the Discrete Markov Probabilistic Model (DMPM), a novel algorithm for discrete data generation. The algorithm operates in the space of bits {0, 1}<sup>d</sup>, where the noising process is a continuous-time Markov chain that can be sampled exactly via a Poissonian clock that flips labels uniformly at random. The time-reversal process, like the forward noise process, is a jump process, with its intensity governed by a discrete analogue of the classical score function. Crucially, this intensity is proven to be the conditional expectation of a function of the forward process, strengthening its theoretical alignment with score-based generative models while ensuring robustness and efficiency. We further establish convergence bounds for the algorithm under minimal assumptions and demonstrate its effectiveness through experiments on low-dimensional Bernoulli-distributed datasets and high-dimensional binary MNIST data. The results highlight its strong performance in generating discrete structures. This work bridges theoretical foundations and practical applications, advancing the development of effective and theoretically grounded discrete generative modeling.

## Introduction

Score-based Generative Models (SGMs) have become a key reference for generating complex data, such as images (see, e.g., Rombach et al., 2022; Ramesh et al., 2022; Saharia et al., 2022), audio (Chen et al., 2020; Kong et al., 2020), and video (Ho et al., 2022; Villegas et al., 2022; Bar-Tal et al., 2024). In continuous time, this approach benefits from a strong theoretical framework, and a scalable, stable learning objective.

By contrast, discrete generative modeling continues to pose significant challenges. Multiple diffusion-based methods have recently been proposed for discrete spaces (Austin et al., 2021; Hoogeboom et al., 2021; Shi et al., 2024; Campbell et al., 2022; Holderrieth et al., 2024; Ren et al., 2024), or spaces of mixed type (Bertazzi et al., 2024); but there is still no consensus on which approach is theoretically sound or most practically efficient. Various formulations rely on complex forward kernels or computationally unstable ratio-based estimators for backward transitions, leading to limited convergence guarantees and high computational costs in high dimensions. Furthermore, recent analyses of discrete diffusions have introduced valuable theoretical tools (Campbell et al., 2022; Holderrieth et al., 2024; Ren et al., 2024), yet most methods remain either overly generic or
require strong assumptions, making them difficult to scale or to deploy with simple, stable training objectives.

**Contributions.** In this paper, we introduce Discrete Markov Probabilistic Models (DMPMs), a new class of generative models for discrete data that bridges these gaps. Our framework specializes the forward noising process to a continuous-time Markov chain on the hypercube $\{0, 1\}^d$. Leveraging theoretical insights on time-reversal Markov dynamics of this process, this choice preserves the key strengths and structure of continuous SGMs, addressing the issues raised in prior work. Our main results are summarized as follows:

* **Forward-Backward Construction.** We provide a principled derivation of the noising (forward) and denoising (backward) processes.

* **Score function and stable estimation.** Our analysis reveals how the time-reversal inherit a score function with an explicit conditional expectation form. By casting learning as an $L^2$ projection, we eliminate the need for numerically unstable ratios of transition probabilities. This leads to a robust training procedure.

* **Theoretical Guarantees.** We prove that DMPMs converge to the underlying data distribution under minimal assumptions, providing non-asymptotic error bounds that underscore the method's reliability.

* **Empirical Performance.** We demonstrate that our approach attains competitive or superior performance on discrete datasets, including binarized MNIST, frequently with fewer function evaluations compared to existing discrete diffusion frameworks (e.g., 2.89 vs 7.34 FID compared to Discrete Flow Matching, Gat et al., 2024, with 2.5x fewer network calls).

**Notation.** Given a measurable space $(E, \mathcal{E})$, we denote by $\mathcal{P}(E)$ the set of probability measures on $E$. Given two probability measures $\mu, \nu \in \mathcal{P}(E)$, the Kullback-Leibler divergence (also called relative entropy) of $\mu$ with respect to $\nu$ is defined as $\mathrm{KL}(\mu|\nu) := \int \log(d\mu/d\nu)d\mu$ if $\mu$ is absolutely continuous with respect to $\nu$, and $\mathrm{KL}(\mu|\nu) = +\infty$ otherwise. The total variation distance of $\mu$ and $\nu$ is defined as $\|\mu - \nu\|_{\mathrm{TV}} = \int |d\mu/dR - d\nu/dR|RdR$ for any $R \in \mathcal{P}(E)$ such that $\mu$ and $\nu$ are absolutely continuous with respect to $R$. Consider a random variable $X$, we denote by $\mathrm{Law}(X)$ the law of $X$.

# 1 Forward and backward process of DMPMs

We introduce a generative modeling framework that adapts classical diffusion-based methods to discrete state spaces. Let $(\overrightarrow{X}_t)_{t \in [0,T]}$ be a forward Markov process, starting from $\overrightarrow{X}_0 \sim \mu^*$, the data distribution of interest, and evolving over the fixed time horizon $T_f > 0$. This forward process smoothly transports $\mu^*$ into a simple base distribution, at time $T_f$. In the continuous setting, a well-known example is the Ornstein-Uhlenbeck process, which converges to a Gaussian base measure.

Given the forward dynamics, we define the corresponding backward process $(\overleftarrow{X}_t)_{t \in [0,T_f]}$ by
$$ \overleftarrow{X}_t := \overrightarrow{X}_{T_f-t}. $$
By construction, $(\overleftarrow{X}_t)$ evolves in reverse from the base distribution (at $t=0$) back to the original data distribution $\mu^*$ (at $t=T_f$). Although $(\overleftarrow{X}_t)$ is also Markovian in many classical settings, its transition rates or drift terms are often intractable, preventing direct simulation. Instead, most generative models leverage the fact that the backward dynamics can be characterized and learned from forward-simulated paths. In the continuous setting, this learning is achieved via score-matching methods.
Following this principle, we define a simple forward continuous time Markov chain (CTMC) on $\{0,1\}^d$ where bits flip according to a Poisson clock, and show its time-reversal remains a tractable CTMC. We derive closed-form expressions for the backward transition rates, which involve conditional expectations over the forward process. This enables, for the first time, an efficient training procedure based on regression in the discrete setting.

## 1.1 Simplest case $X = \{0,1\}$

To introduce the key ideas, we first consider the simplest case $X = \{0,1\}$. The forward process $(\vec{X}_t)_{t \in [0,T_f]}$ starting from $\vec{X}_0 \sim \mu^*$, is defined as follows. Consider the fixed jump times $(T_i)_{i \in \{1, \dots, N\}} | N \stackrel{\text{iid}}{\sim} \text{Unif}([0, T_f])$ of a Poisson process over $[0, T_f]$ where $N \sim \text{Pn}(\lambda T_f)$ is the number of jump, and $\lambda > 0$ is a prescribed jump rate. Without loss of generality, we assume that $0 = T_0 \le T_1 < \dots < T_N$. We define recursively $(\vec{X}_t)_{t \in [0, T_f]}$ over $(T_i, T_{i+1}]$. Suppose that $\vec{X}_{T_i}$ has been defined we set $\vec{X}_t = \vec{X}_{T_i}$ for any $t \in (T_i, T_{i+1}]$ and $\vec{X}_{T_{i+1}} = 1 - \vec{X}_{T_i}$. It is well known that $(\vec{X}_t)_{t \in [0, T_f]}$ is a Markov jump process (Owen, 2013, Section 6) with generator $\vec{q}$ defined for any $x, y \in X$ as

$$ \vec{q}(x, y) := \begin{cases} \lambda, & \text{if } y \neq x, \\ -\lambda, & \text{otherwise.} \end{cases} \tag{1} $$

The transition probability matrix $P(\vec{X}_t = y | \vec{X}_0 = x) = \vec{p}_t^1(x, y)$, for $x, y \in X, 0 \le t \le T_f$, is known to be

$$ \vec{p}_t^1(x, y) = \begin{cases} \frac{1}{2} + \frac{1}{2}e^{-2\lambda t}, & \text{if } x = y, \\ \frac{1}{2} - \frac{1}{2}e^{-2\lambda t}, & \text{otherwise.} \end{cases} \tag{2} $$

The proof of this result is given in the supplementary material B.1. To recover the data distribution, we analyze the time-reversed process, which is denoted by $(\overleftarrow{X}_t)_{t \in [0, T_f]}$, and defined as $\overleftarrow{X}_t = \overleftarrow{X}_{T_f - t}$ for any $t \in [0, T_f]$. Conforti and Léonard (2022, Theorem 2.8) shows that $(\overleftarrow{X}_t)_{t \in [0, T_f]}$ is also a non-homogeneous CTMC, i.e., it is associated with a family of generator matrices $(\overleftarrow{q}_t)_{t \in [0, T_f]}$ which satisfies the time-reversal formula: for any $0 \le t \le T_f$ and $x, y \in X$,

$$ \mu_{T_f-t}(x) \overleftarrow{q}_t(x, y) = \mu_{T_f-t}(y) \overleftarrow{q}(y, x). \tag{3} $$

Since $\overleftarrow{q}$ is symmetric (see (1)) and $\mu_{T_f-t}(x) > 0$ for all $x \in X, t \in [0, T_f)$, we deduce that the backward generator $\overleftarrow{q}_t$ for $0 \le t < T_f$ is given for any $x, y \in X$ by

$$ \overleftarrow{q}_t(x, y) = \overleftarrow{q}(x, y) \frac{\mu_{T_f-t}(y)}{\mu_{T_f-t}(x)}, \tag{4} $$

with $\mu$ the forward marginal distribution satisfies $\mu_0 = \mu^*$. Having access to $(\overleftarrow{q}_t)_{t \in [0, T_f]}$ and $\mu_t$ allows to sample from $\overleftarrow{X}_t$ for any $t \in [0, T_f]$ as follows. Define $s_t: X \to \mathbb{R}$ for any $x \in X$ by

$$ s_t(x) := \frac{\mu_{T_f-t}(x) - \mu_{T_f-t}(1-x)}{\mu_{T_f-t}(x)}. \tag{5} $$
$s_t$ acts as a discrete derivative in $\mathcal{X}$ of $\log \mu_t$, and thus serves as a discrete analogue of the score function in continuous models. With this notation, $\overleftarrow{\mathbf{q}}_t(x, y)$ (4) can be expressed, for any $x, y \in \mathcal{X}$, for $0 \le t < T_f$, as:

$$ \overleftarrow{\mathbf{q}}_t(x, y) = -\lambda \mathbf{1}_{\{y=x\}} + \lambda(1 - s_t(x)) \mathbf{1}_{\{y=1-x\}}. $$

Then, starting from a sample $\overleftarrow{\mathbf{X}}_0 \sim \mu_{T_f}$, and a sequence of i.i.d. random variables distributed according to the exponential distribution with parameter 1, $\{E_i : i \in \mathbb{N}\}$, we can define the jump times $(T_i)_{i \in \mathbb{N}}$ of the backward process and its transition by induction setting $T_0 = 0$. Given $(T_i, \overleftarrow{\mathbf{X}}_{T_i})$, we define the next jump time as $T_{i+1} = T_i + \Delta T_{i+1}$, where $\Delta T_{i+1} = \inf\{t \ge 0 : \int_0^t \bar{\lambda}_{T_i+r}(\overleftarrow{\mathbf{X}}_{T_i})dr \ge E_i\}$, where $\bar{\lambda}_t(x) = \lambda(1 - s_t(x))$. Then, set $\overleftarrow{\mathbf{X}}_t = \overleftarrow{\mathbf{X}}_{T_i}$ for $t \in (T_i, T_{i+1} \wedge T_f)$, and finally if $T_{i+1} < T_f$, set $\overleftarrow{\mathbf{X}}_{T_{i+1}} = 1 - \overleftarrow{\mathbf{X}}_{T_i}$.

Eventhough $\overleftarrow{\mathbf{q}}_t$ is intractable, it can be approximated with the function $s_t$. In particular, we show that $s_t$ can be expressed as a conditional expectation over the forward process. For $x \in \mathcal{X}$ and $t \in [0, T_f)$,

$$ s_t(x) = \mathbb{E} \left[ \frac{2\alpha_{T_f-t}}{1 + \alpha_{T_f-t}} - \frac{4\alpha_{T_f-t}(\overrightarrow{\mathbf{X}}_{T_f-t} - \overleftarrow{\mathbf{X}}_0)^2}{1 - \alpha_{T_f-t}^2} \middle| \overrightarrow{\mathbf{X}}_{T_f-t} = x \right], \quad (6) $$

with

$$ \alpha_t := e^{-2\lambda t}. \quad (7) $$

Indeed, with $\overleftarrow{d}_t : (z, x) \mapsto \mathbb{P}[\overrightarrow{\mathbf{X}}_0 = z | \overrightarrow{\mathbf{X}}_{T_f-t} = x]$:

$$ \begin{aligned} s_t(x) &= \frac{\mu_{T_f-t}(x) - \mu_{T_f-t}(1-x)}{\mu_{T_f-t}(x)} \\ &= \sum_{z \in \mathcal{X}} \left(1 - \frac{p_{T_f-t}^1(z, 1-x)}{p_{T_f-t}^1(z, x)}\right) \frac{\mu_0(z) \overrightarrow{p}_{T_f-t}^1(z, x)}{\mu_{T_f-t}(x)} \\ &= \sum_{z \in \mathcal{X}} \frac{2\alpha_{T_f-t}}{1 - \alpha_{T_f-t}^2} (1 - \alpha_{T_f-t} - 2(x-z)^2) \overleftarrow{d}_t(z, x) \\ &= \mathbb{E} \left[ \frac{2\alpha_{T_f-t}}{1 + \alpha_{T_f-t}} - \frac{4\alpha_{T_f-t}(\overrightarrow{\mathbf{X}}_{T_f-t} - \overleftarrow{\mathbf{X}}_0)^2}{1 - \alpha_{T_f-t}^2} \middle| \overrightarrow{\mathbf{X}}_{T_f-t} = x \right]. \end{aligned} $$

Therefore, the function $s$ is an $\mathbb{L}^2$-projection and its approximation boils down to a regression problem.

## 1.2 General state space $\mathcal{X} = \{0, 1\}^d$

### 1.2.1 Forward noising process

We generalize the previous results for the hypercube in $\mathbb{R}^d$, i.e., the state space is $\mathcal{X} = \{0, 1\}^d$ with $d \in \mathbb{N}^*$. We consider the forward homogeneous Markov process $(\overrightarrow{\mathbf{X}}_t)_{t \in [0, T_f]}$ starting from $\overleftarrow{\mathbf{X}}_0 \sim \mu^*$, defined as follows.

We consider the jump times $(T_i)_{i \in \{1, \dots, N\}} | N \stackrel{\text{iid}}{\sim} \text{Unif}([0, T_f])$ of a Poisson process over $[0, T_f]$ where $N \sim \text{Pn}(\lambda T_f)$ is the number of jump. Without loss of generality, we suppose that $T_0 = 0 \le T_1 < \dots < T_N$. We define recursively $(\overrightarrow{\mathbf{X}}_t)_{t \in [0, T_f]}$ over $(T_i, T_{i+1}]$ as follows. Suppose $\overrightarrow{\mathbf{X}}_{T_i}$ has been defined. We set $\overrightarrow{\mathbf{X}}_t = \overrightarrow{\mathbf{X}}_{T_i}$ for $t \in (T_i, T_{i+1})$, and finally, set $\overrightarrow{\mathbf{X}}_{T_{i+1}} = 1 - \overrightarrow{\mathbf{X}}_{T_i}$, where
$\ell_i \sim \text{Unif}(\{1, \dots, d\})$, with $\ell_i$ independent from the past, and $\vec{X}_{T_{i+1}}^j = \vec{X}_{T_i}^j$ for $j \neq \ell_i$. The process $(\vec{X}_t)_{t \in [0, T_f]}$ is a Markov jump process again, associated with the generator $\vec{q}$ defined for any function $g: X \to \mathbb{R}$ as

$$ \vec{q} g(x) = \lambda\{kg(x) - g(x)\}, \qquad (8) $$

where $\lambda > 0$ is a prescribed jump rate and $k$ is the Markov kernel defined as: for any $x, y \in X$,

$$ k(x, y) := \mathbb{1}_{\|x-y\|_2=1} \cdot 1/d. $$

Similarly to the one-dimensional case, we can establish an explicit expression for the transition probability matrix $\vec{p}_t$ for $0 \le t \le T_f$ as

$$ \vec{p}_t(x, y) = \prod_{i=1}^{d} \vec{p}_t^1(x^i, y^i), \qquad (9) $$

where $\vec{p}_t^1$ is defined in (2) and $x = (x^i)_{i=1}^d \in X$ and $y = (y^i)_{i=1}^d \in X$. The detailed computation is given in the supplementary material B.2.1 The factorization of the transition probability in (9) is of great practical interest, as this tells us that the dynamic of the forward process simply consists in the single-bit forward dynamic applied independently to each component, as described in Section 1.1. As a consequence, the forward marginal distribution $\mu_t$ of $\vec{X}_t$ admits the formula

$$ \mu_t(x) = \sum_{z \in X} \mu_0(z) \prod_{i=1}^{d} \vec{p}_t(z, x). \qquad (10) $$

### 1.2.2 Backward process

Denote by $(\vec{X}_t)_{t \in [0, T_f]}$, the time-reversal process associated with $(\vec{X}_t)_{t \in [0, T_f]}$, and defined as $\vec{X}_t = \vec{X}_{T_f-t}$ for any $t \in [0, T_f]$. As in the case $d=1$, Conforti and Léonard (2022, Theorem 2.8) shows that $(\vec{X}_t)_{t \in [0, T_f]}$ is also a non-homogeneous CTMC, with backward generator matrix $(\vec{q}_t)_{t \in [0, T_f]}$ that satisfies (3) and therefore (4), proceeding as before. As in the case $d=1$, we show that $(\vec{q}_t)_{t \in [0, T_f]}$ depends only on a discrete score function, which we now introduce.

First note that (4) and (8) yield $\vec{q}_t(x, y) = 0$, for $x, y \in X$ satisfying $\|x - y\|^2 \neq 1$ and $x \neq y$. Then, for $0 \le t < T_f$, define $s_t: X \to \mathbb{R}^d$ for any $x \in X$, $s_t(x) = \{s_t^\ell(x)\}_{\ell=1}^d$ as the vector in $\mathbb{R}^d$, with components $\ell \in \{1, \dots, d\}$,

$$ s_t^\ell(x) := \frac{\mu_{T_f-t}(x) - \mu_{T_f-t}(\varphi^{(\ell)}(x))}{\mu_{T_f-t}(x)}, \qquad (11) $$

where $\varphi^{(\ell)}: X \to X$ is defined as $\varphi^{(\ell)}(x) = y$, with $y$ obtained by flipping the $\ell$-th bit of $x$, i.e., $y^\ell = 1 - x^\ell$, and $y^i = x^i$ for $i \neq \ell$. Note that the function $s$ thus defined is an extension to the case $d \ge 1$ of the function $s$ defined for $d=1$ in (5). As a result, $s_t$ is a conditional expectation over the forward process, where each of its components admits an expression similar to the 1d case (6).

**Proposition 1.1.** The score function can be expressed as a conditional expectation:

$$ s_t^\ell(x) = \mathbb{E} \left[ f_t^\ell \left( \vec{X}_0^\ell, \vec{X}_{T_f-t} \right) \middle| \vec{X}_{T_f-t} = x \right], \qquad (12) $$

where $t \in [0, T_f)$, $x \in X$, $\ell = 1, \dots, d$, $s_t^\ell$ is the $\ell$-th component of the score function $s_t$, and

$$ f_t^\ell(\vec{X}_0^\ell, \vec{X}_{T_f-t}) = \frac{2\alpha_{T_f-t}}{1 + \alpha_{T_f-t}} - \frac{4\alpha_{T_f-t}(\vec{X}_{T_f-t}^\ell - \vec{X}_0^\ell)^2}{1 - \alpha_{T_f-t}^2}. \qquad (13) $$
The proof of this result is given in Appendix B.2.2.

Furthermore, for $0 \le t < T_f$, $x, y \in \mathbf{X}$, we can write the backward generator $\overleftarrow{q}_t(x, y)$, as given in (3), as:

$$ \overleftarrow{q}_t(x, y) = -\lambda \mathbf{1}_{\{y=x\}} + \sum_{\ell=1}^{d} \lambda(1 - s_t^{\ell}(x)) \mathbf{1}_{\{y=\varphi^{(\ell)}(x)\}} $$

This defines the non-homogeneous jump rate $\overleftarrow{\lambda}_t$ and jump kernel $\overleftarrow{k}_t$ of the backward process:

$$ \overleftarrow{\lambda}_t(x) = \lambda \sum_{\ell=1}^{d} (1 - s_t^{\ell}(x)) \\ \overleftarrow{k}_t(x, y) = \mathbf{1}_{y=\varphi^{(\ell)}(x)} \cdot \lambda(1 - s_t^{\ell}(x)) / \overleftarrow{\lambda}_t(x) \quad (14) $$

for $x, y \in \mathbf{X}$ and $t \in [0, T_f]$. Thus, having access to $(s_t)_{t \in [0, T_f]}$ and $\mu_{T_f}$, for any $t \in [0, T_f]$, allows to sample from $\overleftarrow{X}_t$, as follows. Starting from a sample $\overleftarrow{X}_0$ from $\mu_{T_f}$ and a sequence of i.i.d. random variables distributed according to the exponential distribution with parameter 1, $\{E_i : i \in \mathbb{N}\}$, we can define the jump times $(T_i)_{i \in \mathbb{N}}$ of the backward process and its transition by induction setting $T_0 = 0$. Given $(T_i, \overleftarrow{X}_{T_i})$, we define the next jump time as $T_{i+1} = T_i + \Delta T_{i+1}$, where $\Delta T_{i+1} = \inf\{t \ge 0 : \int_0^t \overleftarrow{\lambda}_{T_i+r}(\overleftarrow{X}_{T_i})dr \ge E_i\}$. Then, set $\overleftarrow{X}_t = \overleftarrow{X}_{T_i}$ for $t \in (T_i, T_{i+1} \wedge T_f)$, and finally if $T_{i+1} < T_f$, $\overleftarrow{X}_{T_{i+1}} = \varphi^{(\ell)}(\overleftarrow{X}_{T_i})$ for $\ell_i \in \{1, \dots, d\}$ which is distributed according to Cate($\{\overleftarrow{k}_{T_{i+1}}(\overleftarrow{X}_{T_i}, \varphi^{(\ell)}(\overleftarrow{X}_{T_i}))\}_{\ell=1}^d$). Another equivalent procedure to generate the backward is provided in the supplement C.1.

## 1.3 Approximating the backward characteristics

Similarly to common diffusion-based generative models, we aim to sample from the time-reversal process. However, exact simulations are not possible and face similar challenges: a) we do not have access to i.i.d. samples from $\mu_{T_f}$, b) the score function of the forward process defined in (11) is intractable.

**Initialize the backward from the uniform distribution.** we show that $(\overrightarrow{X}_t)_{t \in [0, T_f]}$ converges geometrically to $\gamma^d$, the uniform distribution over $\mathbf{X}$ (see B.2.3 in the supplementary document). This should be put in parallel with diffusion-based models, where the stochastic process at hand, e.g., Ornstein-Uhlenbeck, converges geometrically fast to some Gaussian distribution.

**Training procedure to fit the score function.** To address b), we exploit the conditional expectation structure of the score function, as given in Proposition 1.1.

We approximate $(s_t)_{t \in [0, T_f]}$ using a parameterized family $(t, x) \mapsto s_t^{\theta}(x)\theta \in \Theta$, where the parameter $\theta$ is fitted minimizing an adapted score-matching objective $\mathcal{L}_{L^2}$, defined as the function

$$ \theta \mapsto \int_{0}^{T_f} \mathbb{E} \left[ \| s_{T_f-t}(\overrightarrow{X}_t) - f_{T_f-t}(\overrightarrow{X}_0, \overrightarrow{X}_t) \|^2 \right] dt. \quad (15) $$

Another option to fit $\theta$ is to use the fact that for any $x \in \mathbf{X}$, $t \in [0, T_f]$, $\ell \in \{1, \dots, d\}$, $1 - s_t^{\ell}(x)$ is non-negative. Thus, we introduce the following entropy-based term:

$$ \theta \mapsto \int_{0}^{T_f} \mathbb{E} \left[ \sum_{\ell=1}^{d} (1 - s_{T_f-t}^{\theta, \ell}) h \left( \frac{1 - s_{T_f-t}^{\ell}}{1 - s_{T_f-t}^{\theta, \ell}} \right) (\overrightarrow{X}_t) \right] dt, $$
where $h(a) = a \log(a) - (a - 1)$. Minimizing this function is equivalent to minimizing:

$$ \mathcal{L}_e : \theta \mapsto \int_0^{T_f} \mathbb{E} \left[ \sum_{\ell=1}^d \left( -s_{T_f-t}^{\theta, \ell}(\vec{X}_t) + (f_{T_f-t}^{\ell}(\vec{X}_t) - 1) \log(1 - s_{T_f-t}^{\theta, \ell}(\vec{X}_t)) \right) \right] dt. \quad (16) $$

We further derive a discrete-denoiser structure.

**Proposition 1.2.** The score function admits the following discrete denoiser expression:

$$ s_t^{\ell}(x) = \frac{2\alpha_{T_f-t}}{1 + \alpha_{T_f-t}} - \frac{4\alpha_{T_f-t}d_t^{\ell}(x)}{1 - \alpha_{T_f-t}^2}, \quad (17) $$

where $d_t^{\ell}(x) = \mathbb{P}(\vec{X}_0^{\ell} \neq x^{\ell} | \vec{X}_{T_f-t} = x)$ serves as an optimal classifier, referred to as a *discrete denoiser*.

The proof is given in Appendix C.3. We leverage this structure by reparameterizing our score model:

$$ s_t^{\theta, \ell}(x) = \frac{2\alpha_{T_f-t}}{1 + \alpha_{T_f-t}} - \frac{4\alpha_{T_f-t}d_t^{\theta, \ell}(x)}{1 - \alpha_{T_f-t}^2}. \quad (18) $$

As a result, we modify our objective $\mathcal{L}_{L^2}$ to $\mathcal{L}_{L^2}^{\text{denoiser}}$ to fit the conditional expectation $d_t(x)$ rather than $f_t$, see (37) in Appendix C.4. To fit $d_t^{\theta}(x)$ to $d_t(x)$, we introduce an additional cross-entropy loss $\mathcal{L}_{\text{CE}}$:

$$ \mathcal{L}_{\text{CE}} : \theta \mapsto \int_0^{T_f} \mathbb{E} \left[ \sum_{l=1}^d \mathbb{1}_{\vec{X}_0^{\ell} \neq \vec{X}_{T_f-t}^{\ell}} \log d_t^{\theta, \ell}(\vec{X}_{T_f-t}^{\ell}) + (1 - \mathbb{1}_{\vec{X}_0^{\ell} \neq \vec{X}_{T_f-t}^{\ell}}) \log(1 - d_t^{\theta, \ell}(\vec{X}_{T_f-t}^{\ell})) \right] dt. $$

Based on the previous discussions, we consider a linear combination of the losses $\mathcal{L}_{L^2}^{\text{denoiser}}$, $\mathcal{L}_e$, $\mathcal{L}_{\text{CE}}$, respectively weighted by factors $\varpi_1, \varpi_2, \varpi_3$, which results in the loss $\mathcal{L}_{\varpi}$:

$$ \mathcal{L}_{\varpi} = \varpi_1 \mathcal{L}_{L^2}^{\text{denoiser}} + \varpi_2 \mathcal{L}_e + \varpi_3 \mathcal{L}_{\text{CE}}. \quad (19) $$

The expected value of $d_t^{\ell}$ is given by

$$ w_t = \mathbb{E} \left[ d_t^{\ell}(\vec{X}_{T_f-t}) \right] = (1 - \alpha_{T_f-t})/2, \quad (20) $$

as detailed in Appendix C.4. Thus, we scale losses $\mathcal{L}_{L^2}$, $\mathcal{L}_{\text{CE}}$ by $1/w_t$, ensuring a constant average magnitude across timesteps; see (43) and Figure 4 in Appendix C.4. This leads to the updated loss $\mathcal{L}_{\varpi}^w$ (see (45)). Detailed derivations are provided in Appendix C.4. The final training procedure is outlined in Algorithm 2.

## 1.4 Generative process

Suppose we have access to a score approximation $(s_t^{\theta^*})_{t \in [0, T_f]}$. The generative model can then be sampled analogously to the backward process, replacing replacing $\mu_{T_f}$ with $\gamma^d$ and $(s_t)_{t \in [0, T_f]}$ with $(s^{\theta^*} t)_{t \in [0, T_f]}$, leading to the non-homogeneous jump rate and kernel approximating (14): for $x, y \in \mathbf{X}$ and $t \in [0, T_f]$,

$$ \begin{aligned} \lambda_t^{\theta^*}(x) &= \lambda \sum_{\ell=1}^{d} (1 - s_t^{\theta^*, \ell}(x)) \\ k_t^{\theta^*}(x, y) &= \mathbb{1}_{y=\varphi^{(\ell)}(x)} \cdot \lambda(1 - s_t^{\theta^*, \ell}(x)) / \lambda_t^{\theta^*}(x) \end{aligned} \quad (21) $$
where we denote by $s_t^{\theta^*, \ell}$ the $\ell$-th component of $s_t^{\theta^*}$. For completeness, Algorithm 1 in Appendix C.2 provides the pseudo-code for an ideal, continuous-time approximation of the backward process.

In practice, exact integration of the jump rate is infeasible, requiring time discretization in the backward simulation. We approximate the jump rate and kernel using piecewise constant functions $\hat{\lambda}, \hat{k}$ (21), such that, for $x, y \in X$ and $t \in [t_k, t_{k+1})$,

$$ \hat{\lambda}_t^{\theta^*}(x) = \lambda_{t_k}^{\theta^*}(x) \\ \hat{k}_t^{\theta^*}(x, y) = k_{t_k}^{\theta^*}(x, y), \qquad (22) $$

where the discrete time scheme $\{t_k\}_{k=0}^K$ are associated with step-sizes $\{h_k\}_{k=1}^K$, $t_k = \sum_{i=1}^k h_i$ with $t_0 = 0$ and $t_K = T_f$. Based on this jump rate and jump kernel, we can define a new CTMC which can be sampled in practice, following the same procedure as the ideal backward process and starting from an observation from $\gamma^d$. We denote the resulting process $(\overline{\hat{X}}_t^*)_{t \in [0, T_f]}$. Under some assumptions, the approximations of the backward process on both continuous and discrete time schemes have the final law converging toward the desired data distribution. Details are provided in the next section.

The resulting DMPM sampler used to approximate the backward dynamic is provided in Algorithm 3, Appendix C.5. Time discretization strategies, referred to as time-schedules, are listed in Table 2.

Using the index distribution $k_t^{\theta^*}$, we propose a methodological modification to the previous algorithm by flipping $M_{t_k}$ many bits instead of one at each timestep $t_k$. The sequence $\{M_{t_k}\}_{k=1}^K$ is referred to as a flip-schedule, with two natural choices listed in Table 3. The resulting procedure is given in Algorithm 4, Appendix C.5.

Finally, leveraging the discrete denoiser structure from Proposition 1.2, we introduce a denoise- renoise sampler. This method alternates one-step denoising from $0 = t_0$ to $T_f$, followed by re-noising back to $t_1$ etc., resembling multistep sampling in consistency models (Song et al., 2023). The resulting procedure is given in Algorithm 5, Appendix C.5.

We refer the reader to Appendix C.5 for more details.

## 2 Convergence of DMPMs algorithm

This section provides quantitative error estimates between the generated final distribution Law$(\overline{\hat{X}}_{T_f}^*)$ and our data distribution $\mu^*$ via the Kullback-Leibler divergence KL. To this end, we consider the following assumptions on the parameterized score and the original data distribution:

**Assumption 2.1.** There exists $\epsilon > 0$ such that

$$ \max_{0 \le k \le M} \mathbb{E} \left[ \sum_{\ell=1}^{d} (1 - s_{T_f - t_k}^{\theta^*, \ell}) h \left( \frac{1 - s_{T_f - t_k}^{\ell}}{1 - s_{T_f - t_k}^{\theta^*, \ell}} \right) (\overline{\hat{X}}_{t_k}) \right] \le \epsilon, \qquad (23) $$

with $h(a) := a \log(a) - (a - 1)$ for $a > 0$.

Note that Assumption 2.1 is induced by the entropic term $\mathcal{L}_e$ defined in (16) of the loss function we consider in practice. This condition naturally appears as we bound the KL divergence of the path probability measures corresponding to the approximate score $s^{\theta^*}$ and the ideal one $s$ respectively. Indeed, we prove a Girsanov type theorem which provide an explicit expression of the density between these two measures in Theorem F.3 in the supplement F.1.1. While standard Girsanov
theorem for diffusion implies an L²-type approximation error condition for generative models (see e.g., Conforti et al. (2024); Lee et al. (2023); Chen et al. (2022a)), our result naturally involve the entropic-type condition (23) due to the discrete structure of our noising process.

**Assumption 2.2.** The data distribution does not admit any zero-value, i.e., $\mu^*(x) \in (0, 1)$ for any $x \in X$.

Assumption 2.2 follows that the data distribution has the finite Fisher-like information

$$ \beta_{\gamma^d}(\mu^*) := \mathbb{E} \left[ \sum_{\ell=1}^{d} h \left( e^{g(\bar{X}_0) - g(\varphi^{(\ell)}(\bar{X}_0))} \right) \right] < +\infty, \quad (24) $$

with $g := -\log(d\mu^*/d\gamma^d)$.

Note that Assumption 2.2 is put in parallel with the finite relative Fisher information condition provided by Conforti et al. (2024). However, Assumption 2.2 is much simpler as the state space considered is finite, and the function $h$ is only infinite if $\mu^*$ has not full support.

We are now ready to state the error's bound of the generated data using DMPMs given in Algorithm 3.

**Theorem 2.3.** Under Assumption 2.1 and Assumption 2.2, the following bound holds

$$ \mathrm{KL}(\mu^*|\mathrm{Law}(\bar{X}_{T_f}^*)) \le e^{-T_f} \mathrm{KL}(\mu^*|\gamma^d) + \tau \beta_{\gamma^d}(\mu^*) + \epsilon T_f, \quad (25) $$

with $\tau := \max\{h_k, k = 1, \dots, K\}$.

Theorem 2.3 is one of our distinguishing results, which guarantees the convergence of DMPMs algorithm, and makes it stronger than other algorithms built before for discrete target distribution.

The term $\epsilon T_f$ in (25) appears because the score function $s_t$ is replaced in the discretization by its approximation $s_t^{\theta^*}$ satisfying Assumption 2.1. The term $e^{-T_f} \mathrm{KL}(\mu^*|\gamma^d)$ represents the initialization error, as our backward dynamic starts at $\gamma^d$ instead of $\mu_{T_f}$. Finally, the term $\tau \beta_{\gamma^d}(\mu^*)$ means that the data distribution $\mu^*$ cannot be peculiar, in the sense that $\mu^*$ does not admit any zero-value. The detailed proof of Theorem 2.3 is given in the supplementary material F.3. We deduce the following complexity result for DMPMs to achieve an $\epsilon > 0$ discretization error.

**Corollary 2.4.** Consider fixed step-size $h_k = h$ for any $h > 0$. Assume Assumption 2.1 and Assumption 2.2 hold. If for $\epsilon > 0$, we set the step size and the number of iterations as

$$ h \le \frac{\epsilon}{2\beta_{\gamma^d}(\mu^*)} \\ K_f \ge \frac{\log(2\mathrm{KL}(\mu^*|\gamma^d)/\epsilon)}{h}, $$

then setting the horizon $T_f = hK_f$, it holds $\mathrm{KL}(\mu^*|\mathrm{Law}(\bar{X}_{T_f}^*)) \le \epsilon + \epsilon T_f$.

In our next result, we get rid of Assumption 2.2 using an early stopping strategy.

**Theorem 2.5.** Under Assumption 2.1, for any $\eta \in (0, T_f)$, consider the discrete time scheme $\{\bar{t}_m\}_{m=0}^M$ associated with step-sizes $\{\bar{h}_m\}_{m=1}^M$, $\bar{t}_m = \sum_{i=1}^m \bar{h}_i$ such that $\bar{t}_0 = 0$ and $\bar{t}_M = T_f - \eta$. Then, the following bound holds

$$ \mathrm{KL}(\mu_{\eta}|\mathrm{Law}(\bar{X}_{T_f-\eta}^*)) \le e^{-T_f} \mathrm{KL}(\mu^*|\gamma^d) + \bar{\tau} \beta_{\gamma^d}(\mu_{\eta}) + \epsilon(T_f - \eta), \quad (26) $$

with $\bar{\tau} := \max\{\bar{h}_m, m = 1, \dots, M\}$.
The proof of Theorem 2.5 is a consequence of Theorem 2.3. Note that (10) yields that $\mu_{\eta}$ is always positive for any $\eta \in (0, T_f)$, thus the Fisher-like information $\beta_{\gamma^d}(\mu_{\eta})$ is always finite and Assumption 2.2 holds.

To obtain then a complexity bound for DMPMs on its discretization error without Assumption 2.2, we bound in our next result, the total variation distance between $\mu^*$ and $\mu_{\eta}$ for $\eta > 0$.

**Proposition 2.6.** For any $\eta \in (0, T_f)$, we have

$$ \|\mu_{\eta} - \mu^*\|_{\text{TV}} \le 2 - 2 \left( \frac{1}{2} + \frac{1}{2} e^{-2\lambda\eta} \right)^d \le 2 - 2(1 - \lambda\eta)^d. \quad (27) $$

The proof of Proposition 2.6 is provided in the supplement F.4. Combining Theorem 2.5 and Proposition 2.6, we deduce

**Corollary 2.7.** Consider fixed step-size $\bar{h}_k = \bar{h}$ for any $\bar{h} > 0$. Assume Assumption 2.1 holds. If for $\varepsilon > 0$ sufficiently small, we set the step size and the number of iterations as,

$$ \begin{align} \eta &\le \frac{1 - (1 - \varepsilon/2)^{1/d}}{\lambda} \nonumber \\ \bar{h} &\le \frac{\varepsilon^2 (\lambda\eta)^d}{2^{d+3} d (1 + 2\lambda\eta)^d} \tag{28} \\ \bar{K}_f &\ge \frac{\log(2\text{KL}(\mu^*|\gamma^d)/\varepsilon^2) - \eta}{h}, \nonumber \end{align} $$

then setting the horizon $T_f - \eta = \bar{h}\bar{K}_f$, it holds

$$ \|\mu^* - \text{Law}(\overleftarrow{X}_{T_f - \eta}^*)\|_{\text{TV}} \le 2\varepsilon + \sqrt{2\varepsilon T_f}. $$

The detailed proof of Corollary 2.7 is given in Appendix F.4.

# 3 Existing works on diffusion-based generative models for discrete data

We provide details of existing approaches in discrete generative models. For further details see Appendix A.

A variety of diffusion-based techniques have been adapted for discrete data, often by embedding categorical variables into continuous spaces (Dieleman et al., 2022; Chen et al., 2022b; Richemond et al., 2022). Although this preserves many strengths of continuous diffusion, it can require heavy training regimes and lacks certain theoretical guarantees. Other methods, such as Argmax Flows and Multinomial Diffusion (Hoogeboom et al., 2021), use categorical noise models or argmax transformations to handle discrete tokens, but can impose considerable computational overhead.

Recently, continuous-time Markov chains have become central to modeling discrete diffusion. Campbell et al. (2022) pioneered a CTMC-based approach, on top of which Gat et al. (2024) adapted flow matching to the discrete setting, with correction steps for further performance gains at the expense of sampling cost, while Holderrieth et al. (2024) proposed a more general generator-matching framework that adapts to arbitrary Markov processes at the cost of potential complexity. Masked
diffusion models (Austin et al., 2021; Shi et al., 2024) and stochastic integral formulations (Ren et al., 2024) further exemplify modern attempts at balancing tractability, performance, and theoretical underpinnings. However, many of these methods still lack rigorous error bounds or scale poorly in high dimensions.

Our proposed method takes a step toward bridging these gaps, yielding provable error guarantees and reduced computational burdens, expressing the score function as a conditional expectation and, e.g., avoiding the costly signal-to-noise ratio training used in Shi et al. (2024).

# 4 Experiments

The full experimental details are available in Appendix D. We evaluate our Discrete Markov Path Model (DMPM) on two datasets. The first is a low-dimensional synthetic sawtooth dataset, with dimension $4 \le d \le 16$. The second is binarized MNIST, with $d = 32 \times 32$. We explore various design choices, and compare DMPM against MD4 (masked diffusion) (Shi et al., 2024) and DFM (discrete flow matching) (Gat et al., 2024), two state-of-the-art discrete generative approaches.

## 4.1 Experiments on Small-Dimensional Bernoulli Data

We study a discrete data distribution $p$ such that each component of $X = (X_i)_{i=1}^d \sim p$ is independently distributed as Bernoulli($p_i$). The map $i \mapsto p_i$ forms a sawtooth pattern (see Figure 6). We evaluate performance using a custom Sliced Wasserstein Distance (SWD) between the learned and true distributions (see Appendix D.3). Indeed, the state space size $2^d$ can get too big for traditional histogram-based metrics like KL divergence or Hellinger distance.

<figure><img src="image_2.png" /><figcaption>Figure 1: Comparison of time-schedules (cosine, linear, quadratic) and time horizon ($T_f = 3$ vs. $T_f = 10$).</figcaption></figure>

**Time horizon and time-schedule** We vary the time horizon $T_f$ and consider different time-schedules: uniform, quadratic, cosine (see Table 2) and investigate impact on performance. In Figure 1, we show results for a model trained with the simplest $\mathcal{L}_{\text{L}^2}$ loss with $d = 16$, and evaluated with various reverse steps. The cosine schedule with $T_f = 3$ achieves optimal SWD values,
<table><caption>Figure 2: SWD ↓, in 1e-3, for DMPM, MD4, and DFM across data dimension *d*. Selected the best result with #steps $2 \le K \le 200$ for each method.</caption><thead><tr><th>d</th><th>4</th><th>8</th><th>12</th><th>16</th></tr></thead><tbody><tr><th>DFM</th><td>6.102</td><td>8.864</td><td>5.019</td><td>8.302</td></tr><tr><th>MD4</th><td>9.376</td><td>7.670</td><td>4.045</td><td>8.037</td></tr><tr><th>DMPM</th><td><strong>3.174</strong></td><td><strong>3.308</strong></td><td><strong>2.342</strong></td><td><strong>2.515</strong></td></tr></tbody></table>

outperforming linear/quadratic schedules and longer horizons, and uses less reverse steps. This suggests $T_f = 3$ sufficiently approaches the uniform distribution during forward diffusion, avoiding excessive uniform-state transitions. In light of these observations, we fix this choice in subsequent experiments.

**Comparison with state of the art methods** We compare DMPM (cosine schedule, $T_f = 3$, loss $\mathcal{L}_{L^2}$) against MD4 and DFM. Figure 2 reports SWD with varying data dimension *d*. We find that DMPM outperforms both baselines, with significantly fewer reverse steps required for optimal performance (typically 30 vs 100).

## 4.2 Experiments on Higher-Dimensional Binary MNIST

**DMPM sampler and flip-schedule** We investigate our DMPM sampler with a constant and linear flip-schedule $\{M_{t_k}\}_{k=1}^K$. In practice, we observe optimal performance when we set $\sum_{k=1}^K M_{t_k} = 1000$, thus encouraging as many total bit flips as the data dimension *d*. We dimension each flip-schedule accordingly. Figure 3 illustrate the performance of our two schedule strategies. The performance is stable across the exact choice of reverse steps *K*, as long as total bit flips $\sum_{k=1}^K M_{t_k} = 1000$ stays constant. This facilitates large speedups by reducing reverse steps, thus network calls, for similar performance. Using a model trained with the loss $\mathcal{L}_{1/3,1/3,1/3}^w$, and using the linear flip-schedule, we achieve an optimal FID of 4.77 with only 25 network calls. The linear flip-schedule consistently outperforms the constant one, as fewer bits flipped in the early phases help the model converge to more sensible samples, working similarly to the cosine masking schedule introduced in MD4 (Shi et al., 2024).

**Configuration for the loss function** Through extensive experimentation, the balanced loss $\mathcal{L}_{1,1,1}^w$ yields better performance. The scale factor *w* notably helps to balance magnitudes at each timestep, and improves the synergy of the *l*₂, cross-entropy, and KL components. We refer the reader to Figure 5 for an appropriate illustration of these comments. It should be noted that the simplest losses $\mathcal{L}_{L^2}$, $\mathcal{L}_{L^2}^w$ already yield excellent, close to optimal results.

**Denoise-renoise sampler and comparison with state-of-the-art.** We further exploit the discrete-denoiser structure with our denoise-renoise sampler, as given in Algorithm 5, Appendix C.4. This approach tends to leverage the model’s learned transitions effectively, leading to further improvements in sample quality. We compare our DMPM model, trained with the balanced $\mathcal{L}_{1,1,1}^w$ loss, with the denoise-renoise sampler and the default sampler with a linear flip schedule, to MD4 (masked diffusion) and DFM (discrete flow matching).
<figure><img src="image_3.png" /><figcaption>Figure 3: FID↓ on MNIST, linear vs. constant flip-schedules scaled for 1000 total bit flips, with various loss configurations.</figcaption></figure>

<table><caption>Table 1: FID↓ (first row of each method) and $F_1^{dc}$ ↑ (second row) on MNIST for various total reverse steps. We highlight the best result in **bold**, the 2nd best in *italics*, and *underline* the 3rd best.</caption><thead><tr><th>Method</th><th></th><th>10</th><th>25</th><th>50</th><th>100</th><th>200</th><th>500</th></tr></thead><tbody><tr><td rowspan="2">DFM</td><td>FID</td><td>227.55</td><td>156.26</td><td>88.93</td><td>39.62</td><td>16.26</td><td>7.34</td></tr><tr><td>F<sup>dc</sup><sub>1</sub></td><td>0.00</td><td>0.00</td><td>0.01</td><td>0.14</td><td>0.41</td><td>0.68</td></tr><tr><td rowspan="2">MD4</td><td>FID</td><td>97.97</td><td>33.50</td><td>14.06</td><td>6.83</td><td>4.48</td><td><u>3.43</u></td></tr><tr><td>F<sup>dc</sup><sub>1</sub></td><td>0.04</td><td>0.29</td><td>0.57</td><td>0.76</td><td>0.83</td><td>0.86</td></tr><tr><td rowspan="2">DMPM<sub>flip</sub></td><td>FID</td><td>16.30</td><td>9.98</td><td>11.07</td><td>9.07</td><td>7.80</td><td>10.84</td></tr><tr><td>F<sup>dc</sup><sub>1</sub></td><td>0.64</td><td>0.92</td><td>0.93</td><td>0.93</td><td>0.93</td><td>0.70</td></tr><tr><td rowspan="2">DMPM<sub>denoise</sub></td><td>FID</td><td>78.20</td><td>20.94</td><td>8.62</td><td><u>3.98</u></td><td><b>2.89</b></td><td>4.36</td></tr><tr><td>F<sup>dc</sup><sub>1</sub></td><td>0.13</td><td>0.67</td><td>0.87</td><td><u>0.96</u></td><td><u>1.00</u></td><td><b>1.00</b></td></tr></tbody></table>

For each method, we vary the total number of reverse steps K. We report both the Fréchet Inception Distance (FID) and an $F_1^{dc}$ score, which is a harmonic mean of coverage and density metrics (Naeem et al., 2020). These two metrics complement each other, with FID measuring the global realism of generated samples and $F_1^{dc}$ capturing how well the generated data distribution covers the real data (coverage) while maintaining sample fidelity (density). In other words, $F_1^{dc}$ provides a more localized, distributional perspective; see Appendix D.3 for further details.

As shown in Table 1, the proposed DMPM approaches (rows 3 and 4) consistently outperform the baselines (DFM and MD4) across a range of step counts. At $K = 200$ reverse steps, DMPM (denoise-renoise) achieves the lowest FID of 2.89 (compared to 4.48 for MD4 and 16.26 for DFM), alongside an $F_1^{dc}$ of 1.00. Even at a lower number of steps (e.g., $K = 50$), DMPM (denoise-renoise) attains an FID of 8.62, while preserving a strong $F_1^{dc}$ of 0.87. Meanwhile, DMPM (flip-schedule) shows a similarly favorable trade-off, achieving FID below 10 for $K = 25$ and $F_1^{dc}$ above 0.90, which are remarkable results for this few network calls.

To visually validate generation quality, we refer the reader to Figure 7 for an image grid generated using the DMPM sampler, with 25 reverse steps, and to Figure 8 for an image grid generated using the denoise-renoise sampler, with 200 reverse steps.
## 4.3 Conclusions

Our experiments demonstrate that DMPM matches or outperforms state-of-the-art discrete generative models, on low and high dimensional data. On binarized MNIST, DMPM obtains lower FID and $F_1^{dc}$ than baseline methods, like Discrete Flow Matching, with at least 5× fewer network calls. We trace this success to the discrete diffusion structure, and subsequent sensible design choices. We hope for fast progress and adoption of our DMPM models, thanks to their similarity with existing continuous diffusion approaches.

## Acknowledgements

The work of A. Ocello was funded by the European Union (ERC-2022-SYG-OCEAN-101071601). Views and opinions expressed are however those of the author only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them.

## References

Jacob Austin, Daniel D Johnson, Jonathan Ho, Daniel Tarlow, and Rianne Van Den Berg. Structured denoising diffusion models in discrete state-spaces. *Advances in Neural Information Processing Systems*, 34:17981–17993, 2021.

Dominique Bakry, Ivan Gentil, Michel Ledoux, et al. *Analysis and geometry of Markov diffusion operators*, volume 103. Springer, 2014.

Omer Bar-Tal, Hila Chefer, Omer Tov, Charles Herrmann, Roni Paiss, Shiran Zada, Ariel Ephrat, Junhwa Hur, Guanghui Liu, Amit Raj, et al. Lumiere: A space-time diffusion model for video generation. In *SIGGRAPH Asia 2024 Conference Papers*, pages 1–11, 2024.

Andrea Bertazzi, Dario Shariatian, Umut Simsekli, Eric Moulines, and Alain Durmus. Piecewise deterministic generative models, 2024. URL https://arxiv.org/abs/2407.19448.

Andrew Campbell, Joe Benton, Valentin De Bortoli, Thomas Rainforth, George Deligiannidis, and Arnaud Doucet. A continuous time framework for discrete denoising models. *Advances in Neural Information Processing Systems*, 35:28266–28279, 2022.

Nanxin Chen, Yu Zhang, Heiga Zen, Ron J Weiss, Mohammad Norouzi, and William Chan. Wavegrad: Estimating gradients for waveform generation. *arXiv preprint arXiv:2009.00713*, 2020.

Sitan Chen, Sinho Chewi, Jerry Li, Yuanzhi Li, Adil Salim, and Anru R Zhang. Sampling is as easy as learning the score: theory for diffusion models with minimal data assumptions. *arXiv preprint arXiv:2209.11215*, 2022a.

Ting Chen, Ruixiang Zhang, and Geoffrey Hinton. Analog bits: Generating discrete data using diffusion models with self-conditioning. *arXiv preprint arXiv:2208.04202*, 2022b.

Giovanni Conforti and Christian Léonard. Time reversal of markov processes with jumps under a finite entropy condition. *Stochastic Processes and their Applications*, 144:85–124, 2022.
Giovanni Conforti, Alain Durmus, and Marta Gentiloni Silveri. Kl convergence guarantees for score diffusion models under minimal data assumptions, 2024. URL https://arxiv.org/abs/2308.12240.

Sander Dieleman, Laurent Sartran, Arman Roshannai, Nikolay Savinov, Yaroslav Ganin, Pierre H Richemond, Arnaud Doucet, Robin Strudel, Chris Dyer, Conor Durkan, et al. Continuous diffusion for categorical data. arXiv preprint arXiv:2211.15089, 2022.

Itai Gat, Tal Remez, Neta Shaul, Felix Kreuk, Ricky TQ Chen, Gabriel Synnaeve, Yossi Adi, and Yaron Lipman. Discrete flow matching. arXiv preprint arXiv:2407.15595, 2024.

Ivan Gavrilyuk, Volodymyr Makarov, and Vitalii Vasylyk. Exponentially convergent algorithms for abstract differential equations. Springer Science & Business Media, 2011.

Jonathan Ho, Tim Salimans, Alexey Gritsenko, William Chan, Mohammad Norouzi, and David J Fleet. Video diffusion models. Advances in Neural Information Processing Systems, 35:8633-8646, 2022.

Peter Holderrieth, Marton Havasi, Jason Yim, Neta Shaul, Itai Gat, Tommi Jaakkola, Brian Karrer, Ricky TQ Chen, and Yaron Lipman. Generator matching: Generative modeling with arbitrary markov processes. arXiv preprint arXiv:2410.20587, 2024.

Emiel Hoogeboom, Didrik Nielsen, Priyank Jaini, Patrick Forré, and Max Welling. Argmax flows and multinomial diffusion: Learning categorical distributions. Advances in Neural Information Processing Systems, 34:12454-12465, 2021.

Tero Karras, Miika Aittala, Timo Aila, and Samuli Laine. Elucidating the design space of diffusion-based generative models, 2022. URL https://arxiv.org/abs/2206.00364.

Zhifeng Kong, Wei Ping, Jiaji Huang, Kexin Zhao, and Bryan Catanzaro. Diffwave: A versatile diffusion model for audio synthesis. arXiv preprint arXiv:2009.09761, 2020.

Tuomas Kynkäänniemi, Tero Karras, Samuli Laine, Jaakko Lehtinen, and Timo Aila. Improved precision and recall metric for assessing generative models, 2019. URL https://arxiv.org/abs/1904.06991.

Dawid Laszuk. Python implementation of empirical mode decomposition algorithm. https://github.com/laszukdawid/PyEMD, 2017.

Holden Lee, Jianfeng Lu, and Yixin Tan. Convergence of score-based generative modeling for general data distributions. In International Conference on Algorithmic Learning Theory, pages 946-985. PMLR, 2023.

Christian Léonard. Girsanov theory under a finite entropy condition. In Séminaire de Probabilités XLIV, pages 429-465. Springer, 2012.

Muhammad Ferjad Naeem, Seong Joon Oh, Youngjung Uh, Yunjey Choi, and Jaejun Yoo. Reliable fidelity and diversity metrics for generative models, 2020. URL https://arxiv.org/abs/2002.09797.

Alexander Quinn Nichol and Prafulla Dhariwal. Improved denoising diffusion probabilistic models. In International conference on machine learning, pages 8162-8171. PMLR, 2021.
Marcel Nutz. Introduction to entropic optimal transport. Lecture notes, Columbia University, 2021.

Art B. Owen. Monte Carlo theory, methods and examples. https://artowen.su.domains/mc/, 2013.

Aditya Ramesh, Prafulla Dhariwal, Alex Nichol, Casey Chu, and Mark Chen. Hierarchical text-conditional image generation with clip latents. arXiv preprint arXiv:2204.06125, 1(2):3, 2022.

Yinuo Ren, Haoxuan Chen, Grant M Rotskoff, and Lexing Ying. How discrete and continuous diffusion meet: Comprehensive analysis of discrete diffusion models via a stochastic integral framework. arXiv preprint arXiv:2410.03601, 2024.

Pierre H Richmond, Sander Dieleman, and Arnaud Doucet. Categorical sdes with simplex diffusion. arXiv preprint arXiv:2210.14784, 2022.

Robin Rombach, Andreas Blattmann, Dominik Lorenz, Patrick Esser, and Björn Ommer. High-resolution image synthesis with latent diffusion models. In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pages 10684–10695, 2022.

Chitwan Saharia, William Chan, Saurabh Saxena, Lala Li, Jay Whang, Emily L Denton, Kamyar Ghasemipour, Raphael Gontijo Lopes, Burcu Karagol Ayan, Tim Salimans, et al. Photorealistic text-to-image diffusion models with deep language understanding. Advances in neural information processing systems, 35:36479–36494, 2022.

Dario Shariatian, Umut Simsekli, and Alain Durmus. Denoising lévy probabilistic models, 2024. URL https://arxiv.org/abs/2407.18609.

Jiaxin Shi, Kehang Han, Zhe Wang, Arnaud Doucet, and Michalis K Titsias. Simplified and generalized masked diffusion for discrete data. arXiv preprint arXiv:2406.04329, 2024.

Yang Song, Jascha Sohl-Dickstein, Diederik P. Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score-based generative modeling through stochastic differential equations, 2021. URL https://arxiv.org/abs/2011.13456.

Yang Song, Prafulla Dhariwal, Mark Chen, and Ilya Sutskever. Consistency models, 2023. URL https://arxiv.org/abs/2303.01469.

Nizar Touzi. Optimal stochastic control, stochastic target problems, and backward SDE, volume 29. Springer Science & Business Media, 2012.

Ruben Villegas, Mohammad Babaeizadeh, Pieter-Jan Kindermans, Hernan Moraldo, Han Zhang, Mohammad Taghi Saffar, Santiago Castro, Julius Kunze, and Dumitru Erhan. Phenaki: Variable length video generation from open domain textual descriptions. In International Conference on Learning Representations, 2022.

EUN BI YOON, Keehun Park, Sungwoong Kim, and Sungbin Lim. Score-based generative models with lévy processes. In A. Oh, T. Naumann, A. Globerson, K. Saenko, M. Hardt, and S. Levine, editors, Advances in Neural Information Processing Systems, volume 36, pages 40694–40707. Curran Associates, Inc., 2023. URL https://proceedings.neurips.cc/paper_files/paper/2023/file/8011b23e1dc3f57e1b6211ccad498919-Paper-Conference.pdf.

G Zitkovic. Uniform integrability. course notes: Theorem if Probability II, Fall 2015: Lecture, 12, 2015.
# A Existing works on diffusion-based generative models for discrete data

This section provides details of the recent researches on discrete generative models.

**Embedding discrete structure in the continuous space.** To keep the benefits of continuous representations, Dieleman et al. (2022) and Chen et al. (2022b) mapped discrete structures into Euclidean space, while Richemond et al. (2022) placed them into the simplex, all while continuing to use forward continuous diffusion models. In particular, Dieleman et al. (2022) proposed a continuous diffusion model for categorical data, which has some advantages over autoregressive models, such as the ability to perform arbitrary infilling and a more flexible sampling process. However, this method comes with an expensive training cost and lacks of strong theoretical guarantees.

**Argmax flows and Multinomial Diffusion.** Hoogeboom et al. (2021) introduced two new generative models, Argmax Flows and Multinomial Diffusion, to handle categorical data like text and image segmentation. Argmax Flows connect discrete data with continuous models by using an argmax function combined with a probabilistic inverse, making categorical distributions easy-learning. Multinomial Diffusion process uses a categorical distribution to add noise to discrete data and then trains a model to reverse the process. However, both Argmax Flows and Multinomial Diffusion have some limitations: computational costs increase due to additional steps, and the theoretical guarantee is missing.

**Designing the flow processes over the discrete state space.** Campbell et al. (2022) introduced the first complete continuous-time framework for denoising diffusion models applied to discrete data. They used CTMCs to model the forward noising process and its time-reversal dynamics. While the core idea is similar to ours, their approach is more complex because their method consider generic CTMC and is not specialized to the noising process that we consider. As a result, their method essentially boils down learning density ratios which can be computationally demanding and fail to offer efficient approximation in high dimensions. They also added a correction step to bring the sample distribution closer to the desired one, which increased the practical training cost Gat et al. (2024). By focusing on the random-walk CTMC on X, we were able to provide a discrete counterpart to the score function that is learn in continuous diffusion models and also to establish strong convergence guarantees for our method.

**Generator Matching.** Another recent approach to handle discrete data is generative modeling with arbitrary Markov processes using generator matching, introduced by Holderrieth et al. (2024). In this approach, the authors design an appropriate Markov process that transforms a simple distribution into the desired data one using a generator, which can be efficiently trained with a neural network. This method is quite flexible and can be applied to different state spaces, especially in discrete settings. However, this method being very generic suffer from the same drawback as Campbell et al. (2022).

**Masked diffusion models.** One important step toward more advanced models is the "masked" diffusion process, a discrete diffusion approach first introduced by Austin et al. (2021). Recently, Shi et al. (2024) looked into this model further, simplifying its training objective by expressing it as a signal-to-noise ratio, which helps highlight some useful features. However, despite these improvements, the model still lacks theoretical guarantees and remains expensive to train in practice.

**Discrete Diffusion Models via a Stochastic Integral Framework.** Ren et al. (2024) introduced a new way to analyze discrete diffusion models via Lévy-type stochastic integrals and expanded Poisson random measures. Specifically, they established the stochastic integral expressions of the noising and denoising processes for the categorical data. They provided a unified error analysis
framework and showed the first error bound for their algorithms in KL divergence. However, their results rely on strong assumptions in contrast to our results. Besides, our bounds are simpler and better, in particular with respect to the time horizon.

Our paper takes a step toward bridging these gaps. By clearly describing the forward Markov process, we can express the score function as a conditional expectation, which helps us avoid the costly signal-to-noise ratio training used in Shi et al. (2024). This way, we not only offer a simpler and more affordable training approach, but also provide solid theoretical guarantees for our models in practice.

# B Interpretation of DMPMs

## B.1 The simple case $X = \{0, 1\}$

Detailed calculation of the transition probability in (2). Based on the Kolmogorov equation, the transition matrix $\vec{p}_t$ for $0 \le t \le T_f$ admits the following formula

$$ \vec{p}_t^1 = e^{t} \vec{q}, $$

where $\vec{q}$ is defined in (1). Clearly, the generator $\vec{q}$ admits two eigenvalues 0 and $-2\lambda$ associated with the eigenvectors $(1 \quad 1)^T$ and $(1 \quad -1)^T$ respectively. Then we can diagonalize $\vec{q}$ as

$$ \vec{q} = \begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix} \begin{pmatrix} 0 & 0 \\ 0 & -2\lambda \end{pmatrix} \begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix}^{-1}, $$

and the transition matrix $\vec{p}_t^1$ follows

$$ \vec{p}_t^1 = e^{t} \vec{q} = \begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix} \begin{pmatrix} 1 & 0 \\ 0 & e^{-2\lambda} \end{pmatrix} \begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix}^{-1} = \frac{1}{2} \begin{pmatrix} 1 + e^{-2\lambda} & 1 - e^{-2\lambda} \\ 1 - e^{-2\lambda} & 1 + e^{-2\lambda} \end{pmatrix}. $$

□

## B.2 General state space $X = \{0, 1\}^d$

### B.2.1 Forward transition probability

*Proof of (9).* We start with a note that the generator matrix $\vec{q}$ can be expressed as a sum of matrices $\vec{q}^\ell$ as follows

$$ \vec{q} = \sum_{\ell=1}^{d} \vec{q}^{\ell}, \quad \text{with} \quad \vec{q}^{\ell}(x, y) = \begin{cases} \lambda, & \text{if } x^i = y^i \text{ for } i \neq \ell \text{ and } x^{\ell} \neq y^{\ell}, \\ -\lambda, & \text{if } x = y, \\ 0, & \text{otherwise.} \end{cases} $$

Notice that $\vec{q}^\ell$ also admits the following formula with respect concerning the tensor product

$$ \vec{q}^{\ell} = \underbrace{\mathbb{I} \otimes \mathbb{I} \otimes \dots \otimes \mathbb{I} \otimes \vec{A} \otimes \mathbb{I} \otimes \dots \otimes \mathbb{I}}_{d \text{ times}}, \qquad (29) $$
with $\mathbb{I}$ the 2 × 2 identity matrix and $\vec{A} = \begin{pmatrix} -\lambda & \lambda \\ \lambda & -\lambda \end{pmatrix}$, which is the $\ell^{th}$ matrix in the previous product. Indeed, by the definition of tensor product, for any $x = (x^i)_{i=1}^d, y = (y^i)_{i=1}^d \in \mathcal{X}$, we observe that

$$ (\mathbb{I} \otimes \mathbb{I} \otimes \dots \otimes \mathbb{I} \otimes \underbrace{\vec{A}}_{\ell^{th}} \otimes \mathbb{I} \otimes \dots \otimes \mathbb{I})(x, y) = \mathbb{I}(x^1, y^1)\mathbb{I}(x^2, y^2)\dots\vec{A}(x^\ell, y^\ell)\dots\mathbb{I}(x^d, y^d) \\ = \begin{cases} 1, & \text{if } x^i = y^i \text{ for } i \neq \ell \text{ and } x^\ell \neq y^\ell, \\ -1, & \text{if } x = y, \\ 0, & \text{otherwise.} \end{cases} $$

which is exactly the expression of $\vec{q}^\ell(x, y)$. We now use the Kolmogorov equation combined with the expression of $\vec{q}_t^\ell$ in (29), and apply the formula $e^{\mathbb{I}\otimes A+B\otimes\mathbb{I}} = e^A \otimes e^B$ for any matrix A, B (Gavrilyuk et al., 2011, Appendix) to get

$$ \vec{p}_t = e^{t\vec{q}} = e^{\sum_{\ell=1}^{d} t \vec{q}^{\ell}} = \underbrace{e^{t\vec{A}} \otimes \dots \otimes e^{t\vec{A}}}_{d \text{ times}}. $$

We are thus left with the computation of $e^{t\vec{A}}$. It is clear that the eigenvalues of $\vec{A}$ are 0 and $-2\lambda$, with the corresponding eigenvectors $(1 \quad 1)^T$ and $(1 \quad -1)^T$ respectively. Consequently, we can compute $e^{t\vec{A}}$ as: for any $a, b \in \{0, 1\}$,

$$ \vec{p}_t^1(a, b) := e^{t\vec{A}}(a, b) = \begin{cases} \frac{1}{2} + \frac{1}{2}e^{-2t}, & \text{if } a = b, \\ \frac{1}{2} - \frac{1}{2}e^{-2t}, & \text{if } a \neq b, \end{cases} $$

and the formula of transition probability $\vec{p}_t$ for $0 \le t \le T_f$ follows: for any $x = (x^i)_{i=1}^d$ and $y = (y^i)_{i=1}^d$ in $\mathcal{X}$,

$$ \vec{p}_t(x, y) = \prod_{i=1}^{d} \vec{p}_t^1(x^i, y^i), \quad \text{with} \quad \vec{p}_t^1(x^i, y^i) = \begin{cases} \frac{1}{2} + \frac{1}{2}e^{-2\lambda t}, & \text{if } x^i = y^i, \\ \frac{1}{2} - \frac{1}{2}e^{-2\lambda t}, & \text{otherwise.} \end{cases} $$

□

### B.2.2 Conditional expectation expression of the score function

**Proof of Proposition 1.1.** Fix $x \in \mathcal{X}$ and $\ell = 1, \dots, d$. By the formula of the marginal distribution, we see that

$$ \mu_{T_f-t}(x) - \mu_{T_f-t}(\varphi^{(\ell)}(x)) = \sum_{z \in \mathcal{X}} \mu_0(z) (\vec{p}_{T_f-t}(z, x) - \vec{p}_{T_f-t}(z, \varphi^{(\ell)}(x))). \quad (30) $$
The formula of transition probabilities $\vec{p}_{T_f-t}(z, \varphi^{(\ell)}(x))$ combined with the definition of $\varphi^{(\ell)}(x)$ lead to

$$ 
\begin{aligned}
\vec{p}_{T_f-t}(z, \varphi^{(\ell)}(x)) &= \prod_{i=1}^{d} \vec{p}_{T_f-t}^{1}(z^i, \varphi_{\ell}^{i}(x)) \\
&= \vec{p}_{T_f-t}^{1}(z^{\ell}, \varphi_{\ell}^{\ell}(x)) \prod_{\substack{i=1 \\ i \neq \ell}}^{d} \vec{p}_{T_f-t}^{1}(z^i, x^i) \\
&= \frac{\vec{p}_{T_f-t}^{1}(z^{\ell}, \varphi_{\ell}^{\ell}(x))}{\vec{p}_{T_f-t}^{1}(z^{\ell}, x^{\ell})} \vec{p}_{T_f-t}(z, x).
\end{aligned}
$$

Substituting this into (30) implies

$$ 
\begin{aligned}
\mu_{T_f-t}(x) - \mu_{T_f-t}(\varphi^{(\ell)}(x)) &= \sum_{z \in \mathcal{X}} \mu_0(z) \vec{p}_{T_f-t}^{1}(z, x) \left(1 - \frac{\vec{p}_{T_f-t}^{1}(z^{\ell}, \varphi^{(\ell),\ell}(x))}{\vec{p}_{T_f-t}^{1}(z^{\ell}, x^{\ell})}\right) \\
&= \sum_{z \in \mathcal{X}} \left[ \frac{2e^{-2\lambda(T_f-t)}}{1+e^{-2\lambda(T_f-t)}} - \frac{4e^{-2\lambda(T_f-t)}(x^{\ell} - z^{\ell})^2}{1-e^{-4\lambda(T_f-t)}} \right] \mathbb{P}[\vec{X}_0 = z, \vec{X}_{T_f-t} = x],
\end{aligned}
$$

where the last equality comes from the formula of $\vec{p}_{T_f-t}^{1}$ and the fact that if $z^{\ell} = \varphi^{(\ell),\ell}(x)$ then $z^{\ell} \neq x^{\ell}$. Therefore, the score function in components are

$$ 
\begin{aligned}
s_t^{\ell}(x) &= \frac{\mu_{T_f-t}(x) - \mu_{T_f-t}(\varphi^{(\ell)}(x))}{\mu_{T_f-t}(x)} \\
&= \sum_{z \in \mathcal{X}} \left[ \frac{2e^{-2\lambda(T_f-t)}}{1+e^{-2\lambda(T_f-t)}} - \frac{4e^{-2\lambda(T_f-t)}(x^{\ell} - z^{\ell})^2}{1-e^{-4\lambda(T_f-t)}} \right] \mathbb{P}[\vec{X}_0 = z | \vec{X}_{T_f-t} = x] \\
&= \mathbb{E} \left[ \frac{2\alpha_{T_f-t}}{1+\alpha_{T_f-t}} - \frac{4\alpha_{T_f-t}(\vec{X}_{T_f-t}^{\ell} - \vec{X}_0^{\ell})^2}{1-\alpha_{T_f-t}^2} | \vec{X}_{T_f-t} = x \right],
\end{aligned}
$$

where $\alpha_t = e^{-2\lambda t}$, and we finish the proof of Proposition 1.1. $\square$

### B.2.3 Invariant measure of the forward process

As we have a comprehensive understanding of the forward process, we observe that its invariant measure is the uniform distribution over $\mathcal{X}$, denoted by $\gamma^d$. Indeed, for any $x \in \mathcal{X}$ and $t \in [0, T_f]$,

$$ (\gamma^d \vec{p}_t)(x) = \sum_{z \in \mathcal{X}} \gamma^d(z) \vec{p}_t(z, x) = \frac{1}{2^d} \sum_{z \in \mathcal{X}} \vec{p}_t(z, x) = \frac{1}{2^d} = \gamma^d(x). $$

Furthermore, by formula of $\vec{p}$ given in (9), we have $\vec{p}_t(x, y) \xrightarrow{t \to \infty} \frac{1}{2^d}$ for any $x, y \in \mathcal{X}$. Consequently, the following holds for any $x \in \mathcal{X}$,

$$ \mu_t(x) = \sum_{z \in \mathcal{X}} \mu_0(z) \vec{p}_t(z, x) \xrightarrow{t \to \infty} \frac{1}{2^d} \sum_{z \in \mathcal{X}} \mu_0(z) = \frac{1}{2^d} = \gamma^d(x), $$

meaning that the forward dynamic $(\vec{X}_t)_{t \in [0, T_f]}$ converges geometrically fast to $\gamma^d$.
# C Implementation of DMPMs

## C.1 Alternative ideal backward simulation

Besides the simulation of the backward process provided in Section 1.2.2, we can also use the following procedure to produce the time-reversal dynamic.

The second procedure to sample $(\overleftarrow{X}_t)_{t \in [0, T_f]}$ is to consider a sample $\overleftarrow{X}_0$ from $\mu_{T_f}$ and a sequence of i.i.d. random variables distributed according to the exponential distribution with parameter 1, $\{E_i^\ell : i \in \mathbb{N}, \ell \in \{1, \dots, d\}\}$, we can define the jump times $(T_i)_{i \in \mathbb{N}}$ of the backward process and its transition by induction setting $T_0 = 0$. Given $(T_i, \overleftarrow{X}_{T_i})$, we define the next jump time as $T_{i+1}^j = T_i + \Delta T_{i+1}^j$, where $\Delta T_{i+1}^j = \inf\{t \ge 0 : \int_0^t \lambda(1 - s^j(\overleftarrow{X}_{T_i}))dr \ge E_i^j\}$. Then, set $T_{i+1} = T_{i+1}^{\ell_i}$, where $\ell_i = \arg \min_{j \in \{1, \dots, d\}} T_{i+1}^j$, and $\overleftarrow{X}_t = \overleftarrow{X}_{T_i}$ for $t \in (T_i, T_{i+1} \wedge T_f)$, and finally if $T_{i+1} < T_f$, $\overleftarrow{X}_{T_{i+1}}^{\ell_i} = 1 - \overleftarrow{X}_{T_i}^{\ell_i}$ for $\ell_i \in \{1, \dots, d\}$.

## C.2 Perfect backward approximation

We provide here the pseudo-code of backward approximation sampling in continuous time scheme:

**Algorithm 1 DMPMs Algorithm (Continuous time scheme)**

**Input:** a time horizon $T_f \gg 1$ large enough, a prescribed jump rate $\lambda$, an approximate score function $s^{\theta^*}$

**Backward process:**

Set $T_0 = 0$ and initialize $\overleftarrow{X}_0 \sim \gamma^d$
$i \leftarrow 0$

**while** $T_i \le T_f$ **do**

Draw $E_i \sim \text{Exp}(1)$

Solve $\Delta T_{i+1} = \inf\{t \ge 0 : \int_0^t \lambda_{T_i+r}^{\theta^*}(\overleftarrow{X}_{T_i})dr \ge E_i\}$, with $\lambda_t^{\theta^*}(x) = \lambda \sum_{\ell=1}^d (1 - s_t^{\theta^*, \ell}(x))$

Set $T_{i+1} = T_i + \Delta T_{i+1}$

**if** $T_i < t < \min(T_{i+1}, T_f)$ **then**

Set $\overleftarrow{X}_t = \overleftarrow{X}_{T_i}$

**if** $T_{i+1} < T_f$ **then**

Draw $\ell_i \in \{1, \dots, d\} \sim \text{Cate}(\{\lambda(1 - s_{T_{i+1}}^{\theta^*, \ell}(\overleftarrow{X}_{T_i}))/\lambda_{T_{i+1}}^{\theta^*}(\overleftarrow{X}_{T_i})\}_{\ell=1}^d)$

Set $\overleftarrow{X}_{T_{i+1}} = \varphi^{(\ell_i)}(\overleftarrow{X}_{T_i})$

$i \leftarrow i + 1$

**Output:** $\overleftarrow{X}_{T_f}$

## C.3 Discrete denoiser and score reparameterization

**Discrete-denoiser structure** Recall from Proposition 1.1 that each score component admit the following conditional expectation:

$$s_t^\ell(x) = \mathbb{E} \left[ f_t^\ell \left( \overrightarrow{X}_0^\ell, \overrightarrow{X}_{T_f-t} \right) \middle| \overrightarrow{X}_{T_f-t} = x \right], \quad (31)$$
where

$$f_t^\ell(\vec{X}_0^\ell, \vec{X}_{T_f-t}) = \frac{2\alpha_{T_f-t}}{1+\alpha_{T_f-t}} - \frac{4\alpha_{T_f-t}(\vec{X}_{T_f-t}^\ell - \vec{X}_0^\ell)^2}{1-\alpha_{T_f-t}^2} \quad (32)$$

for $t \in [0, T_f)$, $x \in \mathbf{X}$ and $\ell = 1, \dots, d$. Remark that

$$\mathbb{E}[f_t^\ell(\vec{X}_0^\ell, \vec{X}_{T_f-t})|\vec{X}_{T_f-t} = x] = \frac{2\alpha_{T_f-t}}{1+\alpha_{T_f-t}} - \frac{4\alpha_{T_f-t}\mathbb{E}[(\vec{X}_{T_f-t}^\ell - \vec{X}_0^\ell)^2|\vec{X}_{T_f-t} = x]}{1-\alpha_{T_f-t}^2}. \quad (33)$$

Thus we introduce the function $d_t^\ell$ defined as

$$d_t^\ell : x \mapsto \mathbb{E}[(\vec{X}_{T_f-t}^\ell - \vec{X}_0^\ell)^2 | \vec{X}_{T_f-t} = x], \quad (34)$$

which can be further rewritten as

$$\begin{aligned} d_t^\ell(x) &= \mathbb{E} \left[ \left( \vec{X}_{T_f-t}^\ell - \vec{X}_0^\ell \right)^2 \middle| \vec{X}_{T_f-t} = x \right] \\ &= \mathbb{E} \left[ 1_{\vec{X}_{T_f-t}^\ell \neq \vec{X}_0^\ell} \middle| \vec{X}_{T_f-t} = x \right] \\ &= \mathbb{P} \left( \vec{X}_0^\ell \neq x^\ell \middle| \vec{X}_{T_f-t} = x \right). \end{aligned}$$

In some sense, this is the discrete version of the continuous denoiser $\mathbb{E}[\vec{X}_0|\vec{X}_t]$ approximated by classical diffusion models (Song et al., 2021), as obtained from the score by Tweedie's formula. Thus we call $d_t^\ell(x)$ the discrete denoiser.

**Score reparameterization** Based on the previous derivations, each score component $s_t^\ell(x)$ can be written as a function of $d_t^\ell$:

$$s_t^\ell(x) = \frac{2\alpha_{T_f-t}}{1+\alpha_{T_f-t}} - \frac{4\alpha_{T_f-t}d_t^\ell(x)}{1-\alpha_{T_f-t}^2}, \quad (35)$$

So we can reparameterize our score models $s_t^\theta$ as

$$s_t^{\theta, \ell}(x) = \frac{2\alpha_{T_f-t}}{1+\alpha_{T_f-t}} - \frac{4\alpha_{T_f-t}d_t^{\theta, \ell}(x)}{1-\alpha_{T_f-t}^2}, \quad (36)$$

where $d_t^{\theta, \ell}(x)$ aims to approximate $d_t^\ell(x)$.

## C.4 Objective functions derived from the discrete denoiser structure

Inspired by the previous derivations, we modify our existing $\mathcal{L}_2$ loss function to replace by a denoising loss equivalent. We introduce a cross-entropy loss, and finally propose a scaling of the loss functions, based on the average output magnitude of the discrete denoiser, thus helping with the learning, and improving synergies between loss elements.
**Score-matching objective $\mathcal{L}_{L^2}$** We rewrite the objective function $\mathcal{L}_{L^2}$ to fit the discrete denoiser, considered as a conditional expectation:

$$ \mathcal{L}_{L^2}^{\text{denoiser}} : \theta \mapsto \int_{0}^{T_f} \mathcal{L}_{t,L^2}(\theta)dt, \quad \mathcal{L}_{t,L^2}(\theta) = \mathbb{E} \left[ \|d_{T_f-t}^{\theta}(\vec{X}_t) - (\vec{X}_0 - \vec{X}_t) \odot (\vec{X}_0 - \vec{X}_t)\|^2 \right], \quad (37) $$

where $\odot$ is the element-wise product.

**Cross-entropy objective $\mathcal{L}_{CE}$** Instead of the $\mathcal{L}_{L^2}$ loss suggested by the conditional expectation structure, we can consider a cross-entropy loss to fit our model to the correct distribution: classical derivations from the conditional log-likelihood $\sum_{\ell=1}^d \mathbb{E} \left[ \log p_t^{\theta,\ell}(\vec{X}_{T_f-t} | \vec{X}_0^\ell) \right]$, where

$$ p_t^{\theta,\ell}(x_{T_f-t}|x_0) = \begin{cases} d_t^{\theta,\ell}(x_{T_f-t}) & \text{if } x_{T_f-t} \neq x_0, \\ 1 - d_t^{\theta,\ell}(x_{T_f-t}) & \text{else} \end{cases}, \quad (38) $$

lead to the following cross entropy loss:

$$ \mathcal{L}_{CE}(\theta) = - \int_{0}^{T_f} \mathcal{L}_{t,CE}(\theta)dt, \quad (39) $$

where

$$ \mathcal{L}_{t,CE}(\theta) = \mathbb{E} \left[ \sum_{l=1}^{d} Y_t^l \log d_t^{\theta,l}(\vec{X}_{T_f-t}^l) + (1 - Y_t^l) \log (1 - d_t^{\theta,l}(\vec{X}_{T_f-t}^l)) \right], \quad Y_t^l = \begin{cases} 1 & \text{if } \vec{X}_0^l \neq \vec{X}_{T_f-t}^l, \\ 0 & \text{else.} \end{cases} $$

**Further improvements** To address vanishing gradient problems, we inspect the average magnitude of the loss across the dataset, at each timestep. Indeed, the average value of $d_t^l$ is

$$ w_{T_f-t} = \mathbb{E} \left[ d_t^l(\vec{X}_{T_f-t}) \right] = \mathbb{E} \left[ \mathbb{E} \left[ d_t^l(\vec{X}_{T_f-t}) \middle| \vec{X}_0 \right] \right] \quad (40) $$

$$ = \mathbb{E} \left[ \mathbb{P} \left( \vec{X}_0^l \neq \vec{X}_{T_f-t}^l \middle| \vec{X}_0 \right) \right] \quad (41) $$

$$ = \frac{1}{2} (1 - \alpha_{T_f-t}), \quad (42) $$

as given by the formulas for the transition kernels of the forward process. We can see that the value of $w_t$ is close to zero for small values of $t$, which stalls the learning process. Empirically, we find that dividing the integrand of either loss terms $\mathcal{L}_{L^2}$ or $\mathcal{L}_{CE}$ by $w_t$ yields improvements. As a result, we modify the losses to counterbalance their diminishing magnitude across timesteps: $\mathcal{L}_{t,L^2}^w = \frac{\mathcal{L}_{t,L^2}^{\text{denoiser}}}{w_t}$, $\mathcal{L}_{t,CE}^w = \frac{\mathcal{L}_{t,CE}^{\text{denoiser}}}{w_t}$, and define the associated losses

$$ \mathcal{L}_{L^2}^w(\theta) = \int_0^{T_f} \mathcal{L}_{t,L^2}^w(\theta)dt, \quad \mathcal{L}_{CE}^w(\theta) = - \int_0^{T_f} \mathcal{L}_{t,CE}^w(\theta)dt. \quad (43) $$
**Comparing $\mathcal{L}_{L^2}$, $\mathcal{L}_{CE}$, $\mathcal{L}_{L^2}^w$, $\mathcal{L}_{CE}^w$** In Figure 4, we plot the average loss per timestep, for a trained model on MNIST (following the specifications given in Appendix D.2)). It shows that, on average, the $\mathcal{L}_L$ loss effectively becomes a scaled variant of the cross-entropy objective, which is reflected in similar performance results. This corroborates our derivation that L2 acts as an effective lower bound to the log-likelihood. This also supports its relevancy with respect to the underlying structure of this generative model. It must be noted that both losses still benefit from positive synergies when used together.

Importantly, dividing by $w_t = (1 - \alpha_t)/2$ particularly helps at smaller timesteps, and keeps the loss values at the same magnitude across timesteps, enhancing training dynamics. This is illustrated in Figure 5, where scaling the losses with the w scale factor consistently yields improvements.

<figure><img src="image_4.png" /><figcaption>Figure 4: Comparison of $\mathcal{L}_{L^2}$, $\mathcal{L}_{CE}$, $\mathcal{L}_{L^2}^w$, $\mathcal{L}_{CE}^w$ average losses over timesteps. The two losses become scaled version of one another only when averaged over data, but otherwise benefit from positive synergies when mixed together.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 5: FID↓, on MNIST, for models trained with $\mathcal{L}_{\bar{\omega}}$ and $\mathcal{L}_{\bar{\omega}}^w$ losses, evaluated using 200 reverse steps with the denoise-renoise sampler. Scaling with w yields consistent improvements, with the best loss configuration $\mathcal{L}_{1/3,1/3,1/3}^w$ involving all the methodological improvements we discussed.</figcaption></figure>

**Final objective functions** We choose a linear combination of the previous loss objectives, weighted by positive coefficients $\bar{\omega}_1, \bar{\omega}_2, \bar{\omega}_3$:

$$ \mathcal{L}_{\bar{\omega}} = \bar{\omega}_1 \mathcal{L}_{L^2}^{\text{denoiser}} + \bar{\omega}_2 \mathcal{L}_{\text{e}} + \bar{\omega}_3 \mathcal{L}_{\text{CE}}, \quad (44) $$

and, if we choose their version weighted by $1/w_t$:

$$ \mathcal{L}_{\bar{\omega}}^{w} = \bar{\omega}_{1} \mathcal{L}_{L^{2}}^{w} + \bar{\omega}_{2} \mathcal{L}_{\text{e}} + \bar{\omega}_{3} \mathcal{L}_{\text{CE}}^{w}. \quad (45) $$

## C.5 Generative process and sampling procedures

Once we obtain our neural network $d_t^\theta$ approximating $d_t$, we use it to produce fresh samples that closely mimic the observed data. To do so, we first introduce a DMPM sampler based on the true reverse process. We then propose a slight modification, leveraging the distribution on indices available at each step, by flipping multiple bits instead of just one, using a flip-schedule. Finally, we
<caption>**Algorithm 2** Training Algorithm for DMPM (Reparameterized Score)</caption>

**Require:** Dataset $\mathcal{D}$ of samples $X \in \{0, 1\}^d$;
Time horizon $T_f > 0$ and rate $\lambda > 0$;
Parameterized discrete denoiser model $\{d_t^{\theta, \ell}(x) : \theta \in \Theta\}_{t, \ell, x}$;
Derived score function $s_t^{\theta} := \frac{2\alpha_{T_f-t}}{1+\alpha_{T_f-t}} - \frac{4\alpha_{T_f-t}d_t^{\theta}}{1-\alpha_{T_f-t}^2}$ (score reparameterization (18));
Define $\alpha_t$ as in (7), $f_t$ as in (13);
Loss coefficients $\bar{\omega}_1, \bar{\omega}_2, \bar{\omega}_3 \ge 0$;
1: **while** optimization has not converged **do**
2: &nbsp;&nbsp;&nbsp;&nbsp;Sample a batch $\{X_i\}_{i=1}^B$ from $\mathcal{D}$.
3: &nbsp;&nbsp;&nbsp;&nbsp;Draw $t_1, \dots, t_B \stackrel{\text{iid}}{\sim} \text{Unif}([0, T_f])$
4: &nbsp;&nbsp;&nbsp;&nbsp;**Forward sampling:** fast simulation via $p_{t|0} = (p_{t|0}^1)^{\otimes d}$
5: &nbsp;&nbsp;&nbsp;&nbsp;**for** $i = 1$ to $B$ **do**
6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\vec{X}_{i,0} \leftarrow X_i$
7: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$p_{T_f-t_i} \leftarrow (1 - \alpha_{T_f-t_i})/2$
8: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Compute $\vec{X}_{i,T_f-t_i}$ by flipping each bit of $\vec{X}_{i,0}$ independently with probability $p_{T_f-t_i}$
9: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** Scaling losses with average $d_t$ magnitude **then**
10: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$w_i \leftarrow (1 - \alpha_{T_f-t_i})/2$
11: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**else**
12: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$w_i \leftarrow 1$
13: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\mathcal{L}_{\text{L}^2}(\theta) \leftarrow \frac{1}{B} \sum_{i=1}^B \frac{1}{w_i} \|d_t^{\theta}(\vec{X}_{i,T_f-t_i}) - (\vec{X}_{i,T_f-t_i} - \vec{X}_{i,0}) \odot (\vec{X}_{i,t} - \vec{X}_{i,0})\|^2$
14: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\mathcal{L}_{\text{CE}}(\theta) \leftarrow \frac{1}{Bd} \sum_{i=1}^B \frac{1}{w_i} \sum_{l=1}^d \left( \mathbb{1}_{\vec{X}_0^\ell \ne \vec{X}_{T_f-t_i}^\ell} \log d_{t_i}^{\theta, \ell}(\vec{X}_{T_f-t_i}^\ell) + (1 - \mathbb{1}_{\vec{X}_0^\ell \ne \vec{X}_{T_f-t_i}^\ell}) \log(1 - d_{t_i}^{\theta, \ell}(\vec{X}_{T_f-t_i}^\ell)) \right)$
15: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\mathcal{L}_{\text{e}}(\theta) \leftarrow \frac{1}{B} \sum_{i=1}^B \sum_{\ell=1}^d \left( -s_{T_f-t_i}^{\theta, \ell}(\vec{X}_{t_i}) + (f_{T_f-t_i}^\ell(\vec{X}_{t_i}) - 1) \log(1 - s_{T_f-t_i}^{\theta, \ell}(\vec{X}_{t_i})) \right)$
16: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\mathcal{L}_{\omega}(\theta) \leftarrow \omega_1 \mathcal{L}_{\text{L}^2}^{\text{denoiser}} + \omega_2 \mathcal{L}_{\text{e}} + \omega_3 \mathcal{L}_{\text{CE}}$
17: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Perform a gradient step on $\mathcal{L}_{\omega}(\theta)$ w.r.t. $\theta$.
18: &nbsp;&nbsp;&nbsp;&nbsp;**Return** the final parameter $\theta^*$.

derive a denoise-renoise sampler, solely based on the discrete denoiser structure of the problem, as inspired by similar lines of work in conitnuous diffusion.

**DMPM sampler** A first sampling procedure is given in Algorithm 3. It is designed to be as close as possible to the true backward process, while enabling efficient parallelization when implemented. It consists in a piecewise-approximation of the functions of interest, parameterized by the choice of a time discretization grid $0 = t_0 < t_1 < \dots < t_K = T_f$, which we call a **time-schedule**.

In Table 2, we give the different time-schedules we experiment with. We draw inspiration from numerous lines of work on continuous and discrete diffusion (Shi et al., 2024; Karras et al., 2022), in which these are common choices.

**DMPM sampler with flip-schedule** In Algorithm 4, we further take advantage of the specific structure of our backward process, by leveraging the distribution over indices given by the learned score model at each timestep $t$. Instead of flipping a single bit per timestep $t_k$, we flip a total of $M_{t_k}$
bits sampled without replacements from the given distribution. We call the sequence {$M_t$}$_{0\le t\le T_f}$ the flip-schedule. When a time-schedule {$t_k$}$_{k=1}^K$ has been chosen, we also call the corresponding discrete sequence {$M_{t_k}$}$_{k=1}^K$ a flip-schedule.

In Table 3, we give the two flip-schedules we explore in this paper. The choice for the linear schedule is inspired from the philosophy of the masking schedule introduced in the context of masked diffusion by Shi et al. (2024).

<table><caption>Table 2: Different time schedules ($t_k$)$_{k=1}^K$ used in our experiments. $T_f$ denotes the final time, and $K$ is the number of reverse steps.</caption><thead><tr><td>Time-schedule</td><td>Value of t<sub>k</sub></td></tr></thead><tbody><tr><td>Linear</td><td>T<sub>f</sub> (K/2)<sup>k</sup></td></tr><tr><td>Quadratic</td><td>T<sub>f</sub> (K/2)<sup>2</sup></td></tr><tr><td>Cosine</td><td>T<sub>f</sub> cos((1-k/K)π/2)</td></tr></tbody></table>

<table><caption>Table 3: Different flip schedules ($M_t$)$_{0\le t\le T_f}$ used in our experiments. In both schedules, $M$ is a constant to be fixed and controls the total number of bits flipped during generation.</caption><thead><tr><td>Flip-schedule</td><td>Value of M<sub>t</sub></td></tr></thead><tbody><tr><td>Constant</td><td>M</td></tr><tr><td>Linear</td><td>M t/T<sub>f</sub></td></tr></tbody></table>

**Denoise-renoise sampler** In Algorithm 5, we introduce the following denoise/renoise cycle, interpreting the model output $d_t^\theta$ as the probability that each bit should be flipped at timestep $t$ to reach timestep 0. After doing a full denoise pass (from time $T_f \to 0$), we noise the sample with the transition kernel of the forward process (from time $0 \to T_f - \Delta$). Then we can do another denoise pass from $(T_f - \Delta) \to 0$, etc.
<caption>**Algorithm 3 Backward sampling of DMPM with piecewise-constant score**</caption>

**Require:** Time horizon $T_f > 0$ and rate $\lambda > 0$;
$K > 0$ number of reverse steps and time-schedule $0 = t_0 < t_1 < \dots < t_K = T_f$;
Flip-schedule, i.e., sequence of positive integers $\{M_{t_k}\}_{k=1}^K$;
Discrete denoiser model $d^\theta$;
Derived score function $s_t^\theta := \frac{2\alpha_{T_f-t}}{1+\alpha_{T_f-t}} - \frac{4\alpha_{T_f-t}d_t^\theta}{1-\alpha_{T_f-t}^2}$ (score reparameterization (18));
Define $\alpha_t$ as in (7);
1: $\overleftarrow{X}_0^\theta \sim \text{Unif}(0, 1)^{\otimes d}$
2: $E \sim \mathcal{E}(1)$
3: $\Lambda \leftarrow 0$
4: **for** $k = 0$ to $K-1$ **do**
5: &nbsp;&nbsp;&nbsp;&nbsp; $\bar{\lambda}_{t_k} \leftarrow \lambda \sum_{l=1}^d (1 - s_{t_k}^{\theta,l})$
6: &nbsp;&nbsp;&nbsp;&nbsp; $\Delta t_k \leftarrow t_{k+1} - t_k$
7: &nbsp;&nbsp;&nbsp;&nbsp; $\Lambda \leftarrow \Lambda + \bar{\lambda}_{t_k} \Delta t_k$
8: **if** $\Lambda > E$ **then**
9: &nbsp;&nbsp;&nbsp;&nbsp; $\ell^* \sim \text{Cate}] \left( \left\{ \frac{\lambda (1 - s_{t_k}^{\theta,l})}{\bar{\lambda}_{t_k}} \right\}_{l=1}^d \right)$
10: &nbsp;&nbsp;&nbsp;&nbsp; $\overleftarrow{X}_{t_k}^{\theta,l^*} \leftarrow 1 - \overleftarrow{X}_{t_k}^{\theta,l^*}$
11: &nbsp;&nbsp;&nbsp;&nbsp; $\Lambda \leftarrow 0$
12: &nbsp;&nbsp;&nbsp;&nbsp; $E \sim \mathcal{E}(1)$
13: &nbsp;&nbsp;&nbsp;&nbsp; $\overleftarrow{X}_{t_{k+1}}^\theta \leftarrow \overleftarrow{X}_{t_k}^\theta$
**Output:** $\overleftarrow{X}_{T_f}^\theta$
**Algorithm 4 Backward sampling of DMPM with piecewise-constant score and flip-schedule**

**Require:** Time horizon $T_f > 0$ and rate $\lambda > 0$;
$K > 0$ number of reverse steps and time-schedule $0 = t_0 < t_1 < \cdots < t_K = T_f$;
Flip-schedule, i.e., sequence of positive integers $\{M_{t_k}\}_{k=0}^K$;
Discrete denoiser model $d^\theta$;
Derived score function $s_t^\theta := \frac{2\alpha_{T_f-t}}{1+\alpha_{T_f-t}} - \frac{4\alpha_{T_f-t}d_t^\theta}{1-\alpha_{T_f-t}^2}$ (score reparameterization (18));
Define $\alpha_t$ as in (7);
1: $\overleftarrow{X}_0^\theta \sim \text{Unif}(0, 1)^{\otimes d}$
2: $E \sim \mathcal{E}(1)$
3: $\Lambda \leftarrow 0$
4: **for** $k = 0$ to $K-1$ **do**
5: &nbsp;&nbsp;&nbsp;&nbsp; $\bar{\lambda}_{t_k} \leftarrow \lambda \sum_{l=1}^d (1 - s_{t_k}^{\theta,l})$
6: &nbsp;&nbsp;&nbsp;&nbsp; $\Delta t_k \leftarrow t_{k+1} - t_k$
7: &nbsp;&nbsp;&nbsp;&nbsp; $\Lambda \leftarrow \Lambda + \bar{\lambda}_{t_k} \Delta t_k$
8: &nbsp;&nbsp;&nbsp;&nbsp; **if** $\Lambda > E$ **then**
9: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $[\ell_1^*, \dots, \ell_M^*] \sim \text{Hypergeometric} \left( \left\{ \frac{\lambda (1 - s_{t_k}^{\theta,l})}{\bar{\lambda}_{t_k}} \right\}_{l=1}^d, M_{t_k} \right)$
10: &nbsp;&nbsp;&nbsp;&nbsp; **for** $i = 1$ to $M_{t_k}$ **do**
11: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\overleftarrow{X}_{t_k}^{\theta,l_i^*} \leftarrow 1 - \overleftarrow{X}_{t_k}^{\theta,l_i^*}$
12: &nbsp;&nbsp;&nbsp;&nbsp; $\Lambda \leftarrow 0$
13: &nbsp;&nbsp;&nbsp;&nbsp; $E \sim \mathcal{E}(1)$
14: &nbsp;&nbsp;&nbsp;&nbsp; $\overleftarrow{X}_{t_{k+1}}^\theta \leftarrow \overleftarrow{X}_{t_k}^\theta$
**Output:** $\overleftarrow{X}_{T_f}^\theta$

**Algorithm 5 Denoise–Noise Cycling with a Discrete Denoiser Model**

**Require:** Time horizon $T_f > 0$ and rate $\lambda > 0$;
$K > 0$ number of reverse steps and time-schedule $0 = t_0 < t_1 < \cdots < t_K = T_f$;
Discrete denoiser model $d^\theta$;
1: $\overleftarrow{X}_0^\theta \sim \text{Unif}(0, 1)^{\otimes d}$ {initial sample in $\{0, 1\}^d$}
2: **for** $k = 0$ to $K-1$ **do**
3: **Denoise phase:**
4: &nbsp;&nbsp;&nbsp;&nbsp; $d_{t_k} \leftarrow d_{t_k}^\theta (\overleftarrow{X}_{t_k}^\theta)$
5: &nbsp;&nbsp;&nbsp;&nbsp; Compute $\overleftarrow{X}_{T_f}^\theta$ by flipping each component $l$ of $\overleftarrow{X}_{t_k}^\theta$ with probability $d_{t_k}^l$
6: **Noise phase:**
7: &nbsp;&nbsp;&nbsp;&nbsp; Sample $\overleftarrow{X}_{t_{k+1}}^\theta \sim p_{T_f-t_{k+1}|0}(\cdot|\overleftarrow{X}_{T_f}^\theta)$, as in Algorithm 2
**Output:** $\overleftarrow{X}_{T_f}^\theta$
# D Experiments

All experiments are conducted using PyTorch. All the training and experiments are conducted on four NVIDIA RTX8000 GPU.

We use the score parameterization introduced in (36):

$$s_t^{\theta, \ell}(x) = \frac{2\alpha_{T_f-t}}{1 + \alpha_{T_f-t}} - \frac{4\alpha_{T_f-t}d_t^{\theta, \ell}(x)}{1 - \alpha_{T_f-t}^2}, \quad (46)$$

where the neural network $d_t^{\theta, \ell}(x)$ aims to approximate $d_t^{\ell}(x) = \mathbb{P}(\vec{X}_0^{\ell} \neq x^{\ell} | \vec{X}_{T_f-t} = x)$. Since the output of the neural network is $d_t^{\theta}(x) \in (0, 1)^d$, we add a sigmoid activation function at the last layer.

We consider various loss configurations $\mathcal{L}_{\omega}$, $\mathcal{L}_{\omega}^w$ as introduced in (19), (45), with 6 choices of coefficients $(\omega_1, \omega_2, \omega_3)$ normalized in the 2-simplex $\Delta_2 \subset \mathbb{R}^3$. We test all $2^3 - 1 = 7$ possible non-empty combinations, minus the single $\mathcal{L}_e$ loss combination $(\omega_2 = 1)$, as the latter only acts as entropic regularization and does not perform well by itself. This lets us study the synergies between the different loss terms.

## D.1 Small dimension data

We first conduct experiments on a discrete data distribution $p$ supported on $\{0, 1\}^d$. Each component of $X = (X_i)_{i=1}^d \sim p$ is independently distributed as a Bernoulli distribution with parameter $p_i$:

$$p(x) = \prod_{i=1}^{d} p_i(x_i), \quad (47)$$

where the map $i \mapsto p_i$ forms a sawtooth-like pattern, oscillating linearly between 0.05 and 0.95, as can be seen in Figure 6.

<figure><img src="image_6.png" /><figcaption>Figure 6: Sawtooth pattern used to define $i \mapsto p_i$, plotted with $d = 16$. Values oscillate linearly from 0.05 to 0.95, and back.</figcaption></figure>

For training, we use 20 000 datapoints resampled at each epoch, and a batch size of 1024. We train each model for 300 epochs, using AdamW with a learning rate of 1e-3. We employ a network composed of multiple MLP blocks: 4 residual blocks, each consisting of two feed-forward layers of width 256; layer normalization and SiLU activations in each block; a feed-forward embedding for the timesteps, mapping $\mathbb{R}$ to a hidden dimension of 256, whose output is then injected into each residual block by an additional MLP of dimension $256 \times 256$.
For evaluation, we estimate each distribution with 20,000 samples, and draw 1000 vectors uniformly on the simplex $\Delta_d$ to compute our SWD metric (see Appendix D.3).

## D.2 Image data

We work on the binarized MNIST dataset, which we scale from 28 × 28 to 32 × 32 in order to fit in the U-Net architecture. We set the pixel value to 0 if its intensity is below 0.5, and to 1 otherwise.

We compare DMPM to MD4 (masked diffusion, as in Shi et al. (2024)) and DFM (discrete flow matching, as in Gat et al. (2024)). We reimplement MD4 with the cosine schedule and the algorithms given in Appendix F of Shi et al. (2024). We implement DFM based on the Pytorch implementation in https://github.com/gle-bellier/discrete-fm, and we use corrector sampling for better results.

For DMPM, we are using the cosine time-schedule and time horizon $T_f = 3$. For both MD4 and DFM, we set the mask value to the integer 2.

To establish a fair comparison, we use the same network model for every method. We use a U-Net following the implementation of Nichol and Dhariwal (2021) available in https://github.com/openai/improved-diffusion. We dimension the network as follows.

The first layer is an embedding layer of output dimension 32 and input dimension $d_{input}$, where $d_{input} = 2$ for DMPM (input values are either 0 and 1) and $d_{input} = 3$ for MD4 and DFM (input values are either 0, 1 or the mask value 2).

We set the hidden layers to [128, 256, 256, 256], fix the number of residual blocks to 2 at each level, and add self-attention block at resolution $16 \times 16$, using 4 heads. We use an exponential moving average with a rate of 0.99. We use the silu activation function at every layer. Timestep $t$ is fed to the model through the Transformer sinusoidal position embedding.

For DMPM and MD4 we set the number of output channels to 1 and add a sigmoid activation at the last layer. For DFM, we set the output channels to 3 and apply softmax channel-wise.

The optimizer is AdamW with learning rate 5e-4. We use the StepLR scheduler which scales the learning rate by $\gamma = .99$ every 400 steps. We train on MNIST for 120 000 steps with batch size 256. A single training run on MNIST takes approximately 6 hours per GPU, and requires about 6-12GB of VRAM for our settings.

To assess the quality of our generative models, we compute our metrics between 4 000 real images and 4 000 generated images. Generating 4 000 images with 1 000 reverse steps takes approximately 2 hours on one GPU.

## D.3 Metrics

For low-dimensional data, we use a custom sliced Wasserstein metric. For image data, in addition to the classical FID metric, we use a $F_1^{DC}$ summary score, based on the density and coverage metrics.

**$F_1^{DC}$ as summary metric of density-coverage** The density and coverage metrics are introduced in the setting of generative models by Naeem et al. (2020). They assess the overlap of sample distributions using local geometric structures. Density measures how much the generated distribution is contained in the original data distribution (measuring quality), and coverage measures how much of the original data distribution is covered by the generated distribution (diversity).

These metrics are improvements of the precision and recall metrics for generative models (Kynkäänniemi et al., 2019). They offer different measures to characterize the performance of
generative models. For instance they can decorrelate the negative effect of mode collapse from the negative effect of noisy/blurry generations, each of them decreasing respectively coverage and density, and have been of importance in recent studies, e.g., in heavy-tailed generative modeling (Shariatian et al., 2024; YOON et al., 2023).

We consider a single summary $F_1^{DC}$ score, which we define as the harmonic mean of these two values:

$$F_1^{DC} = 2 \cdot \frac{\text{density} \cdot \text{coverage}}{\text{density} + \text{coverage}}. \qquad (48)$$

**Sliced Wasserstein metric SWD** Since the state space of our dataset over $\{0, 1\}^d$ is of size $2^d$, we cannot work with histogram-based metrics, which would require exponentially many samples when $d$ increases.

We address this issue with our sliced Wasserstein metric SWD. This metric is defined between distributions $\mu, \nu$ on $\{0, 1\}^d$ as:

$$\text{SWD}(\mu, \nu) = \int_{\Delta_d} W(u_{\#}\mu, u_{\#}\nu) \, du, \qquad (49)$$

where, for $u \in \Delta_d$, the pushforward $u_{\#}$ is derived from the function

$$x \in \{0, 1\}^d \mapsto \langle u, x \rangle \in [0, 1]. \qquad (50)$$

Simple Monte-Carlo averages are used to evaluate the integral with respect to the uniform distribution over the simplex $\Delta_d$, and we compute the Wasserstein distance between the pushforward measures with the `pyemd` package (Laszuk, 2017).

# E Additional results

In this section, we give grid images of generated samples for DMPM models trained on binarized MNIST, with the loss $\mathcal{L}_{1/3,1/3,1/3,1/3}^w$.
<figure><img src="image_7.png" /><figcaption>Figure 7: Default DMPM sampler, 25 reverse steps, cosine time-schedule, linear flip-schedule dimensioned for 1000 total bit flips.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 8: Denoise-renoise sampler, 200 reverse steps, cosine time-schedule.</figcaption></figure>

# F Convergence of DMPMs

The proof of DMPMs' convergence requires understanding the backward dynamic under the canonical process point of view equivalent to the transition matrix point of view we provided in Section 1.

## F.1 Canonical process point of view

We now want to give a description of the time reversal process as the solution of an optimal control process like in the continuous setting in Conforti et al. (2024). To this purpose, we consider the following canonical setting. Let $\Omega = D([0, T_f]; \mathcal{X})$ be the canonical space of all càdlàg (right continuous and left limited) paths from $[0, T_f]$ to $\mathcal{X}$. With abuse of notation, we denote as $(\vec{X}_t)_{t \in [0, T_f]}$ the canonical process defined by

$$ \vec{X}_t(\omega) = \omega_t, \quad \text{for } t \in [0, T_f], (\omega_s)_{s \in [0, T_f]} \in \Omega. $$

Denote by $\mu_t = \text{Law}(X_t)$ the law of $\vec{X}_t$, $\Delta\vec{X}_t = \vec{X}_t - \vec{X}_{t-}$ the size of jump at time $t$ for $t \in [0, T_f]$. Consider one-hot vectors $(e_i)_{i=1}^d$. Denote by $Q_d$ the set of all effective jumps, i.e.

$$ Q_d = \{(-1)^p e_i \quad \text{for } i = 1, \dots, d \text{ and } p = 0, 1\}. $$

Denote $|q| := (|q^j|)_{j=1}^d$ for $q \in Q_d$. Endow $\Omega$ with the $\sigma$-field $\sigma(X_t; t \in [0, T_f])$ generated by the canonical projections and consider the associated canonical filtration. We denote the set

$$ A := \{(x, q) \in \mathcal{X} \times Q_d \text{ such that } (x + q) \in \mathcal{X}\}, $$
and the rate of jump $\epsilon_t$ by

$$ \epsilon_t(\vec{X}_{t-}, \vec{X}_{t-} + q) := \begin{cases} 1, & \text{if } (\vec{X}_{t-}, q) \in \mathcal{A}, \\ 0, & \text{else.} \end{cases} $$

We then denote the kernel $\bar{L}$ of the canonical process $(\vec{X}_t)_{t \in [0, T_f]}$ as

$$ \begin{aligned} \bar{L}_{\omega}(dqdt) &:= dt \sum_{i=1}^{d} (1_{\{\vec{X}_{t-}^i=0\}} \epsilon_t(\vec{X}_{t-}, \vec{X}_{t-} + q)\delta_{\{e_i\}}(dq) + 1_{\{\vec{X}_{t-}^i=1\}} \epsilon_t(\vec{X}_{t-}, \vec{X}_{t-} + q)\delta_{\{-e_i\}}(dq)) \\ &= dt \sum_{i=1}^{d} (1_{\{\vec{X}_{t-}^i=0\}} \delta_{\{e_i\}}(dq) + 1_{\{\vec{X}_{t-}^i=1\}} \delta_{\{-e_i\}}(dq)) \\ &=: L_{\omega}(t, dq)dt. \end{aligned} $$

We now follow the approach for time reversal used in Léonard (2012) to characterize the evolution of $X$ with the use of the following martingale problem.

**Definition F.1 (Martingale problem).** We say that $\mathbb{P} \in \mathcal{P}(\Omega)$ solves the martingale problem $MP(\bar{L})$, and write $\mathbb{P} \in MP(\bar{L})$, if the integrability assumption

$$ \mathbb{E}_{\mathbb{P}} \int_{[0, T_f]} \sum_{q \in \mathcal{Q}_d} (|q|^2 \wedge 1) \bar{L}(dqdt) < \infty $$

holds, and if the process

$$ \sum_{x \in \mathcal{X}} f(x) \mu_t(x) - \sum_{x \in \mathcal{X}} f(x) \mu_0(x) - \int_0^t \sum_{x \in \mathcal{X}} \sum_{q \in \mathcal{Q}_d} (f(x+q) - f(x)) \mu_{s-}(x) \bar{L}(dsdq) $$

is a local $\mathbb{P}$-martingale for any function $f: \mathcal{X} \to \mathbb{R}$, where $\mu_t$ is the marginal distribution of $\vec{X}_t$.

Note that $\mathbb{P} \in MP(\bar{L})$ is equivalent to $\mathbb{P} \in \mathcal{P}(\Omega)$ such that

$$ \mathbb{E}_{\mathbb{P}} \sum_{0 \le s \le t: \Delta \vec{X}_s \ne 0} g(s, \Delta \vec{X}_s) = \mathbb{E}_{\mathbb{P}} \int_{[0, T_f]} \sum_{q \in \mathcal{Q}_d} g(s, q) \bar{L}(dqds), $$

for all measurable functions $g$.

**Definition F.2 (Condition (U)).** Léonard (2012) One says that $\mathbb{P} \in MP(\bar{L})$ satisfies the uniqueness condition (U) if for any probability measure $\mathbb{P}'$ on $\Omega$ such that the initial laws $\mathbb{P}'_0 = \mathbb{P}_0$ are equal, $\mathbb{P}' \ll \mathbb{P}$ and $\mathbb{P}' \in MP(\bar{L})$, we have $\mathbb{P} = \mathbb{P}'$.

For the sake of this paper, let $\vec{R}$ be an invariant probability measure on $\Omega$ that fulfills the uniqueness condition (U), and solves the martingale problem $MP(\bar{L})$ w.r.t. the canonical process $(\vec{X}_t)_{t \in [0, T_f]}$.

Following (Léonard, 2012, Corollary 2.7), the process $(\vec{X}_t)_{t \in [0, T_f]}$ can be decomposed as

$$ \vec{X} = \vec{X}_0 + q \odot \mu_L^X, \quad \vec{R} - \text{a.s.}, $$

with $\mu_L^X := \sum_{t \in [0, T_f]: \Delta \vec{X}_t \ne 0} \delta_{(t, \Delta \vec{X}_t)}$.
### F.1.1 Girsanov's theorem

From Léonard (2012) we see that, for jump processes, the relative entropy of two path measures can be decomposed with the help of the function $\varrho(a) := e^a - a - 1$, for $a \in \mathbb{R}$, and its convex conjugate $\varrho^*(b) = (b+1)\log(b+1) - b$ for $b > -1$ with convention $\varrho^*(-1) = 1$ and $\varrho^*(b) = \infty$ for $b < -1$. This is proven by the following theorem.

**Theorem F.3 (Girsanov's theorem).** Let $\mathbb{P} \in \mathcal{P}(\Omega)$ verifying $\text{KL}(\mathbb{P}|\vec{\mathbb{R}}) < \infty$. Then, there exists a unique predictable non-negative process $u: \Omega \times [0, T_f] \times \mathcal{Q}_d \to [0, \infty)$ satisfying the integrability condition

$$ \mathbb{E}_{\mathbb{P}} \int_{[0, T_f]} \sum_{q \in \mathcal{Q}_d} \varrho^*(|u-1|) d\bar{L} < \infty, $$

and $\mathbb{P} \in \text{MP}(u\bar{L})$. Moreover, we have that

$$ \text{KL}(\mathbb{P}|\vec{\mathbb{R}}) = \text{KL}(\mathbb{P}_0|\vec{\mathbb{R}}_0) + \mathbb{E}_{\mathbb{P}} \int_{[0, T_f]} \sum_{q \in \mathcal{Q}_d} h(u_t(\vec{X}_t, q)) dt, $$

$$ \text{with } h(a) := \varrho^*(a-1) = a \log a - a + 1. $$

The proof of Theorem F.3 relies on several following technical lemmas. Let us introduce their framework first. Let $\mathbb{P} \in \mathcal{P}(\Omega)$ such that $\mathbb{P} \in \text{MP}(\bar{K})$, and $\chi$ a $\mathbb{R}$-valued predictable process on $\Omega \times [0, T_f] \times \mathcal{Q}_d$ such that $\int \sum_{\mathcal{Q}_d} \varrho(\chi) d\bar{K} < \infty$. We define

$$ Z_t := \exp \left( \chi \odot \tilde{\mu}_t^K - \int_{[0, T_f]} \sum_{\mathcal{Q}_d} \varrho(\chi) d\bar{K} \right), \quad \text{for } t \in [0, T_f], $$

and the stopping time for $k, j \ge 1$,

$$ \sigma_j^k := \inf \left\{ t \in [0, T_f] ; \int_{[0, T_f]} \sum_{q \in \mathcal{Q}_d} \varrho(\chi) d\bar{K} \ge k \text{ or } (\Delta \vec{X}_t \ne 0 \text{ and } \rho_t(\Delta \vec{X}_t) \notin [-j, k]) \right\}. $$

**Lemma F.4.** Assume further that $\chi$ satisfies

$$ \mathbb{E}_{\mathbb{P}} \int_{[0, T_f]} \sum_{q \in \mathcal{Q}_d} \varrho(\rho_t(q)) \bar{K}(dq dt) < \infty. \quad (51) $$

Then $\chi \odot \tilde{\mu}^K$ is a $\mathbb{P}$-martingale with $\tilde{\mu}^K := \mu_K^X - \bar{K} = \sum_{t \in [0, T_f]: \Delta \vec{X}_t \ne 0} \delta_{(t, \Delta \vec{X}_t)} - \bar{K}$. Moreover, the process $Z$ defined as above is a local $\mathbb{P}$-martingale and a positive $\mathbb{P}$-supermartingale, which satisfies

$$ dZ_t = Z_{t-} \left[ (e^{\chi(q)} - 1) \odot d\tilde{\mu}_t^K \right]. $$

*Proof of Lemma F.4.* This result is an adaptation of Lemma 6.1 in Léonard (2012). From its definition, we have that $\tilde{\mu}^K$ is a $\mathbb{P}$-martingale measure. Therefore the stochastic integral

$$ M_t^\chi := \chi \odot \tilde{\mu}_t^K = \int_{[0, T_f]} \sum_{q \in \mathcal{Q}_d} \lambda_s(q) d\tilde{\mu}^K $$
is a local $\mathbb{P}$-martingale. Denote $Y_t := M_t^\rho - \int_{[0,T_f]} \beta_s ds$ with $\beta_t := \sum_{q \in \mathcal{Q}_d} \varrho(\rho_t(q)) K(t,q)$. Recalling Itô's formula for the jump process $(Y_t)_{t \in [0,T_f]}$, we have

$$ df(Y_t) = \left[ \sum_{q \in \mathcal{Q}_d} [f(Y_{t-} + \chi(q)) - f(Y_{t-})] K(t, q) \right] dt + \nabla f(Y_{t-}) \cdot \beta_t dt + dM_t, \quad \mathbb{P}\text{-a.s.}, $$

for $f$ a measurable function, with $M$ a local $\mathbb{P}$-martingale. Using this formula for $f(y) = e^y$, we obtain

$$ \begin{aligned} \mathrm{d}e^{Y_t} &= \left[ \sum_{q \in \mathcal{Q}_d} (e^{Y_{t-}+\chi(q)} - e^{Y_{t-}}) K(t,q) \right] \mathrm{d}t - e^{Y_{t-}} \beta_t \mathrm{d}t + \mathrm{d}M_t \\ &= e^{Y_{t-}} \beta_t \mathrm{d}t - e^{Y_{t-}} \beta_t \mathrm{d}t + \mathrm{d}M_t = \mathrm{d}M_t, \end{aligned} $$

$\mathbb{P}$-a.s., with $M$ a local $\mathbb{P}$-martingale. This implies $Z = e^Y$ is a local $\mathbb{P}$-martingale and, since $Z$ is positive, we can conclude that $Z$ is a $\mathbb{P}$-supermartingale. Moreover, we have

$$ \begin{aligned} \mathrm{d}e^{Y_t} &= e^{Y_t-} [\varrho(\Delta Y_t) + \mathrm{d}Y_t] \\ &= e^{Y_t-} \left[ \varrho(\chi(q)) \odot \mathrm{d}\tilde{\mu}_t^K + \left( \sum_{q \in \mathcal{Q}_d} \varrho(\rho_t(q)) K(t,q) \right) \mathrm{d}t - \beta_t \mathrm{d}t + \chi(q) \odot \mathrm{d}\tilde{\mu}_t^K \right] \\ &= e^{Y_t-} \left[ \varrho(\chi(q)) \odot \mathrm{d}\tilde{\mu}_t^K + \chi(q) \odot \mathrm{d}\tilde{\mu}_t^K \right] \\ &= e^{Y_t-} \left[ (e^{\chi(q)} - 1) \odot \mathrm{d}\tilde{\mu}_t^K \right], \end{aligned} $$

which finishes the proof of Lemma F.4. $\square$

**Lemma F.5.** For all $j, k \ge 1$, $Z_j^{\sigma_j^k}$ is a genuine $\mathbb{P}$-martingale and the measure

$$ Q_j^k := Z_1^{\sigma_j^k} \mathbb{P}_j^k $$

is a probability measure on $\Omega$ which satisfies

$$ Q^k \in \mathrm{MP}(\mathbb{1}_{[0,\sigma_j^k]} e^{\chi} \bar{K}). $$

*Proof of Lemma F.5.* Fix $j, k \ge 1$. We have

$$ Z_j^{\sigma_j^k} = \exp(\chi_j^k \odot \tilde{\mu}^K - \int_{[0,T_f]} \sum_{q \in \mathcal{Q}_d} \varrho(\chi_j^k) d\bar{K}), $$

where $\chi_j^k = \mathbb{1}_{[0,\sigma_j^k]} \chi$ is predictable since $\chi$ is predictable and $\mathbb{1}_{[0,\sigma_j^k]}$ is left continuous, and $\chi_j^k \odot \tilde{\mu}^K := \int_{[0,T_f]} \sum_{q \in \mathcal{Q}_d} \rho_s(q) d\tilde{\mu}^K$. For simplicity, we write $\chi = \mathbb{1}_{\{\lambda_j^k \in [-j,k]\}} \chi_j^k$ and $Z = Z_j^{\sigma_j^k}$ in the rest of the proof. From the definition of $\sigma_j^k$, we obtain

$$ \int_{[0,T_f]} \sum_{q \in \mathcal{Q}_d} \varrho(\chi) d\bar{K} \le k, \quad \text{for } \chi \in [-j, k], \quad \mathbb{P}_j^k \text{-a.s.} \qquad (52) $$
First, we prove that $Z$ is a $\mathbb{P}_j^k$-martingale. From Lemma F.4, $Z$ is a local martingale, it is enough to show that

$$ \mathbb{E}_{\mathbb{P}_j^k} Z_T^p < \infty \quad \text{for some } p > 1. $$

For all $p \ge 0$, we have

$$ Z^p = \exp \left( p\chi \odot \tilde{\mu}^K - p \int_{[0, T_f]} \sum_{q \in \mathcal{Q}_d} \varrho(\chi) d\bar{K} \right) \le \exp \left( p\chi \odot \tilde{\mu}^K \right), $$

and

$$ \exp \left( p\chi \odot \tilde{\mu}^K - \int_{[0, T_f]} \sum_{q \in \mathcal{Q}_d} \varrho(p\chi) d\bar{K} \right) \ge \exp \left( p\chi \odot \tilde{\mu}^K \right) / C(k, p), $$

for some finite deterministic constant $C(k, p) > 0$ since $\varrho(p\chi) \le c(k, p)\varrho(\chi)$ holds for all $\chi \le k$ and some constant $0 < c(k, p) < \infty$. Then we have

$$ \begin{aligned} \exp \left( \int_{[0, T_f]} \sum_{q \in \mathcal{Q}_d} \varrho(p\chi) d\bar{K} \right) &\le \exp \left( \int_{[0, T_f]} \sum_{q \in \mathcal{Q}_d} c(k, p) \varrho(\chi) d\bar{K} \right) \\ &\le \exp(kc(k, p)) =: C(k, p). \end{aligned} $$

This implies

$$ Z^p \le \exp(p\chi \odot \tilde{\mu}^K) \le C(k, p) \exp \left( p\chi \odot \tilde{\mu}^K - \int_{[0, T_f]} \sum_{q \in \mathcal{Q}_d} \varrho(p\chi) d\bar{K} \right). $$

Applying Lemma F.4 for $p\chi$, we deduce that $\exp(p\chi \odot \tilde{\mu}_t^K - \int_{[0, T_f]} \sum_{q \in \mathcal{Q}_d} \varrho(p\chi) d\bar{K})$ is a $\mathbb{P}_j^k$-supermartingale. This yields

$$ \mathbb{E}_{\mathbb{P}_j^k} \exp \left( p\chi \odot \tilde{\mu}^K - \int_{[0, T_f]} \sum_{q \in \mathcal{Q}_d} \varrho(p\chi) d\bar{K} \right) \le \exp \left( p\chi \odot \tilde{\mu}_0^K - \int_{[0,0]} \sum_{q \in \mathcal{Q}_d} \varrho(p\chi) d\bar{K} \right) = 1. $$

Therefore,

$$ Z^p \le C(k, p) < \infty, $$

which allow us to conclude that $Z$ is a $\mathbb{P}_j^k$-martingale (see, e.g, Zitkovic, 2015). Then since $\mathbb{E}_{\mathbb{P}_j^k}(Z_T) = Z_0 = 1$, we have $1 = \int Z_1 \mathbb{P}_j^k = \int Q_j^k$, i.e. $Q_j^k$ is a probability measure on $\Omega$. Now, we show that

$$ Q_j^k \in \mathrm{MP}(\mathbb{1}_{[0, \sigma_j^k]} e^{\chi} \bar{K}). $$

Let $\tau$ be a finitely valued stopping time which will be specified later, and for any measurable function $f$, we denote $g(t, \Delta \vec{X}_t) := f(\vec{X}_{t-} + \Delta \vec{X}_t) - f(\vec{X}_{t-})$ with $g(t, 0) = 0$ for all $t \in [0, T_f]$. Denote

$$ F_t := \sum_{0 \le s \le t \wedge \tau} g(s, \Delta \vec{X}_s). $$
By Lemma F.4, the martingale Z satisfies

$$dZ_t = 1_{[0, \sigma_j^k]}(t)Z_{t-} \left[ (e^\chi - 1) \odot \tilde{\mu}^K \right],$$

and

$$dF_t = 1_{[0, \tau]}(t)g(t, \Delta \vec{X}_t).$$

Therefore,

$$d[Z, F]_t = 1_{[0, \sigma_j^k \wedge \tau]}(t)Z_{t-} (e^{\chi(\Delta \vec{X}_t)} - 1)g(t, \Delta \vec{X}_t), \quad \mathbb{P}_j^k \text{ a.s.}$$

Combining these above with the fact that Z is $\mathbb{P}_j^k$-martingale, we get

$$
\begin{align*}
& \mathbb{E}_{Q_j^k} \sum_{0 \le s \le t \wedge \tau} g(s, \Delta \vec{X}_s) \\
&= \mathbb{E}_{P_j^k} (Z_{t \wedge \tau} F_{t \wedge \tau} - Z_0 F_0) \\
&= \mathbb{E}_{P_j^k} \int_{[0, t \wedge \tau]} (F_s dZ_s + Z_s dF_s + d[Z, F]_s) \\
&= \mathbb{E}_{P_j^k} \left[ \int_{[0, \tau]} F_s dZ_s + \sum_{0 \le s \le t \wedge \tau} Z_{s-} g(t, \Delta \vec{X}_s) + \sum_{0 \le s \le t \wedge \tau} Z_{s-} (e^{\chi(s, \Delta \vec{X}_s)} - 1) g(s, \Delta \vec{X}_s) \right] \\
&= \mathbb{E}_{P_j^k} \sum_{0 \le s \le t \wedge \tau} Z_{t-} e^{\chi(s, \Delta \vec{X}_s)} g(s, \Delta \vec{X}_s) \\
&= \mathbb{E}_{P_j^k} \int_{[0, t \wedge \tau]} \sum_{q \in Q_d} Z_{s-} g(s, q) e^{\chi(s, q)} \bar{K}(dq ds) \\
&= \mathbb{E}_{Q_j^k} \int_{[0, t \wedge \tau]} \sum_{q \in Q_d} g(s, q) e^{\chi(s, q)} \bar{K}(dq ds),
\end{align*}
$$

where we used $\mathbb{P} \in LK(\bar{K})$ in the last equality. Choosing $\tau$ such that the above terms are meaningful, we can conclude that $Q_j^k \in MP(e^{\chi_j^k} \bar{K})$. $\square$

**Proof of Theorem F.3.** This proof is an adaptation of Theorem 2.6 in Léonard (2012) based on technical lemmas provided above. By Lemma F.4, we have $0 < \mathbb{E}_{\bar{R}} Z_t \le 1$ for all $\chi$ satisfies integrability condition (51). Note that the relative entropy admits the following variational representation

$$\mathrm{KL}(\mathbb{P}|\bar{R}) = \sup \left\{ \int u d\mathbb{P} - \log \int e^u d\bar{R}; \ u : \int e^u d\bar{R} < \infty \right\}, \quad (53)$$

for any $\mathbb{P} \in \mathcal{P}(\Omega)$ such that $\mathrm{KL}(\mathbb{P}|\bar{R}) < \infty$. Using this expression of $\mathrm{KL}(\mathbb{P}|\bar{R})$, we have

$$\mathbb{E}_{\mathbb{P}}(\chi \odot \tilde{\mu}_{T}^{L} - \int_{[0, T_f]} \sum_{q \in Q_d} \varrho(\lambda_t(q)) \bar{L}(dq dt)) \le \mathrm{KL}(\mathbb{P}|\bar{R}),$$
for any $\chi$ satisfying (51). Therefore,

$$ \mathbb{E}_{\mathbb{P}}(\chi \odot \tilde{\mu}^L) \leq \text{KL}(\mathbb{P}|\vec{R}) + \int_{[0,T_f]} \sum_{q \in \mathcal{Q}_d} \varrho(\chi) d\bar{L}. $$

Consider $\|\cdot\|_{\varrho}$ defined as

$$ \|\chi\|_{\varrho} := \inf \left\{ a > 0; \mathbb{E}_{\mathbb{P}} \int_{[0,T_f]} \sum_{q \in \mathcal{Q}_d} \varrho(\chi/a) \bar{L}(dqdt) \le 1 \right\}. $$

This norm is the Luxemburg norm of the small Orlicz space

$$ S_{\varrho} := \left\{ \chi : \Omega \times [0, T_f] \times \mathcal{Q}_d \to \mathbb{R}; \text{ measurable s.t. } \mathbb{E}_{\mathbb{P}} \int_{[0,T_f]} \sum_{q \in \mathcal{Q}_d} \varrho(b|\chi|) d\bar{L} < \infty, \forall b \ge 0 \right\}, $$

Taking $\phi := \frac{\chi}{\|\chi\|_{\varrho}}$ implies

$$ \mathbb{E}_{\mathbb{P}}(\phi \odot \tilde{\mu}^L) \le (\text{KL}(\mathbb{P}|\vec{R}) + 1)\|\phi\|_{\varrho}, \quad \forall \phi. \quad (54) $$

Consider now the space $\mathcal{B}$ of all bounded processes such that

$$ \mathbb{E}_{\mathbb{P}} \int_{[0,T_f]} \sum_{q \in \mathcal{Q}_d} \varrho(|\phi|) d\bar{L} < \infty, $$

respectively its subspace $\mathcal{H} \subset \mathcal{B}$ of the predictable processes. Since $\mathcal{B} \subset S_{\varrho}$ and any $\phi \in \mathcal{H}$ satisfies (51), Lemma F.4 entails (54) for all $\phi \in \mathcal{H}$, as $\text{KL}(\mathbb{P}|\vec{R}) < \infty$. This implies the linear mapping $\phi \mapsto \mathbb{E}_{\mathbb{P}}(\phi \odot \tilde{\mu}_T^L)$ is continuous on $\mathcal{H}$ equipped with the norm $\|\cdot\|_{\varrho}$.

Note that the convex conjugate of the Young function $\varrho(|a|)$ is $\varrho^*(|b|)$. Thus, the dual space of $(S_{\varrho}, \|\cdot\|_{\varrho})$ is isomorphic to the space

$$ L_{\varrho^*} := \left\{ k : \Omega \times [0, T_f] \times \mathcal{Q}_d \to \mathbb{R}; \text{ measurable s.t. } \mathbb{E}_{\mathbb{P}} \int_{[0,T_f]} \sum_{q \in \mathcal{Q}_d} \varrho^*(|k|) d\bar{L} < \infty \right\}. $$

Using the Riesz theorem, there exists a function $k \in L_{\varrho^*}$ such that

$$ \mathbb{E}_{\mathbb{P}}(\phi \odot \tilde{\mu}^L) = \mathbb{E}_{\mathbb{P}} \int_{[0,T_f]} \sum_{q \in \mathcal{Q}_d} k\phi d\bar{L}, \quad \text{for } \phi \in \mathcal{H}. \quad (55) $$

We now prove the uniqueness and predictability. Introduce the predictable projection of $k \in L_{\varrho^*}$ as $k^{pr} := \mathbb{E}_{\mathbb{P}}(k|X_{[0,t)})$, for $t \in [0, T_f]$. Since $\mathcal{B}$ is dense in $S_{\varrho}$, $\mathcal{H}$ is dense in the subspace of all the predictable processes in $S_{\varrho}$. Then, any two functions $g, k \in L_{\varrho^*}$ satisfying (55) must share the same projection, i.e. $g^{pr} = k^{pr}$. It follows that there exists a unique predictable process $k$ in the space

$$ \mathcal{K}(P) := \left\{ k : \Omega \times [0, T_f] \times \mathcal{Q}_d \to \mathbb{R}; \text{ predictable s.t. } \mathbb{E}_{\mathbb{P}} \int_{[0, T_f]} \sum_{q \in \mathcal{Q}_d} \varrho^*(|k|) d\bar{L} < \infty \right\}, $$
which satisfies (55). Moreover,

$$ \phi \odot \tilde{\mu}^L - \phi \odot k\bar{L} = \phi \odot (\mu^X - \bar{L} - k\bar{L}) = \phi \odot (\mu^X - (k+1)\bar{L}) = \phi \odot (\mu^X - u\bar{L}), \quad \text{for } \phi \in \mathcal{H}, $$

with $u = k + 1$. This means that the equation (55) is equivalent to

$$ \mathbb{E}_{\mathbb{P}}[\phi \odot (\mu^X - u\bar{L})] = 0, \quad \text{for } \phi \in \mathcal{H}. \qquad (56) $$

Thus, $u\bar{L}$ is a positive measure and $u$ is nonnegative. Furthermore, $\mathbb{P} \in \text{MP}(l\bar{L})$ since equation (56) implies

$$ \mathbb{E}_{\mathbb{P}} \left[ \sum_{0 \le s \le t} \phi(s, \Delta \vec{X}_s) \right] = \mathbb{E}_{\mathbb{P}} \left[ \int_{[0, T_f]} \sum_{q \in \mathcal{Q}_d} \phi(s, q) u \bar{L}(dq ds) \right], \quad \forall \phi \in \mathcal{H}. $$

We now prove the second claim regarding the expression of the relative entropy $\text{KL}(\mathbb{P}|\vec{R})$. When $\mathbb{P} \sim \vec{R}$, we define the stopping time $\tau_j^k$ as

$$ \tau_j^k := \inf \left\{ t \in [0, T_f] ; \int_{[0, T_f]} \sum_{q \in \mathcal{Q}_d} \varrho(\log u) d\bar{L} \ge k \text{ or } (\Delta \vec{X}_t \ne 0 \text{ and } \log u_t(\Delta \vec{X}_t) \notin [-j, k]) \right\}. $$

By conditioning w.r.t. $\vec{X}_0$, we can assume without loss of generality that $\vec{R}_0 = \mathbb{P}_0$, i.e. $\frac{d\mathbb{P}_0}{d\vec{R}_0}(\vec{X}_0) = 1$. Applying Lemma F.5 with $\mathbb{P} \in \text{MP}(u\bar{L})$ and $\chi = -\log u$, we obtain

$$ \begin{aligned} Q_{\tau_j^k}^k &:= \exp \left( -\log u \odot \tilde{\mu}_1^{uL} - \int_{[0, T_f]} \sum_{q \in \mathcal{Q}_d} \varrho(-\log u) u \bar{L}(dq ds) \right) \mathbb{P}_{\tau_j^k}^k \\ &\in \text{MP}(\mathbb{1}_{[0, \tau_j^k]} e^{-\log u} u \bar{L}) = \text{MP}(\mathbb{1}_{[0, \tau_j^k]} \bar{L}). \end{aligned} \qquad (57) $$

Since $\vec{R} \tau_j^k$ fulfills the uniqueness condition (U), from (57), we derive

$$ Q_{\tau_j^k}^k = \vec{R} \tau_j^k. $$

First, we apply F.5 with $\vec{R} \in \text{MP}(\bar{L})$ and $\chi = \log u$ to get

$$ \begin{aligned} \tilde{\mathbb{P}}_{\tau_j^k}^k &:= \exp \left( \log u \odot \tilde{\mu}_1^{uL} - \int_{[0, T_f]} \sum_{q \in \mathcal{Q}_d} \varrho(\log u) \bar{L}(dq ds) \right) \vec{R}_{\tau_j^k}^k \\ &\in \text{MP}(\mathbb{1}_{[0, \tau_j^k]} e^{\log u} u \bar{L}) = \text{MP}(\mathbb{1}_{[0, \tau_j^k]} l \bar{L}). \end{aligned} $$

Secondly, applying Lemma F.5 with $\tilde{\mathbb{P}}_{\tau_j^k}^k \in \text{MP}(\mathbb{1}_{[0, \tau_j^k]} u \bar{L})$ and $\chi = -\log u$ yields

$$ \begin{aligned} \tilde{Q}_{\tau_j^k}^k &:= \exp \left( -\log u \odot \tilde{\mu}_1^{uL} - \int_{[0, T_f]} \sum_{q \in \mathcal{Q}_d} \varrho(-\log u) u \bar{L}(dq ds) \right) \tilde{\mathbb{P}}_{\tau_j^k}^k \\ &\in \text{MP}(\mathbb{1}_{[0, \tau_j^k]} e^{-\log u} u \bar{L}) = \text{MP}(\mathbb{1}_{[0, \tau_j^k]} \bar{L}). \end{aligned} $$
From the uniqueness condition (U) satisfied by $\vec{R}^{\tau_j^k}$, it follows that $\tilde{Q}^{\tau_j^k} = \vec{R}^{\tau_j^k}$. Combining it with $Q^{\tau_j^k} = \vec{R}^{\tau_j^k}$ implies

$$Q^{\tau_j^k} = \tilde{Q}^{\tau_j^k},$$

which means

$$\begin{aligned} & \exp\left(-\log u \odot \tilde{\mu}_1^{uL} - \int_{[0,T_f]} \sum_{q \in Q_d} \varrho(-\log u) u \bar{L}(dqds)\right) \mathbb{P}^{\tau_j^k} \\ & = \exp\left(-\log u \odot \tilde{\mu}_1^{uL} - \int_{[0,T_f]} \sum_{q \in Q_d} \varrho(-\log u) u \bar{L}(dqds)\right) \tilde{\mathbb{P}}^{\tau_j^k}. \end{aligned}$$

Notice that $\exp(-\log u \odot \tilde{\mu}_1^{uL} - \int_{[0,T_f]} \sum_{q \in Q_d} \varrho(-\log u) u \bar{L}(dqds)) > 0$, we finally conclude that $\mathbb{P}^{\tau_j^k} = \tilde{\mathbb{P}}^{\tau_j^k}$ i.e.

$$\mathbb{1}_{[0, \tau_j^k \wedge 1]} \frac{d\mathbb{P}}{d\vec{R}} = \mathbb{1}_{[0, \tau_j^k \wedge 1]} \frac{d\mathbb{P}_0}{d\vec{R}_0}(\vec{X}_0) \exp\left(\left(\mathbb{1}_{[0, \tau_j^k \wedge 1]} \log u\right) \odot \tilde{\mu}^L - \int_{[0, \tau_j^k \wedge 1]} \sum_{q \in Q_d} \varrho(\log u) d\bar{L}\right).$$

Letting $k$ and $j$ tend to infinity, since $\tau := \lim_{k,j \to \infty} \tau_j^k = \infty$, we get

$$\frac{d\mathbb{P}}{d\vec{R}} = \frac{d\mathbb{P}_0}{d\vec{R}_0}(\vec{X}_0) \exp\left(\int_{[0,T_f]} \sum_{q \in Q_d} \log u_t(q) \tilde{\mu}^L(dtdq) - \int_{[0,T_f]} \sum_{q \in Q_d} \varrho(\log u_t(q)) \bar{L}(dtdq)\right), \quad \mathbb{P}\text{-a.s.}$$

We now extend the result above to the case when $\mathbb{P}$ might not be equivalent to $\vec{R}$. The idea is to approximate $\mathbb{P}$ by a sequence $(\mathbb{P}_n)$, which satisfies $\mathbb{P}_n \sim \vec{R}$ for all $n \ge 1$. Denoting

$$\mathbb{P}_n = \left(1 - \frac{1}{n}\right)\mathbb{P} + \frac{\vec{R}}{n}, \quad n \ge 1,$$

we have $\mathbb{P}_n \sim \vec{R}$ and $\lim_{n \to \infty} \text{KL}(\mathbb{P}|\mathbb{P}_n) = 0$. For simplicity, write $\chi = \log u$ and $\chi^n = \log u^n$, well-defined $\mathbb{P}$-a.s. From (53) and using $\mathbb{P} \in \text{MP}(u\bar{L})$ combined with Lemma F.4, we obtain

$$\text{KL}(\mathbb{P}|\mathbb{P}_n) \ge \mathbb{E}_{\mathbb{P}}((\chi - \chi^n) \odot \tilde{\mu}^{u^n L} - \int_{[0,T_f]} \sum_{q \in Q_d} \varrho(\chi - \chi^n) u^n d\bar{L}).$$

By definition, we have

$$\tilde{\mu}^{u^n L} = \mu^X - u^n \bar{L} = \mu^X - u \bar{L} + (u - u^n) \bar{L} = \tilde{\mu}^{uL} + (u - u^n) \bar{L}, \quad \mathbb{P}\text{-a.s.},$$

which yields

$$\begin{aligned} \text{KL}(\mathbb{P}|\mathbb{P}_n) &\ge \mathbb{E}_{\mathbb{P}}((\chi - \chi^n) \odot (\tilde{\mu}^{uL} + \bar{L}(u - u^n))) - \int_{[0,T_f]} \sum_{q \in Q_d} \left(\frac{u}{u^n} - \log \frac{u}{u^n} - 1\right) u^n d\bar{L} \\ &= \mathbb{E}_{\mathbb{P}}((\chi - \chi^n) \odot \tilde{\mu}^{uL} + \int_{[0,T_f]} \sum_{q \in Q_d} \log \frac{u}{u^n} u d\bar{L} - \int_{[0,T_f]} \sum_{q \in Q_d} \left(\frac{u}{u^n} - 1\right) u^n d\bar{L}). \end{aligned}$$
Since $\mathbb{P} \in \text{MP}(u\bar{L})$, we deduce that the stochastic integral $(\chi - \chi^n) \odot \tilde{\mu}^{u\bar{L}}$ is a local $\mathbb{P}$-martingale. Therefore,

$$
\begin{align*}
\text{KL}(\mathbb{P}|\mathbb{P}_n) &\ge \mathbb{E}_{\mathbb{P}}\left(\int_{[0,T_f]} \sum_{q \in \mathcal{Q}_d} (u^n - u - u \log \frac{u^n}{u}) d\bar{L}\right) \\
&= \mathbb{E}_{\mathbb{P}}\left(\int_{[0,T_f]} \sum_{q \in \mathcal{Q}_d} \left(\frac{u^n}{u} - \log \frac{u^n}{u} - 1\right) ud\bar{L}\right) \\
&= \mathbb{E}_{\mathbb{P}}\left(\int_{[0,T_f]} \sum_{q \in \mathcal{Q}_d} \varrho(\chi^n - \chi) d(u\bar{L})\right).
\end{align*}
$$

Since $\lim_{n \to \infty} \text{KL}(\mathbb{P}|\mathbb{P}_n) = 0$, we obtain

$$
\lim_{n \to \infty} \mathbb{E}_{\mathbb{P}} \int_{[0, T_f]} \sum_{q \in \mathcal{Q}_d} \varrho(\chi^n - \chi) du \bar{L} = 0. \qquad (58)
$$

Moreover, the fact that $\mathbb{P}_n \sim \bar{\mathbb{R}}$ yields

$$
\frac{d\mathbb{P}_n}{d\bar{\mathbb{R}}} = \frac{d\mathbb{P}_{n,0}}{d\bar{\mathbb{R}}_0}(\bar{X}_0) \exp \left( \chi^n \odot \tilde{\mu}^L - \int_{[0,T_f]} \sum_{q \in \mathcal{Q}_d} \varrho(\chi^n) d\bar{L} \right).
$$

Taking the limit for $n \to \infty$, and using dominated convergence theorem, we get

$$
\text{KL}(\mathbb{P}|\vec{R}) = \text{KL}(\mathbb{P}_0|\vec{R}_0) + \lim_{n \to \infty} \mathbb{E}_{\mathbb{P}}(\chi^n \odot \tilde{\mu}^L - \int_{[0,T_f]} \sum_{q \in \mathcal{Q}_d} \varrho(\chi^n) d\bar{L}).
$$

By direct computation, we can show

$$
\begin{align*}
\mathbb{E}_{\mathbb{P}}(\chi^n \odot \tilde{\mu}^L - \int_{[0,T_f]} \sum_{q \in \mathcal{Q}_d} \varrho(\chi^n) d\bar{L}) &= \mathbb{E}_{\mathbb{P}}(\chi \odot \tilde{\mu}^L - \int_{[0,T_f]} \sum_{q \in \mathcal{Q}_d} \varrho(\chi) d\bar{L}) \\
&\quad - \mathbb{E}_{\mathbb{P}} \int_{[0,T_f]} \sum_{q \in \mathcal{Q}_d} \varrho(\chi^n - \chi) du \bar{L}.
\end{align*}
$$

Taking $n \to \infty$ and using (58), we deduce that

$$
\text{KL}(\mathbb{P}|\vec{R}) = \text{KL}(\mathbb{P}_0|\vec{R}_0) + \mathbb{E}_{\mathbb{P}}\left(\int_{[0,T_f]} \sum_{q \in \mathcal{Q}_d} \log u \, du \tilde{\mu}^L - \int_{[0,T_f]} \sum_{q \in \mathcal{Q}_d} \varrho(\log u) \, d\bar{L}\right).
$$

Applying (56) to the function $\phi = \log u$, we get

$$
\begin{align*}
\text{KL}(\mathbb{P}|\vec{R}) &= \text{KL}(\mathbb{P}_0|\vec{R}_0) + \mathbb{E}_{\mathbb{P}}\left(\int_{[0,T_f]} \sum_{q \in \mathcal{Q}_d} (l-1) \log u \, d\bar{L} - \int_{[0,T_f]} \sum_{q \in \mathcal{Q}_d} \varrho(\log u) \, d\bar{L}\right) \\
&= \text{KL}(\mathbb{P}_0|\vec{R}_0) + \mathbb{E}_{\mathbb{P}}\left(\int_{[0,T_f]} \sum_{q \in \mathcal{Q}_d} [(u-1) \log u - u + \log u + 1] \, d\bar{L}\right) \\
&= \text{KL}(\mathbb{P}_0|\vec{R}_0) + \mathbb{E}_{\mathbb{P}}\left(\int_{[0,T_f]} \sum_{q \in \mathcal{Q}_d} (u \log u - u + 1) \, d\bar{L}\right).
\end{align*}
$$
Replacing the formula of $\bar{L}$ and extending the function $u$ to $u_t(\vec{X}_{t-}, q) = 0$ if $(\vec{X}_{t-}, q) \notin \mathcal{A}$, the previous equality can be rewritten as

$$ 
\begin{aligned}
\mathrm{KL}(\mathbb{P}|\vec{R}) &= \mathrm{KL}(\mathbb{P}_0|\vec{R}_0) + \mathbb{E}_{\mathbb{P}} \int_{[0,T_f]} \sum_{q \in \mathcal{Q}_d} (u \log u - u + 1)(\vec{X}_{t-}, q)dt \\
&= \mathrm{KL}(\mathbb{P}_0|\vec{R}_0) + \mathbb{E}_{\mathbb{P}} \int_{[0,T_f]} \sum_{q \in \mathcal{Q}_d} (u \log u - u + 1)(\vec{X}_t, q)dt,
\end{aligned}
$$

where the last equality above relies on the fact that $\vec{X}_{t-} = X_t$ for Lebesgue almost all $t \in [0, T_f]$. We conclude that the relative entropy admits the following expression

$$ \mathrm{KL}(\mathbb{P}|\vec{R}) = \mathrm{KL}(\mathbb{P}_0|\vec{R}_0) + \mathbb{E}_{\mathbb{P}} \int_{[0,T_f]} \sum_{q \in \mathcal{Q}_d} h(u(\vec{X}_t, q))dt, $$

with $h(a) := \varrho^*(a-1) = a \log a - a + 1$ for $a > 0$. The proof of Theorem F.3 is then finished. $\square$

## F.1.2 Optimal Control problem of the time reversal process

In the continuous case, Conforti et al. (2024) showed that the time reversal process satisfies an optimal control problem, which characterizes its evolution and provides us a power tool to prove the convergence of the algorithm simulating the backward process. In this section, using Girsanov's theorem F.3, we aim to characterize this entropic optimization problem. Let $\vec{R}$ be a reversible path probability measure on path space, i.e. $\vec{R} = \vec{R}^*$. Let $\vec{P}^* \in \vec{P}^*$ be the forward probability measure on $[0, T_f]$ started at $\mu^*$, and denote by $\vec{P}^*$ the backward probability measure with the final distribution $\mu^*$. In Appendix B.2.3, we showed that the forward invariant measure is $\gamma^d$.

**Proposition F.6.** Let $\mathbb{P} \in \mathcal{P}(\Omega)$ verifying $\mathrm{KL}(\mathbb{P}|\vec{R}) < \infty$. We then have that its time reversal $\vec{P}^* \in \vec{P}^*$ satisfies the following optimization problem

$$ \vec{P}^* = \arg \min_{\mathbb{P} \text{ s.t. } \vec{P}_{T_f} = \mu^*} (\mathrm{KL}(\mathbb{P}|\vec{R}) + \int g d\mathbb{P}_{T_f}), \quad \text{with } g = -\log \frac{d\mu^*}{d\gamma^d}. $$

*Proof of Proposition F.6.* From Léonard (2012, Proposition 3.1), we have that the relative entropy $\mathrm{KL}(\mathbb{P}|\vec{R})$ admits the following variational representation

$$ \mathrm{KL}(\mathbb{P}|\vec{R}) = \sup_{f \text{ s.t. } \int e^f d\vec{R} < \infty} \left( \int f d\mathbb{P} - \log \int e^f d\vec{R} \right), \quad \text{for } \mathbb{P} \in \mathcal{P}(\Omega). $$

Taking $f(X.) = -g(\vec{X}_{T_f})$, for all $\mathbb{P} \in \mathcal{P}(\Omega)$ satisfying the finite relative entropy condition, we get

$$ \mathrm{KL}(\mathbb{P}|\vec{R}) \ge \int -g d\mathbb{P}_{T_f} - \log \int e^{-g} d\vec{R} = - \int g d\mathbb{P}_{T_f} - \log \int d\mu^* = - \int g d\mathbb{P}_{T_f}, $$

since $\int d\mu^* = 1$. As $\vec{P}^*$ is the backward process ended at $\mu^*$ and $\vec{R}$ is a reversible path probability measure on $[0, T_f]$, i.e. $\vec{R} = \vec{R}^*$, we have
$$ \frac{d\overleftarrow{\mathbb{P}}^{\mu^*}}{d\overrightarrow{R}} = \frac{d\overleftarrow{\mathbb{P}}^{\mu^*}}{d\overrightarrow{R}} = \frac{d\mu^*(\overleftarrow{X}_{T_f})}{d\gamma^d} = e^{-g(\overleftarrow{X}_{T_f})}. $$

This implies

$$ \mathrm{KL}(\overleftarrow{\mathbb{P}}^{\mu^*} | \overrightarrow{R}) = \mathrm{KL}(\mu^* | \gamma^d) = \int \log \frac{d\mu^*}{d\gamma^d} d\mu^* = - \int g d\mu^* = - \int g d\mathbb{P}_{T_f}. $$

Combining the previous results, we obtain that the time reversal $\overleftarrow{\mathbb{P}}^{\mu^*}$ is the optimal solution to the following problem

$$ \overleftarrow{\mathbb{P}}^{\mu^*} = \underset{\mathbb{P} \text{ s.t. } \mathbb{P}_{T_f} = \mu^*}{\arg \min} (\mathrm{KL}(\mathbb{P}|\overrightarrow{R}) + \int g d\mathbb{P}_{T_f}), $$

which is the desired conclusion.

Now, using the expression of $\mathrm{KL}(\mathbb{P}|\overrightarrow{R})$ in Girsanov's theorem F.3, we can derive the following Optimal Control problem.

**Definition F.7 (Admissible controls).** We say that a measurable function $u: \Omega \times [0, T_f] \times \mathcal{X} \times \mathcal{Q}_d \to \mathbb{R}^+$ is an admissible control if $u_t(x, q) = 0$ if $(x, q) \in \mathcal{X} \setminus \mathcal{A}$. We denote by $\mathcal{U}$ the set of all admissible controls.

**Theorem F.8.** $\overleftarrow{\mathbb{P}}^{\mu^*}$ is the law of $\overleftarrow{X}^{u^*}$ with $u^*$ is the optimal solution to

$$ V(0, \overleftarrow{X}_0) = \inf_{u \in \mathcal{U}} \mathbb{E} \left[ \int_{[0, T_f]} \sum_{q \in \mathcal{Q}_d} h(u_t(\overleftarrow{X}_t^u, q)) dt + g(\overleftarrow{X}_{T_f}^u) \right], \quad (59) \\ \text{s.t.} \quad \begin{cases} \overleftarrow{X}^u = \overleftarrow{X}_0^u + q \odot \mu_{uL}^X, \\ \overleftarrow{X}_{T_f}^u \sim \mu^*. \end{cases} $$

*Proof of Theorem F.8.* Theorem F.8 is a consequence of Theorem F.3 and Proposition F.6 together with the note that $\overleftarrow{X}_{t-}^u = \overleftarrow{X}_t^u$ for Lebesgue almost all $t \in [0, T_f]$. □

### F.1.3 Hamilton–Jacobi–Bellman equation

The goal of this section is to characterize the previous optimization problem via the Hamilton–Jacobi–Bellman (HJB) equation. To this purpose, we first consider the generalization of the previous control problem. Let $J$ be the following cost

$$ J(t, x, u) := \mathbb{E} \left[ \int_{[t, T_f]} \sum_{q \in \mathcal{Q}_d} h(u_s(\overleftarrow{X}_s^{t,x,u}, q)) ds + g(\overleftarrow{X}_{T_f}^{t,x,u}) \right], \\ \text{s.t.} \quad \begin{cases} \overleftarrow{X}_{t,x,u} = \overleftarrow{X}_t^{t,x,u} + q \odot \mu_{uL}^X, & \text{for } (x, t, u) \in \mathcal{X} \times [0, T_f] \times \mathcal{U}. \\ \overleftarrow{X}_t^{t,x,u} = x, & \end{cases} $$
Consider $V(t, x)$ to be the value function of the previous cost function, i.e.,

$$V(t, x) := \inf_{u \in \mathcal{U}} J(t, x, u).$$

The following Dynamic Programming Principle is the main tool to derive the HJB equation.

**Lemma F.9.** For any stopping time $\kappa \in [t, T_f]$, the Dynamic Programming Principle (DPP) implies

$$V(t, x) = \inf_{u \in \mathcal{U}} \mathbb{E} \left[ \int_{(t, \kappa]} \sum_{q \in \mathcal{Q}_d} h(u_s(\overleftarrow{X}_s^{t,x,u}, q))ds + V(\kappa, \overleftarrow{X}_{\kappa}^{t,x,u}) \right]. \quad (60)$$

*Proof of Lemma F.9.* Refer to Touzi (2012, Section 2.2). □

The expression of $V$ given in Lemma F.9 leads us to the following HJB equation, which is a characterization of the optimal control to the problem (59).

**Theorem F.10.** Assume that $V$ is continuously differentiable in time. Then, the optimal control $u^*$ to the problem (59) is

$$u_t^*(x, q) = \begin{cases} e^{V(t,x)-V(t,x+q)}, & \text{if } (x, q) \in \mathcal{A}, \\ 0, & \text{otherwise}, \end{cases}$$

with $V$ satisfies the following HJB equation

$$\left\{ \begin{array}{l} \partial_t V(t, x) + \inf_{u \in \mathcal{U}} \sum_{q \in \mathcal{Q}_d} [h(u_t(x, q)) + [V(t, x + q) - V(t, x)]u_t(x, q)] = 0, \\ V(T_f, x) = g(x), \end{array} \right. \quad \text{for } (t, x) \in [0, T_f] \times \mathbf{X}.$$

The HJB equation can be rewritten as

$$\left\{ \begin{array}{l} \partial_t V(t, x) - \sum_{q \in \mathcal{Q}_d} (e^{V(t,x)} - V(t,x+q)) \mathbb{1}_{\mathcal{A}}(x, q) + d = 0, \\ V(T_f, x) = g(x), \end{array} \right. \quad \text{for } (t, x) \in [0, T_f] \times \mathbf{X}. \quad (61)$$

*Proof of Theorem F.10.* The DPP formula (60) for $\kappa = t + \alpha$ with $\alpha > 0$ leads to

$$\mathbb{E} \left[ \int_{(t, t+\alpha]} \sum_{q \in \mathcal{Q}_d} h(u_s(\overleftarrow{X}_{s-}^{t,x,u}, q))ds + V(t+\alpha, \overleftarrow{X}_{t+\alpha}^{t,x,u}) - V(t, x) \right] \geq 0,$$

for any admissible control $u \in \mathcal{U}$. Using Itô's formula on the process $\overleftarrow{X}^{t,x,u}$, we get

$$\mathbb{E} \left[ \int_{(t, t+\alpha]} \sum_{q \in \mathcal{Q}_d} h(u_s(\overleftarrow{X}_{s-}^{t,x,u}, q))ds \right. \\ \left. + \int_{(t, t+\alpha]} (\partial_t V(s, \overleftarrow{X}_{s-}^{t,x,u}) + \sum_{q \in \mathcal{Q}_d} (V(s, \overleftarrow{X}_{s-}^{t,x,u} + q) - V(s, \overleftarrow{X}_{s-}^{t,x,u}))u_s(\overleftarrow{X}_{s-}^{t,x,u}, q))ds \right] \geq 0.$$
Multiplying the both hand sides by $\frac{1}{\alpha}$ and pushing $\alpha \to 0$, we derive that

$$ \sum_{q \in \mathcal{Q}_d} h(u_t(x, q)) + \partial_t V(t, x) + \sum_{q \in \mathcal{Q}_d} [V(t, x + q) - V(t, x)] u_t(x, q) \geq 0, $$

for any admissible control $u \in \mathcal{U}$. Taking the infimum w.r.t. $u$, we get

$$ \partial_t V(t, x) + \inf_{u \in \mathcal{U}} \sum_{q \in \mathcal{Q}_d} [h(u_t(x, q)) + [V(t, x + q) - V(t, x)]u_t(x, q)] \geq 0, \quad \text{for } (t, x) \in [0, T_f] \times \mathbf{X}. $$

We prove next the equality by contradiction. Assume that there exists $(t_0, x_0) \in [0, T_f] \times \mathbf{X}$ such that

$$ \partial_t V(t_0, x_0) + \inf_{u \in \mathcal{U}} \sum_{q \in \mathcal{Q}_d} [h(u_{t_0}(x_0, q)) + [V(t_0, x_0 + q) - V(t_0, x_0)]u_{t_0}(x_0, q)] > 0. $$

Denote $\Delta V(t_0, x_0, q) := V(t_0, x_0 + q) - V(t_0, x_0)$. The previous inequality implies that there exists $\varepsilon > 0$ such that

$$ \partial_t V(t_0, x_0) + \inf_{u \in \mathcal{U}} \sum_{q \in \mathcal{Q}_d} [h(u) + u\Delta V](t_0, x_0, q) \geq \varepsilon > 0. \quad (62) $$

Take $\xi > 0$ small enough such that

$$ \sum_{q \in \mathcal{Q}_d} (e^{-\Delta V + \xi} - e^{-\Delta V})(t_0, x_0, q) < \frac{\varepsilon}{2}, \quad (63) $$

and define the function $\varphi \leq V$ as

$$ \varphi(t, x) := V(t, x) - \xi [|T_f - t_0|^2 + \delta_{\{x_0\}}(x)], \quad \text{for } (t, x) \in [0, T_f] \times \mathbf{X}. $$

It is clear that

$$ \varphi(t_0, x_0) = V(t_0, x_0), \quad \partial_t \varphi(t_0, x_0) = \partial_t V(t_0, x_0), \quad \text{and} \quad \varphi(t_0, x) - V(t_0, x) = -\xi \text{ for } x \neq x_0. $$

Therefore,

$$ \begin{aligned} & \partial_t \varphi(t_0, x_0) + \inf_{u \in \mathcal{U}} \sum_{q \in \mathcal{Q}_d} [h(u_{t_0}(x_0, q)) + [\varphi(t_0, x_0 + q) - \varphi(t_0, x_0)]u_{t_0}(x_0, q)] \\ &= \partial_t V(t_0, x_0) + \inf_{u \in \mathcal{U}} \sum_{q \in \mathcal{Q}_d} [h(u_{t_0}(x_0, q)) + [V(t_0, x_0 + q) - V(t_0, x_0) - \xi] u_{t_0}(x_0, q)] \\ &= \partial_t V(t_0, x_0) + \inf_{u \in \mathcal{U}} \sum_{q \in \mathcal{Q}_d} [h(u) + (\Delta V - \xi)u](t_0, x_0, q). \end{aligned} $$
The minimum above is attained at $u \in \mathcal{U}$ such that $u(t_0, x_0, q) = e^{-\Delta V+\xi}(t_0, x_0, q)$, thus

$$
\begin{align*}
& \partial_t \varphi(t_0, x_0) + \inf_{u \in \mathcal{U}} \sum_{q \in \mathcal{Q}_d} [h(u_{t_0}(x_0, q)) + [\varphi(t_0, x_0 + q) - \varphi(t_0, x_0)]u_{t_0}(x_0, q)] \\
&= \partial_t V(t_0, x_0) + \sum_{q \in \mathcal{Q}_d} [h(e^{-\Delta V+\xi}) + (\Delta V - \xi)e^{-\Delta V+\xi}] (t_0, x_0, q) \\
&= \partial_t V(t_0, x_0) + \sum_{q \in \mathcal{Q}_d} (1 - e^{-\Delta V+\xi})(t_0, x_0, q) \\
&= \partial_t V(t_0, x_0) + \sum_{q \in \mathcal{Q}_d} (1 - e^{-\Delta V})(t_0, x_0, q) + \sum_{q \in \mathcal{Q}_d} (e^{-\Delta V} - e^{-\Delta V+\xi})(t_0, x_0, q) \\
&> \varepsilon - \frac{\varepsilon}{2} = \frac{\varepsilon}{2} > 0,
\end{align*}
$$

where the last inequality relies on (63) and (62) with $u = e^{-\Delta V}$. Therefore, we obtain

$$
\partial_t \varphi(t_0, x_0) + \inf_{u \in \mathcal{U}} \sum_{q \in \mathcal{Q}_d} [h(u_{t_0}(x_0, q)) + [\varphi(t_0, x_0 + q) - \varphi(t_0, x_0)]u_{t_0}(x_0, q)] > 0.
$$

From the continuity in time of the Hamiltonian, the previous inequality yields that

$$
\begin{equation}
\begin{split}
& \partial_t \varphi(t, x) + \inf_{u \in \mathcal{U}} \sum_{q \in \mathcal{Q}_d} [h(u_t(x, q)) + [\varphi(t, x+q) - \varphi(t, x)]u_t(x, q)] \ge 0, \\
& \hspace{2em} \text{for } (t, x) \in (t_0 - r, t_0 + r) \times \{x_0\},
\end{split}
\tag{64}
\end{equation}
$$

for some $r > 0$. Defining the stopping time $\kappa^u$ as

$$
\kappa^u := \inf \left\{ t \in (t_0, T_f] : \overleftarrow{X}_{t_0, x_0, u}^{t_0, x_0, u} \neq x_0 \right\} \wedge (t_0 + r),
$$

for an arbitrary control $u$, we have

$$
\varphi(\kappa^u, \overleftarrow{X}_{\kappa^u}^{t_0, x_0, u}) = \begin{cases} \varphi(t_0 + r, \overleftarrow{X}_{t_0+r}^{t_0, x_0, u}), & \text{if } \overleftarrow{X}_{\kappa^u-}^{t_0, x_0, u} = x_0, \\ \varphi(\kappa^u, \overleftarrow{X}_{\kappa^u}^{t_0, x_0, u}), & \text{if } \overleftarrow{X}_{\kappa^u-}^{t_0, x_0, u} \neq x_0. \end{cases}
$$

This implies that

$$
\begin{align*}
\varphi(\kappa^u, \overleftarrow{X}_{\kappa^u}^{t_0, x_0, u}) - V(\kappa^u, \overleftarrow{X}_{\kappa^u}^{t_0, x_0, u}) &= \begin{cases} -\xi r^2, & \text{if } \overleftarrow{X}_{\kappa^u-}^{t_0, x_0, u} = x_0, \\ -\xi(|\kappa^u - t_0|^2 + 1), & \text{if } \overleftarrow{X}_{\kappa^u-}^{t_0, x_0, u} \neq x_0, \end{cases} \\
&\le -\xi r^2.
\end{align*}
$$

Therefore,

$$
\begin{align*}
& \mathbb{E} \left[ \int_{(t_0, \kappa^u]} \sum_{q \in \mathcal{Q}_d} h(u_s(\overleftarrow{X}_{s-}^{t_0, x_0, u}, q)) ds + V(\kappa^u, \overleftarrow{X}_{\kappa^u}^{t_0, x_0, u}) \right] \\
\ge & \mathbb{E} \left[ \int_{(t_0, \kappa^u]} \sum_{q \in \mathcal{Q}_d} h(u_s(\overleftarrow{X}_{s-}^{t_0, x_0, u}, q)) ds + \varphi(\kappa^u, \overleftarrow{X}_{\kappa^u}^{t_0, x_0, u}) + \xi r^2 \right] \\
= & \mathbb{E} \left[ \int_{(t_0, \kappa^u]} \sum_{q \in \mathcal{Q}_d} h(u_s(\overleftarrow{X}_{s-}^{t_0, x_0, u}, q)) ds + \varphi(\kappa^u, \overleftarrow{X}_{\kappa^u}^{t_0, x_0, u}) - \varphi(t_0, x_0) \right] + \varphi(t_0, x_0) + \xi r^2.
\end{align*}
$$
Using Itô's formula and the fact that $V(t_0, x_0) = \varphi(t_0, x_0)$, we obtain

$$
\begin{align*}
& \mathbb{E} \left[ \int_{(t_0, \kappa^u]} \sum_{q \in \mathcal{Q}_d} h(u_s(\overleftarrow{X}_{s-}^{t_0, x_0, u}, q)) ds + V(\kappa^u, \overleftarrow{X}_{\kappa^u}^{t_0, x_0, u}) \right] \\
&= \mathbb{E} \int_{(t_0, \kappa^u]} \left[ \partial_t \varphi(s, \overleftarrow{X}_{s-}^{t_0, x_0, u}) + \sum_{q \in \mathcal{Q}_d} h(u_s(\overleftarrow{X}_{s-}^{t_0, x_0, u}, q)) \right. \\
& \qquad \left. + (\varphi(s, \overleftarrow{X}_{s-}^{t_0, x_0, u} + q) - \varphi(s, \overleftarrow{X}_{s-}^{t_0, x_0, u})) u_s(\overleftarrow{X}_{s-}^{t_0, x_0, u}, q) \right] ds + V(t_0, x_0) + \xi r^2.
\end{align*}
$$

This together with (64) yields

$$
\mathbb{E} \left[ \int_{(t_0, \kappa^u]} \sum_{q \in \mathcal{Q}_d} h(u_s(\overleftarrow{X}_{s-}^{t_0, x_0, u}, q)) ds + V(\kappa^u, \overleftarrow{X}_{\kappa^u}^{t_0, x_0, u}) \right] \geq V(t_0, x_0) + \xi r^2. \quad (65)
$$

Since the above control $u$ is arbitrary, (65) is indeed a contradiction to DPP formula (60).
Consequently, we can deduce the following HJB equation satisfied by the value function

$$
\left\{
\begin{array}{l}
\displaystyle \partial_t V(t,x) + \inf_{u \in \mathcal{U}} \sum_{q \in \mathcal{Q}_d} [h(u_t(x,q)) + [V(t,x+q) - V(t,x)]u_t(x,q)] = 0, \quad \text{for } (t,x) \in [0, T_f] \times \mathbf{X}. \\
\displaystyle V(T_f,x) = g(x),
\end{array}
\right.
\qquad (66)
$$

We now find the minimum in (66) as done before. By direct computation, we easily obtain the
optimal solution as follows

$$
u_t^*(x, q) = \begin{cases} e^{V(t,x)-V(t,x+q)}, & \text{for } (x,q) \in \mathcal{A}, \\ 0, & \text{otherwise.} \end{cases}
$$

Hence, replacing $u^*$ into (66) yields

$$
\left\{
\begin{array}{ll}
\partial_t V(t,x) - \sum_{q \in \mathcal{Q}_d} (e^{V(t,x)} - V(t,x+q)) \mathbb{1}_{\mathcal{A}}(x,q) + d = 0, & \text{for } (t,x) \in [0, T_f] \times \mathbf{X}. \\
V(T_f,x) = g(x), &
\end{array}
\right.
$$

and we complete the proof of Theorem F.10.

<figure><img src="image_9.png" /></figure>

The previous HJB equation will be instrumental in the proof of our convergence bound. To do
this, we first consider the following martingale and monotone property.

**Proposition F.11.** With all the notations above, $\sum_{q \in \mathcal{Q}_d} u_t^*(\overleftarrow{X}_t^{u^*}, q)$ is a martingale. Furthermore, for any $0 \le s \le t \le T_f$, it holds

$y_s \le e^{-4(t-s)} y_t,$

$$
\text{with } y_r := \mathbb{E} \left[ \sum_{q \in \mathcal{Q}_d} h(u_r^*(\overleftarrow{X}_r^{u^*}, q)) \right], \text{ for } r \in [0, T_f].
$$
Proof of Proposition F.11. Fix $t \in [0, T_f]$, we have $\overline{\mathbb{P}}^{\mu_0}(\overleftarrow{X}_{t-}^{u^*} = \overleftarrow{X}_t^{u^*}) = 1$. Applying Itô's formula on

$$ \varphi(t, \overleftarrow{X}_t^{u^*}) := \sum_{q \in \mathcal{Q}_d} u_t^*(\overleftarrow{X}_t^{u^*}, q) = \sum_{q \in \mathcal{Q}_d} e^{V(t, \overleftarrow{X}_t^{u^*}) - V(t, \overleftarrow{X}_t^{u^*} + q)}, $$

we obtain that the process

$$ \varphi(t, \overleftarrow{X}_t^{u^*}) - \varphi(0, \overleftarrow{X}_0^{u^*}) + \int_0^t \left[ \partial_t \varphi(s, \overleftarrow{X}_s^{u^*}) + \sum_{\bar{q} \in \mathcal{Q}_d} \left[ \varphi(s, \overleftarrow{X}_{s-}^{u^*} + \bar{q}) - \varphi(s, \overleftarrow{X}_{s-}^{u^*}) \right] u_s^*(\overleftarrow{X}_{s-}^{u^*}, \bar{q}) \right] ds $$

is a martingale. Since $\overleftarrow{X}_t = \overleftarrow{X}_{t-}$ for Lebesgue almost all $t \in [0, T_f]$, then the process

$$ \varphi(t, \overleftarrow{X}_t^{u^*}) - \varphi(0, \overleftarrow{X}_0^{u^*}) + \int_0^t \left[ \partial_t \varphi(s, \overleftarrow{X}_s^{u^*}) + \sum_{\bar{q} \in \mathcal{Q}_d} \left[ \varphi(s, \overleftarrow{X}_s^{u^*} + \bar{q}) - \varphi(s, \overleftarrow{X}_s^{u^*}) \right] u_s^*(\overleftarrow{X}_s^{u^*}, \bar{q}) \right] ds $$

is a martingale. Denote

$$ b_s := \partial_t \varphi(s, \overleftarrow{X}_s^{u^*}) + \sum_{\bar{q} \in \mathcal{Q}_d} \left[ \varphi(s, \overleftarrow{X}_s^{u^*} + \bar{q}) - \varphi(s, \overleftarrow{X}_s^{u^*}) \right] u_s^*(\overleftarrow{X}_s^{u^*}, \bar{q}), \quad \text{for } s \in [0, T_f]. $$

By the definition of $\varphi$ and the HJB equation (61), we get that

$$ \begin{aligned} b_s &= \sum_{q \in \mathcal{Q}_d} u_s^*(\overleftarrow{X}_s^{u^*}, q) \left[ \partial_t V(s, \overleftarrow{X}_s^{u^*}) - \partial_t V(s, \overleftarrow{X}_s^{u^*} + q) \right] \\ &\quad + \sum_{\bar{q} \in \mathcal{Q}_d} \left[ \sum_{q \in \mathcal{Q}_d} u_t^*(\overleftarrow{X}_t^{u^*} + \bar{q}, q) - \sum_{q \in \mathcal{Q}_d} u_t^*(\overleftarrow{X}_t^{u^*}, q) \right] u_s^*(\overleftarrow{X}_s^{u^*}, \bar{q}) \\ &= \sum_{q \in \mathcal{Q}_d} u_s^*(\overleftarrow{X}_s^{u^*}, q) \left[ \sum_{\bar{q} \in \mathcal{Q}_d} u_s^*(\overleftarrow{X}_s^{u^*}, \bar{q}) - \sum_{\bar{q} \in \mathcal{Q}_d} u_t^*(\overleftarrow{X}_t^{u^*} + q, \bar{q}) \right] \\ &\quad + \sum_{\bar{q} \in \mathcal{Q}_d} \left[ \sum_{q \in \mathcal{Q}_d} u_t^*(\overleftarrow{X}_t^{u^*} + \bar{q}, q) - \sum_{q \in \mathcal{Q}_d} u_t^*(\overleftarrow{X}_t^{u^*}, q) \right] u_s^*(\overleftarrow{X}_s^{u^*}, \bar{q}). \end{aligned} $$

Swapping $q$ with $\bar{q}$ in the first term implies $b_s = 0$. Therefore $\sum_{q \in \mathcal{Q}_d} u_t^*(\overleftarrow{X}_t^{u^*}, q)$ is a martingale. Combining this with the fact that $h$ is a convex function, we thus obtain $\sum_{q \in \mathcal{Q}_d} h(u_t^*(\overleftarrow{X}_t^{u^*}, q))$ is a submartingale and the monotonicity follows.

We now want to go further to see the monotonicity in details. Define

$$ f(t, \overleftarrow{X}_t^{u^*}) := \sum_{q \in \mathcal{Q}_d} h(u_t^*(\overleftarrow{X}_t^{u^*}, q)) = \sum_{q \in \mathcal{Q}_d} h(e^{V(t, \overleftarrow{X}_t^{u^*}) - V(t, \overleftarrow{X}_t^{u^*} + q)}) \mathbb{1}_{\mathcal{A}}(\overleftarrow{X}_t^{u^*}, q). $$

Applying Itô's formula on $f(t, \overleftarrow{X}_t^{u^*})$, we get that the process

$$ f(t, \overleftarrow{X}_t^{u^*}) - f(0, \overleftarrow{X}_0^{u^*}) - \int_0^t \partial_t f(s, \overleftarrow{X}_s^{u^*}) + \sum_{\bar{q} \in \mathcal{Q}_d} \left[ f(s, \overleftarrow{X}_{s-}^{u^*} + \bar{q}) - f(s, \overleftarrow{X}_{s-}^{u^*}) \right] u_s^*(\overleftarrow{X}_{s-}^{u^*}, \bar{q}) ds $$
is a martingale. Since $\overline{X}_t = \overline{X}_{t-}$ for Lebesgue almost all $t \in [0, T_f]$, then the process

$$f(t, \overline{X}_t^{u^*}) - f(0, \overline{X}_0^{u^*}) - \int_0^t \partial_t f(s, \overline{X}_s^{u^*}) + \sum_{\bar{q} \in \mathcal{Q}_d} [f(s, \overline{X}_s^{u^*} + \bar{q}) - f(s, \overline{X}_s^{u^*})] u_s^*(\overline{X}_s^{u^*}, \bar{q}) \mathrm{d}s$$

is a martingale. Denote

$$c_s := \partial_t f(s, \overline{X}_s^{u^*}) + \sum_{\bar{q} \in \mathcal{Q}_d} [f(s, \overline{X}_s^{u^*} + \bar{q}) - f(s, \overline{X}_s^{u^*})] u_s^*(\overline{X}_s^{u^*}, \bar{q}), \quad \text{for } s \in [0, T_f].$$

By definition of $f$ and $h$, we have that

$$c_s = \sum_{q \in \mathcal{Q}_d} h'(u_s^*(\overline{X}_s^{u^*}, q)) u_s^*(\overline{X}_s^{u^*}, q) [\partial_t V(s, \overline{X}_s^{u^*}) - \partial_t V(s, \overline{X}_s^{u^*} + q)] \\ + \sum_{\bar{q} \in \mathcal{Q}_d} \sum_{q \in \mathcal{Q}_d} [h(u_s^*(\overline{X}_s^{u^*} + \bar{q}, q)) - h(u_s^*(\overline{X}_s^{u^*}, q))] u_s^*(\overline{X}_s^{u^*}, \bar{q}).$$

Using the HJB equation (61) and the definition of $h$ and $u$, we obtain

$$c_s = \sum_{q \in \mathcal{Q}_d} \log u_s^*(\overline{X}_s^{u^*}, q) u_s^*(\overline{X}_s^{u^*}, q) \left[ \sum_{\bar{q} \in \mathcal{Q}_d} u_s^*(\overline{X}_s^{u^*}, \bar{q}) - \sum_{\bar{q} \in \mathcal{Q}_d} u_s^*(\overline{X}_s^{u^*} + q, \bar{q}) \right] \\ + \sum_{q \in \mathcal{Q}_d} \left[ u_s^*(\overline{X}_s^{u^*} + q, -q) \log u_s^*(\overline{X}_s^{u^*} + q, -q) - u_s^*(\overline{X}_s^{u^*} + q, -q) \right. \\ \left. - u_s^*(\overline{X}_s^{u^*}, q) \log u_s^*(\overline{X}_s^{u^*}, q) + u_s^*(\overline{X}_s^{u^*}, q) \right] u_s^*(\overline{X}_s^{u^*}, q) \\ = \sum_{q \in \mathcal{Q}_d} \left[ \log u_s^*(\overline{X}_s^{u^*}, q) (u_s^*(\overline{X}_s^{u^*}, q))^2 - \log u_s^*(\overline{X}_s^{u^*}, q) u_s^*(\overline{X}_s^{u^*}, q) u_s^*(\overline{X}_s^{u^*} + q, -q) \right. \\ \left. + \left[ u_s^*(\overline{X}_s^{u^*} + q, -q) \log u_s^*(\overline{X}_s^{u^*} + q, -q) - u_s^*(\overline{X}_s^{u^*} + q, -q) \right. \right. \\ \left. \left. - u_s^*(\overline{X}_s^{u^*}, q) \log u_s^*(\overline{X}_s^{u^*}, q) + u_s^*(\overline{X}_s^{u^*}, q) \right] u_s^*(\overline{X}_s^{u^*}, q) \right].$$

Note that for $(\overline{X}_t^{u^*}, q) \in \mathcal{A}$, we have

$$u_s^*(\overline{X}_s^{u^*} + q, -q) = e^{V(s, \overline{X}_s^{u^*} + q) - V(s, \overline{X}_s^{u^*})} = \frac{1}{e^{V(s, \overline{X}_s^{u^*}) - V(s, \overline{X}_s^{u^*})}} = \frac{1}{u_s^*(\overline{X}_s^{u^*}, q)}.$$
Therefore,

$$
\begin{align*}
c_s &= \sum_{q \in \mathcal{Q}_d} \left[ \log u_s^*(\overleftarrow{X}_s^{u^*}, q) (u_s^*(\overleftarrow{X}_s^{u^*}, q))^2 - \log u_s^*(\overleftarrow{X}_s^{u^*}, q) - \log u_s^*(\overleftarrow{X}_s^{u^*}, q) \right. \\
&\qquad \left. - 1 - (u_s^*(\overleftarrow{X}_s^{u^*}, q))^2 \log u_s^*(\overleftarrow{X}_s^{u^*}, q) + (u_s^*(\overleftarrow{X}_s^{u^*}, q))^2 \right] \\
&= \sum_{q \in \mathcal{Q}_d} \left[ (u_s^*(\overleftarrow{X}_s^{u^*}, q))^2 - 2 \log u_s^*(\overleftarrow{X}_s^{u^*}, q) - 1 \right].
\end{align*}
$$

For simplicity, we write $u_s^*(\overleftarrow{X}_s^{u^*}, q) = u$, then $c_s$ can be rewritten as

$$
c_s = \sum_{q \in \mathcal{Q}_d} (u^2 - 2 \log u - 1).
$$

We now show that $c_s \ge 4f(s, \overleftarrow{X}_s^{u^*})$, meaning that $u^2 - 2\log u - 1 \ge 4(u\log u - u + 1)$. To this
purpose, we consider a function

$$
g(u) = u^2 - 2 \log u - 1 - 4(u \log u - u + 1), \quad \text{for } u > 0.
$$

We have

$$
g'(u) = 2u - \frac{2}{u} - 4 \log u \quad \text{and} \quad g''(u) = 2 + \frac{2}{u^2} - \frac{4}{u} = 2\left(\frac{1}{u} - 1\right)^2 \ge 0.
$$

Thus, $g$ is a convex function. Combining this with the fact that $u=1$ is a solution of $g'$ implies
that the minimum of $g$ is attained at $u=1$. Therefore $g(u) \ge g(1) = 0$, for any $u > 0$. This yields
$c_s \ge 4f(s, \overline{X}_s^{u^*})$, which is our desired estimate.

Since $f(t, \overleftarrow{X}_t^{u^*}) - f(0, \overleftarrow{X}_0^{u^*}) - \int_0^t c_s ds$ is a martingale, for $0 \le s \le t \le T_f$, we have

$$
\mathbb{E} \left[ f(t, \overleftarrow{X}_t^{u^*}) - f(0, \overleftarrow{X}_0^{u^*}) - \int_0^t c_r dr \right] = \mathbb{E} \left[ f(s, \overleftarrow{X}_s^{u^*}) - f(0, \overleftarrow{X}_0^{u^*}) - \int_0^s c_r dr \right].
$$

It is equivalent to

$$
\mathbb{E} [f(t, \overleftarrow{X}_t^{u^*})] = \mathbb{E} [f(s, \overleftarrow{X}_s^{u^*})] + \mathbb{E} \left[ \int_s^t c_r dr \right].
$$

Denote $y_t := \mathbb{E}[f(t, \overleftarrow{X}_t^{u^*})]$. By the dominated convergence theorem and the fact that $c_r \ge 4f(r, \overleftarrow{X}_r^{u^*})$, we get

$$
y_t \ge y_s + \int_s^t 4y_r dr,
$$

which means $y'_t \ge 4y_t$. Combining this with Grönwall's inequality implies our claim

$$
y_t \ge e^{4(t-s)} y_s \quad \text{or} \quad y_s \le e^{-4(t-s)} y_t, \quad \text{for } 0 \le s \le t \le T_f.
$$

□
## F.2 Connection between the transition matrix and canonical process point of view

As we see in previous sections, the time reversal process can be understood not only via the backward transition matrix but also via the process corresponding to the optimal control problem. The transition matrix point of view provides an approximation of the score to simulate the backward process, which is very useful in practice. In parallel, the canonical process point of view gives us a better understanding of the evolution of the time reversal process, which allows us to show a theoretical guarantee on our algorithm. These two points of view in fact have a strong relation, which will be specified in this section.

**Proposition F.12.** We see that the optimal control $u^*$ satisfies the following relation with respect to the score function defined in (12) as

$$ u_t^*(x, q) = \begin{cases} 1 - s_t^\ell(x), & \text{if } (x, q) \in \mathcal{A}, \\ 0, & \text{otherwise}, \end{cases} \quad \text{with } \ell = 1, \dots, d \text{ such that } q^\ell \neq 0. \quad (67) $$

*Proof of Proposition F.12.* The proof relies on the fact that the function $\psi(t, x) = -\log \frac{d\mu_{T_f-t}}{d\gamma^d}(x)$, for $(t, x) \in [0, T_f] \times \mathcal{X}$ is a solution to the HJB equation (61). Indeed, for $(t, x) \in [0, T_f] \times \mathcal{X}$, we have that

$$ \begin{aligned} \partial_t \psi(t, x) - \sum_{q \in \mathcal{Q}_d} e^{\psi(t,x) - \psi(t,x+q)} \mathbb{1}_{\mathcal{A}}(x, q) &= \frac{\partial_t \mu_{T_f-t}(x)}{\mu_{T_f-t}(x)} - \frac{1}{\mu_{T_f-t}(x)} \sum_{q \in \mathcal{Q}_d} \mu_{T_f-t}(x+q) \mathbb{1}_{\mathcal{A}}(x, q) \\ &= \frac{1}{\mu_{T_f-t}(x)} (\partial_t \mu_{T_f-t}(x) - \sum_{q \in \mathcal{Q}_d} \mu_{T_f-t}(x+q) \mathbb{1}_{\mathcal{A}}(x, q)). \end{aligned} $$

Using the formula of marginal distribution in (10), we can compute the numerator as follows

$$ \begin{aligned} & \partial_t \mu_{T_f-t}(x) - \sum_{q \in \mathcal{Q}_d} \mu_{T_f-t}(x+q) \mathbb{1}_{\mathcal{A}}(x, q) \\ &= \sum_{z \in \mathcal{X}} \mu_0(z) \sum_{i=1}^{d} \left( \partial_t \overrightarrow{p}_{T_f-t}^{1}(z^i, x^i) \prod_{\substack{j=1 \\ j \neq i}}^{d} \overrightarrow{p}_{T_f-t}^{1}(z^j, x^j) \right) - \sum_{z \in \mathcal{X}} \sum_{i=1}^{d} \mu_0(z) \overrightarrow{p}_{T_f-t}(z, \varphi_i(x)), \end{aligned} $$

with the map $\varphi_i: \mathcal{X} \to \mathcal{X}$ defined for any $x \in \mathcal{X}$ as the vector where the $i$-bit is flipped compared to $x$. Therefore

$$ \begin{aligned} & \partial_t \mu_{T_f-t}(x) - \sum_{q \in \mathcal{Q}_d} \mu_{T_f-t}(x+q) \mathbb{1}_{\mathcal{A}}(x, q) \\ &= \sum_{z \in \mathcal{X}} \mu_0(z) \sum_{i=1}^{d} \left( \overrightarrow{\partial_t p}_{T_f-t}^{1}(z^i, x^i) \overrightarrow{p}_{T_f-t}(z, x) - \overrightarrow{\partial_t p}_{T_f-t}^{1}(z^i, \varphi_i^i(x)) \overrightarrow{p}_{T_f-t}^{1}(z, x) \right) \\ &= \sum_{z \in \mathcal{X}} \mu_0(z) \overrightarrow{p}_{T_f-t}^{1}(z, x) \sum_{i=1}^{d} \frac{\partial_t \overrightarrow{p}_{T_f-t}^{1}(z^i, x^i) - \overrightarrow{p}_{T_f-t}^{1}(z^i, \varphi_i^i(x))}{\overrightarrow{p}_{T_f-t}^{1}(z^i, x^i)}. \end{aligned} $$
We know further that either $z^i = x^i$ or $z^i = \varphi_i^i(x)$, and in both cases, together with the formula of $\vec{p}_{T_f-t}^1$, we obtain that

$$ \frac{\partial_t \vec{p}_{T_f-t}^1(z^i, x^i) - \vec{p}_{T_f-t}^1(z^i, \varphi_i^i(x))}{\vec{p}_{T_f-t}^1(z^i, x^i)} = -1. $$

Replacing this into the previous equality, we get

$$ \partial_t \mu_{T_f-t}(x) - \sum_{q \in \mathcal{Q}_d} \mu_{T_f-t}(x+q) \mathbb{1}_{\mathcal{A}}(x, q) = -d\mu_{T_f-t}(x), $$

and therefore

$$ \partial_t \psi(t, x) - \sum_{q \in \mathcal{Q}_d} e^{\psi(t,x) - \psi(t,x+q)} \mathbb{1}_{\mathcal{A}}(x, q) = -d. $$

Moreover, $\psi$ also satisfies the final condition

$$ \psi(T_f, x) = -\log \frac{d\mu^*}{d\gamma^d}(x) = g(x), $$

meaning that $\varphi$ solves the HJB equation (61). Consequently, the optimal control is

$$ u_t^*(x, q) = \begin{cases} \frac{\mu_{T_f-t}(x+q)}{\mu_{T_f-t}(x)}, & \text{if } (x, q) \in \mathcal{A}, \\ 0, & \text{otherwise}. \end{cases} $$

In other words, the optimal control admits the following formula

$$ u_t^*(x, q) = \begin{cases} \frac{\mu_{T_f-t}(\varphi^{(\ell)}(x))}{\mu_{T_f-t}(x)}, & \text{if } (x, q) \in \mathcal{A}, \\ 0, & \text{otherwise}, \end{cases} \quad \text{with } \ell = 1, \dots, d \text{ such that } q^\ell \neq 0. $$

This implies the relation between the optimal control and the score function as

$$ u_t^*(x, q) = \begin{cases} 1 - s_t^\ell(x), & \text{if } (x, q) \in \mathcal{A}, \\ 0, & \text{otherwise}, \end{cases} \quad \text{with } \ell = 1, \dots, d \text{ such that } q^\ell \neq 0, $$

which means the transition matrix and the canonical process point of view are equivalent. $\square$

## F.3 Convergence of DMPMs

Based on the canonical process point of view, we can characterize the backward evolution using martingale and optimal control problems. These tools equip us to prove the error's bound in Theorem 2.3.
Proof of Theorem 2.3. We show first the bound for the "distance" between the backward path measure in the continuous time $\overleftarrow{\mathbb{P}}^{\mu^*}$ of $\overleftarrow{\mathbb{X}}^{u^*}$ (we write $\overleftarrow{\mathbb{P}}$ for short) and the path measure $\overleftarrow{\mathbb{P}}^{\star}$ of the simulated backward process $\overleftarrow{\mathbb{X}}^{\star}$ generated in Algorithm 1. By Girsanov's theorem F.3, we have

$$ \frac{d\overleftarrow{\mathbb{P}}}{d\overrightarrow{\mathbb{R}}} = \frac{d\overleftarrow{\mathbb{P}}_0}{d\overrightarrow{\mathbb{R}}_0}(\overleftarrow{\mathbb{X}}_0) \exp \left( \log u_t^* \odot \tilde{\mu}^L - \int_{[0, T_f]} \sum_{q \in \mathcal{Q}_d} \varrho(\log u_t^*) d\bar{L} \right). $$

With a partition $0 = t_0 < \dots < t_K = T_f$ for $M \ge 1$ of $[0, T_f]$ and $\tau = \max\{h_k, k = 1, \dots, K\}$, the previous expression implies

$$ \frac{d\overleftarrow{\mathbb{P}}}{d\overrightarrow{\mathbb{R}}} = \frac{d\overleftarrow{\mathbb{P}}_0}{d\overrightarrow{\mathbb{R}}_0}(\overleftarrow{\mathbb{X}}_0) \exp \sum_{k=0}^{K-1} \left( \int_{[t_k, t_{k+1}]} \sum_{q \in \mathcal{Q}_d} \log u_t^* d\tilde{\mu}^L - \int_{[t_k, t_{k+1}]} \sum_{q \in \mathcal{Q}_d} \varrho(\log u_t^*) d\bar{L} \right). $$

Apply Girsanov's theorem F.3 again for the path measure $\overleftarrow{\mathbb{P}}^{\star}$ of the process $\overleftarrow{\mathbb{X}}^{\star}$ in (1), we obtain

$$ \frac{d\overleftarrow{\mathbb{P}}^{\star}}{d\overrightarrow{\mathbb{R}}} = \frac{d\overleftarrow{\mathbb{P}}_0^{\star}}{d\overrightarrow{\mathbb{R}}_0}(\overleftarrow{\mathbb{X}}_0) \exp \sum_{k=0}^{K-1} \left( \int_{[t_k, t_{k+1}]} \sum_{q \in \mathcal{Q}_d} \log u_{t_k}^{\theta^*} d\tilde{\mu}^L - \int_{[t_k, t_{k+1}]} \sum_{q \in \mathcal{Q}_d} \varrho(\log u_{t_k}^{\theta^*}) d\bar{L} \right). $$

Combining the previous quantities, we see that

$$ \frac{d\overleftarrow{\mathbb{P}}}{d\overrightarrow{\mathbb{P}^{\star}}} = \frac{d\overleftarrow{\mathbb{P}}_0}{d\overrightarrow{\mathbb{P}}_0^{\star}}(\overleftarrow{\mathbb{X}}_0) \exp \sum_{k=0}^{K-1} \left( \int_{[t_k, t_{k+1}]} \sum_{q \in \mathcal{Q}_d} (\log u_t - \log u_{t_k}^{\theta^*}) d\tilde{\mu}^L - \int_{[t_k, t_{k+1}]} \sum_{q \in \mathcal{Q}_d} (\varrho(\log u_t) - \varrho(\log u_{t_k}^{\theta^*})) d\bar{L} \right). $$

This leads to the following expression of the relative entropy

$$ \mathrm{KL}(\overleftarrow{\mathbb{P}} | \overrightarrow{\mathbb{P}}^{\star}) = \mathrm{KL}(\mu_{T_f} | \gamma^d) + \sum_{k=0}^{K-1} \mathrm{E}_{\overleftarrow{\mathbb{P}}} \left[ \int_{[t_k, t_{k+1}]} \sum_{q \in \mathcal{Q}_d} (\log u_t - \log u_{t_k}^{\theta^*}) d\tilde{\mu}^L - \int_{[t_k, t_{k+1}]} \sum_{q \in \mathcal{Q}_d} (\varrho(\log u_t) - \varrho(\log u_{t_k}^{\theta^*})) d\bar{L} \right]. $$
Using equation (55) and the definition of $\rho$, we derive

$$
\begin{align}
\mathrm{KL}(\overleftarrow{\mathbb{P}} | \overrightarrow{\mathbb{P}}^*) &= \mathrm{KL}(\mu_{T_f} | \gamma^d) + \sum_{k=0}^{K-1} \mathbb{E}_{\overleftarrow{\mathbb{P}}} \left[ \int_{[t_k, t_{k+1}]} \sum_{q \in \mathcal{Q}_d} \left( \begin{aligned}[t] & (\log u_t - \log u_{t_k}^{\theta^*})(u_t - 1) \\ & - (u_t - \log u_t - 1) + (u_{t_k}^{\theta^*} - \log u_{t_k}^{\theta^*} - 1) \end{aligned} \right) \mathrm{d}\bar{L} \right] \nonumber \\
&= \mathrm{KL}(\mu_{T_f} | \gamma^d) + \sum_{k=0}^{K-1} \mathbb{E}_{\overleftarrow{\mathbb{P}}} \left[ \int_{[t_k, t_{k+1}]} \sum_{q \in \mathcal{Q}_d} \left( \frac{u_t}{u_{t_k}^{\theta^*}} \log \frac{u_t}{u_{t_k}^{\theta^*}} - \frac{u_t}{u_{t_k}^{\theta^*}} + 1 \right) u_{t_k}^{\theta^*} \mathrm{d}\bar{L} \right] \nonumber \\
&= \mathrm{KL}(\mu_{T_f} | \gamma^d) + \sum_{k=0}^{K-1} \mathbb{E}_{\overleftarrow{\mathbb{P}}} \left[ \int_{[t_k, t_{k+1}]} \sum_{q \in \mathcal{Q}_d} u_{t_k}^{\theta^*} h \left( \frac{u_t}{u_{t_k}^{\theta^*}} \right) \mathrm{d}t \right] \nonumber \\
&= \mathrm{KL}(\mu_{T_f} | \gamma^d) + \sum_{k=0}^{K-1} \mathbb{E}_{\overleftarrow{\mathbb{P}}} \left[ \int_{[t_k, t_{k+1}]} \sum_{q \in \mathcal{Q}_d} u_{t_k} h \left( \frac{u_t}{u_{t_k}} \right) \mathrm{d}t \right] \nonumber \\
&\qquad + \sum_{k=0}^{K-1} \mathbb{E}_{\overleftarrow{\mathbb{P}}} \left[ \int_{[t_k, t_{k+1}]} \sum_{q \in \mathcal{Q}_d} \left( u_{t_k}^{\theta^*} h \left( \frac{u_t}{u_{t_k}^{\theta^*}} \right) - u_{t_k} h \left( \frac{u_t}{u_{t_k}} \right) \right) \mathrm{d}t \right]. \tag{68}
\end{align}
$$

By definition of the function $h$ and the tower property, the last term can be computed as

$$
\begin{align*}
I := & \sum_{k=0}^{K-1} \mathbb{E}_{\overleftarrow{\mathbb{P}}} \left[ \int_{[t_k, t_{k+1}]} \sum_{q \in \mathcal{Q}_d} \left( u_{t_k}^{\theta^*} h \left( \frac{u_t}{u_{t_k}^{\theta^*}} \right) - u_{t_k} h \left( \frac{u_t}{u_{t_k}} \right) \right) \mathrm{d}t \right] \\
= & \sum_{k=0}^{K-1} \mathbb{E}_{\overleftarrow{\mathbb{P}}} \left[ \int_{[t_k, t_{k+1}]} \sum_{q \in \mathcal{Q}_d} \left( u_t \log \frac{u_{t_k}}{u_{t_k}^{\theta^*}} + u_{t_k}^{\theta^*} - u_{t_k} \right) \mathrm{d}t \right] \\
= & \sum_{k=0}^{K-1} \mathbb{E}_{\overleftarrow{\mathbb{P}}} \left[ \int_{[t_k, t_{k+1}]} \mathbb{E}_{\overleftarrow{\mathbb{P}}} \left[ \sum_{q \in \mathcal{Q}_d} u_t \log \frac{u_{t_k}}{u_{t_k}^{\theta^*}} \middle| \mathcal{F}_{t_k} \right] + \sum_{q \in \mathcal{Q}_d} (u_{t_k}^{\theta^*} - u_{t_k}) \mathrm{d}t \right],
\end{align*}
$$

with $\mathcal{F}_{t_k}$ the $\sigma$-algebra of $\overleftarrow{X}_{t_k}$. Proposition F.11 implies that $\sum_{q \in \mathcal{Q}_d} u_t$ is a $\overleftarrow{\mathbb{P}}$-martingale, hence

$$
\begin{align*}
I &= \sum_{k=0}^{K-1} \mathbb{E}_{\overleftarrow{\mathbb{P}}} \left[ \int_{[t_k, t_{k+1}]} \sum_{q \in \mathcal{Q}_d} \left( u_{t_k} \log \frac{u_{t_k}}{u_{t_k}^{\theta^*}} + u_{t_k}^{\theta^*} - u_{t_k} \right) \mathrm{d}t \right] \\
&= \sum_{k=0}^{K-1} \mathbb{E}_{\overleftarrow{\mathbb{P}}} \left[ \int_{[t_k, t_{k+1}]} \sum_{q \in \mathcal{Q}_d} u_{t_k}^{\theta^*} h \left( \frac{u_{t_k}}{u_{t_k}^{\theta^*}} \right) \mathrm{d}t \right] \\
&= \sum_{k=0}^{K-1} (t_{k+1} - t_k) \mathbb{E}_{\overleftarrow{\mathbb{P}}} \left[ \sum_{q \in \mathcal{Q}_d} u_{t_k}^{\theta^*} h \left( \frac{u_{t_k}}{u_{t_k}^{\theta^*}} \right) \right].
\end{align*}
$$
Assumption 2.1 combined with the fact that the backward partition $(t_k)_{k=0}^K$ is associated with the forward sampling imply

$$I \le \epsilon \sum_{k=0}^{K-1} (t_{k+1} - t_k) = \epsilon(t_N - t_0) = \epsilon T_f,$$

since $\sum_{k=0}^{K-1} (t_{k+1} - t_k)$ is a telescoping sum. Replacing this into (68) yields

$$\mathrm{KL}(\widehat{\mathbb{P}}|\widehat{\mathbb{P}}^*) \le \mathrm{KL}(\mu_{T_f}|\gamma^d) + \sum_{k=0}^{K-1} \mathrm{E}_{\widehat{\mathbb{P}}} \left[ \int_{[t_k, t_{k+1}]} \sum_{q \in \mathcal{Q}_d} u_{t_k} h \left( \frac{u_t}{u_{t_k}} \right) dt \right] + \epsilon T_f.$$

We now use the tower property of conditional expectations combined with the monotonicity showed in Proposition F.11 to bound the second term above as follows

$$\begin{align*}
\mathrm{KL}(\widehat{\mathbb{P}}|\widehat{\mathbb{P}}^*) &\le \mathrm{KL}(\mu_{T_f}|\gamma^d) + \sum_{k=0}^{K-1} \mathrm{E}_{\widehat{\mathbb{P}}} \left[ \int_{[t_k, t_{k+1}]} \mathrm{E}_{\widehat{\mathbb{P}}} \left[ \sum_{q \in \mathcal{Q}_d} u_{t_k} h \left( \frac{u_t}{u_{t_k}} \right) \middle| \mathcal{F}_{t_k} \right] dt \right] + \epsilon T_f \\
&\le \mathrm{KL}(\mu_{T_f}|\gamma^d) + \sum_{k=0}^{K-1} \mathrm{E}_{\widehat{\mathbb{P}}} \left[ \int_{[t_k, t_{k+1}]} \mathrm{E}_{\widehat{\mathbb{P}}} \left[ \sum_{q \in \mathcal{Q}_d} u_{t_k} h \left( \frac{u_{t_{k+1}}}{u_{t_k}} \right) \middle| \mathcal{F}_{t_k} \right] dt \right] + \epsilon T_f \\
&\le \mathrm{KL}(\mu_{T_f}|\gamma^d) + \sum_{k=0}^{K-1} (t_{k+1} - t_k) \mathrm{E}_{\widehat{\mathbb{P}}} \left[ \sum_{q \in \mathcal{Q}_d} u_{t_k} h \left( \frac{u_{t_{k+1}}}{u_{t_k}} \right) \right] + \epsilon T_f \\
&\le \mathrm{KL}(\mu_{T_f}|\gamma^d) + \tau \sum_{k=0}^{K-1} \mathrm{E}_{\widehat{\mathbb{P}}} \left[ \sum_{q \in \mathcal{Q}_d} u_{t_k} h \left( \frac{u_{t_{k+1}}}{u_{t_k}} \right) \right] + \epsilon T_f,
\end{align*}$$

where the last inequality comes from the fact that $\tau = \max\{h_k, k = 1, \dots, K\}$. By definition of the function $h$, we have

$$h \left( \frac{u_{t_{k+1}}}{u_{t_k}} \right) = \frac{1}{u_{t_k}} \left( h(u_{t_{k+1}}) - h(u_{t_k}) - (u_{t_k} - u_{t_{k+1}}) \log(u_{t_k}) \right),$$

which leads to

$$\mathrm{KL}(\widehat{\mathbb{P}}|\widehat{\mathbb{P}}^*) \le \mathrm{KL}(\mu_{T_f}|\gamma^d) + \tau \sum_{k=0}^{K-1} \left( \mathrm{E}_{\widehat{\mathbb{P}}} \left[ \sum_{q \in \mathcal{Q}_d} h(u_{t_{k+1}}) \right] - \mathrm{E}_{\widehat{\mathbb{P}}} \left[ \sum_{q \in \mathcal{Q}_d} h(u_{t_k}) \right] \right) \\ - \tau \sum_{k=0}^{K-1} \left( \mathrm{E}_{\widehat{\mathbb{P}}} \left[ \sum_{q \in \mathcal{Q}_d} (u_{t_k} - u_{t_{k+1}}) \log(u_{t_k}) \right] \right) + \epsilon T_f.$$

Using the tower property and the fact that $\sum_{q \in \mathcal{Q}_d} u_t$ is a $\widehat{\mathbb{P}}$-martingale proved in Proposition F.11, we obtain that

$$\mathrm{E}_{\widehat{\mathbb{P}}} \left[ \sum_{q \in \mathcal{Q}_d} (u_{t_k} - u_{t_{k+1}}) \log(u_{t_k}) \right] = \mathrm{E}_{\widehat{\mathbb{P}}} \left[ \left( \sum_{q \in \mathcal{Q}_d} u_{t_k} \log(u_{t_k}) - \mathrm{E}_{\widehat{\mathbb{P}}} \left[ \sum_{q \in \mathcal{Q}_d} u_{t_{k+1}} \log(u_{t_k}) \middle| \mathcal{F}_{t_k} \right] \right) \right] = 0.$$
Combine this with the fact that $\gamma^d$ is the invariant measure of the forward process, that satisfies a log-Sobolev inequality and therefore an exponential entropy decays (Bakry et al., 2014, Theorem 5.2.1), we get

$$KL(\overleftarrow{\mathbb{P}} | \overleftarrow{\mathbb{P}}^*) \le e^{-T_f} KL(\mu^* | \gamma^d) + \tau \sum_{k=0}^{K-1} \left( \mathbb{E}_{\overleftarrow{\mathbb{P}}} \left[ \sum_{q \in \mathcal{Q}_d} h(u_{t_{k+1}}) \right] - \mathbb{E}_{\overleftarrow{\mathbb{P}}} \left[ \sum_{q \in \mathcal{Q}_d} h(u_{t_k}) \right] \right) + \epsilon T_f.$$

We obtain a telescoping sum on the right hand side, which yields

$$KL(\overleftarrow{\mathbb{P}} | \overleftarrow{\mathbb{P}}^*) \le e^{-T_f} KL(\mu^* | \gamma^d) + \tau \left( \mathbb{E}_{\overleftarrow{\mathbb{P}}} \left[ \sum_{q \in \mathcal{Q}_d} h(u_{t_M}) \right] - \mathbb{E}_{\overleftarrow{\mathbb{P}}} \left[ \sum_{q \in \mathcal{Q}_d} h(u_{t_0}) \right] \right) + \epsilon T_f \\ = e^{-T_f} KL(\mu^* | \gamma^d) + \tau \left( \mathbb{E}_{\overleftarrow{\mathbb{P}}} \left[ \sum_{q \in \mathcal{Q}_d} h(u_{T_f}) \right] - \mathbb{E}_{\overleftarrow{\mathbb{P}}} \left[ \sum_{q \in \mathcal{Q}_d} h(u_{t_0}) \right] \right) + \epsilon T_f.$$

Since $h \ge 0$, we obtain

$$KL(\overleftarrow{\mathbb{P}} | \overleftarrow{\mathbb{P}}^*) \le e^{-T_f} KL(\mu^* | \gamma^d) + \tau \mathbb{E}_{\overleftarrow{\mathbb{P}}} \left[ \sum_{q \in \mathcal{Q}_d} h(u_{T_f}) \right] + \epsilon T_f \\ = e^{-T_f} KL(\mu^* | \gamma^d) + \tau \beta_{\gamma^d}(\mu^*) + \epsilon T_f. \qquad (69)$$

Finally, notice that $\mu^* = \text{Law}(\overleftarrow{X}_T^{u^*})$, therefore

$$KL(\mu^*|\text{Law}(\overleftarrow{X}_{T_f}^{u^*})) = KL(\text{Law}(\overleftarrow{X}_{T_f}^{u^*})|\text{Law}(X_{T_f}^{u^*})) \le KL(\text{Law}(\overleftarrow{X}_{T_f}^{u^*})|\text{Law}(\overleftarrow{X}_{T_f}^{u^*})) = KL(\overleftarrow{\mathbb{P}}|\overleftarrow{\mathbb{P}}^*),$$

where the inequality is known as *Data processing inequality* for relative entropy (Nutz, 2021, Lemma 1.6). Combining this with (69), we can conclude that

$$KL(\mu^*|\text{Law}(\overleftarrow{X}_{T_f}^{u^*})) \le e^{-T_f} KL(\mu^*|\gamma^d) + \tau \beta_{\gamma^d}(\mu^*) + \epsilon T_f,$$

and finish the proof of the convergence bound. □

## F.4 Convergence of DMPMs with early stopping strategy

*Proof of Proposition 2.6.* Recall the definition of total variation distance of $\mu_\eta$ and $\mu^*$ for any $\eta > 0$ is

$$\|\mu_\eta - \mu^*\|_{\text{TV}} = \sum_{x \in \mathcal{X}} |\mu_\eta(x) - \mu^*(x)|.$$

By the triangle inequality, we obtain

$$\|\mu_\eta - \mu^*\|_{\text{TV}} \le \sum_{x \in \mathcal{X}} |\mu_\eta(x) - \mu^*(x) \overrightarrow{p}_\eta(x, x)| + |\mu^*(x) - \mu^*(x) \overrightarrow{p}_\eta(x, x)|,$$

where the transition probability $\overrightarrow{p}_\eta$ is defined in (9). The two terms above are nonnegative as

$$\mu_\eta(x) = \sum_{z \in \mathcal{X}} \mu^*(z) \overrightarrow{p}_\eta(z, x) \ge \mu^*(x) \overrightarrow{p}_\eta(x, x),$$
and the transition probability $\vec{p}_{\eta}(x, x) \le 1$ for any $x \in \mathcal{X}$. This together with the formula of $\vec{p}_{\eta}$ in (9) yield

$$ \begin{aligned} \|\mu_{\eta} - \mu^{*}\|_{\text{TV}} &\le \sum_{x \in \mathcal{X}} \left[ \mu_{\eta}(x) + \mu^{*}(x) - 2\mu^{*}(x) \left( \frac{1}{2} + \frac{1}{2}e^{-2\lambda\eta} \right)^{d} \right] \\ &\le 2 - 2 \left( \frac{1}{2} + \frac{1}{2}e^{-2\lambda\eta} \right)^{d}. \end{aligned} $$

To simplify this upper bound, we use the exponential inequality for $e^{-2\lambda\eta}$ as follows

$$ \|\mu_{\eta} - \mu^{*}\|_{\text{TV}} \le 2 - 2 \left( \frac{1}{2} + \frac{1}{2}(-2\lambda\eta + 1) \right)^{d} = 2 - 2(1 - \lambda\eta)^{d}, $$

and the proof is completed.

<figure><img src="image_10.png" /></figure>

*Proof of Corollary 2.7.* We evaluate first the behavior of the Fisher-like information at time $\eta$ as follows

$$ \begin{aligned} \beta_{\gamma^d}(\mu_{\eta}) &= \mathbb{E} \left[ \sum_{\ell=1}^{d} h \left( e^{-\log\left(\frac{d\mu_{\eta}(\vec{X}_{\eta})}{d\gamma^d}\right)} + \log\left(\frac{d\mu_{\eta}(\varphi^{(\ell)}(\vec{X}_{\eta}))}{d\gamma^d}\right) \right) \right] \\ &= \mathbb{E} \left[ \sum_{\ell=1}^{d} h \left( e^{\log\left(\frac{d\mu_{\eta}(\varphi^{(\ell)}(\vec{X}_{\eta}))}{d\mu_{\eta}(\vec{X}_{\eta})}\right)} \right) \right] \\ &= \mathbb{E} \left[ \sum_{\ell=1}^{d} h \left( \frac{d\mu_{\eta}(\varphi^{(\ell)}(\vec{X}_{\eta}))}{d\mu_{\eta}(\vec{X}_{\eta})} \right) \right]. \end{aligned} $$

By definition of the function $h$, we have

$$ \begin{aligned} \beta_{\gamma^d}(\mu_{\eta}) &= \mathbb{E} \left[ \sum_{\ell=1}^{d} \left( \frac{d\mu_{\eta}(\varphi^{(\ell)}(\vec{X}_{\eta}))}{d\mu_{\eta}(\vec{X}_{\eta})} \log \frac{d\mu_{\eta}(\varphi^{(\ell)}(\vec{X}_{\eta}))}{d\mu_{\eta}(\vec{X}_{\eta})} - \frac{d\mu_{\eta}(\varphi^{(\ell)}(\vec{X}_{\eta}))}{d\mu_{\eta}(\vec{X}_{\eta})} + 1 \right) \right] \\ &\le \mathbb{E} \left[ \sum_{\ell=1}^{d} \left( \frac{1}{d\mu_{\eta}(\vec{X}_{\eta})} \log \frac{1}{d\mu_{\eta}(\vec{X}_{\eta})} + 1 \right) \right] \\ &= d + d \sum_{x \in \mathcal{X}} \log \frac{1}{d\mu_{\eta}(x)}. \end{aligned} $$

Using the inequality $\log \frac{1}{d\mu_{\eta}(x)} \le \frac{1}{d\mu_{\eta}(x)} - 1$, we obtain

$$ \beta_{\gamma^d}(\mu_{\eta}) \le d + d \sum_{x \in \mathcal{X}} \left( \frac{1}{d\mu_{\eta}(x)} - 1 \right) \le d \sum_{x \in \mathcal{X}} \frac{1}{d\mu_{\eta}(x)}. $$

Replacing the marginal distribution $\mu_{\eta}$ above by (10) and note that $\vec{p}_{t}^{1}(x, y) \ge \frac{1}{2} - \frac{1}{2}e^{-2\lambda t}$ for any $t \in (0, T_f)$ and $x, y \in \{0, 1\}$, we get

$$ \beta_{\gamma^d}(\mu_{\eta}) \le d \sum_{x \in \mathcal{X}} \frac{1}{\sum_{z \in \mathcal{X}} \mu^*(z) \left(\frac{1}{2} - \frac{1}{2}e^{-2\lambda\eta}\right)^d} = d \sum_{x \in \mathcal{X}} \frac{1}{\left(\frac{1}{2} - \frac{1}{2}e^{-2\lambda\eta}\right)^d} = \frac{2^{2d}d}{(1 - e^{-2\lambda\eta})^d}. $$
We now use the exponential inequality to arrive at

$$ \beta_{\gamma^d}(\mu_\eta) \le \frac{2^{2d} d}{\left(\frac{2\lambda\eta}{2\lambda\eta+1}\right)^d} = \frac{2^d d(2\lambda\eta + 1)^d}{(\lambda\eta)^d}. $$

Multiplying both hand sides with the fixed step size $\bar{h}$ conditioned by (28) implies

$$ \bar{h}\beta_{\gamma^d}(\mu_\eta) \le \frac{\varepsilon^2}{2}. \qquad (70) $$

Return to the main estimate of Corollary 2.7 and apply the triangle inequality, we have

$$ \begin{aligned} \|\mu^* - \text{Law}(\overleftarrow{X}_{T_f-\eta}^*)\|_{\text{TV}} &\le \|\mu^* - \mu_\eta\|_{\text{TV}} + \|\mu_\eta - \text{Law}(\overleftarrow{X}_{T_f-\eta}^*)\|_{\text{TV}} \\ &\le \|\mu^* - \mu_\eta\|_{\text{TV}} + \sqrt{2\text{KL}(\mu_\eta|\text{Law}(\overleftarrow{X}_{T_f-\eta}^*))}, \end{aligned} $$

where the last line used Pinsker's inequality. Then using Theorem 2.5 and Proposition 2.6 implies

$$ \|\mu^* - \text{Law}(\overleftarrow{X}_{T_f-\eta}^*)\|_{\text{TV}} \le 2 - 2(1 - \lambda\eta)^d + \sqrt{2e^{-T_f} \text{KL}(\mu^*|\gamma^d)} + \sqrt{2\bar{h}\beta_{\gamma^d}(\mu_\eta)} + \sqrt{2\epsilon(T_f - \eta)}. \qquad (71) $$

Furthermore, the conditions of $\eta$ and $\bar{K}_f$ in (28) follow

$$ 2 - 2(1 - \lambda\eta)^d \le \varepsilon \quad \text{and} \quad 2e^{-T_f} \text{KL}(\mu^*|\gamma^d) \le \frac{\varepsilon}{2}. \qquad (72) $$

Finally, substituting estimates (70) and (72) into (71) yields the desired claim

$$ \|\mu^* - \text{Law}(\overleftarrow{X}_{T_f-\eta}^*)\|_{\text{TV}} \le 2\varepsilon + \sqrt{2\epsilon T_f}, $$

and the proof is finished. $\square$