Discrete Diffusion Schrödinger Bridge Matching for Graph Transformation

## Abstract

Abstract Transporting between arbitrary distributions is a fundamental goal in generative modeling.
Recently proposed diffusion bridge models provide a potential solution, but they rely on a joint distribution that is difficult to obtain in practice.
Furthermore, formulations based on continuous domains limit their applicability to discrete domains such as graphs.
To overcome these limitations, we propose Discrete Diffusion Schrödinger Bridge Matching (DDSBM), a novel framework that utilizes continuous-time Markov chains to solve the SB problem in a high-dimensional discrete state space.
Our approach extends Iterative Markovian Fitting to discrete domains, and we have proved its convergence to the SB.
Furthermore, we adapt our framework for the graph transformation and show that our design choice of underlying dynamics characterized by independent modifications of nodes and edges can be interpreted as the entropy-regularized version of optimal transport with a cost function described by the graph edit distance.
To demonstrate the effectiveness of our framework, we have applied DDSBM to molecular optimization in the field of chemistry.
Experimental results demonstrate that DDSBM effectively optimizes molecules’ property-of-interest with minimal graph transformation, successfully retaining other features.

## 1 Introduction

Transporting an initial distribution to a target distribution is a foundational concept in modern generative modeling.
Denoising diffusion models (DDMs) have been highly influential in this area, with a primary focus on generating data distributions from simple prior (Sohl-Dickstein et al., ; Song & Ermon, ; Ho et al., ; Song et al., ; Kim et al., ).
Despite their promising results, setting the initial distribution as a simple prior makes DDMs hard to work in tasks where the initial distribution becomes a data distribution, such as image-to-image translation.
To tackle this, diffusion bridge models (DBMs) extend DDMs to transport data between arbitrary distributions (Liu & Wu, ; Liu et al., ; Zhou et al., ).
However, training DBMs requires a coupling between the initial and target distributions, which is often difficult to obtain in practice.
The Schrödinger Bridge (SB) provides an attractive framework for constructing a joint distribution of two data distributions while aligning with the underlying stochastic dynamics (Pariset et al., ; Kim et al., ; Dong et al., ).

Formally, the SB problem seeks the stochastic process that connects two distributions and is closest to a reference process, measured by the Kullback-Leibler (KL) divergence (Schrödinger, ).
The SB problem can be described as an entropy-regularized optimal transport (EOT) problem, which introduces an entropy term to the optimal transport (OT) objective, resulting in randomness in the transport process (Léonard, ).
Here, the transportation cost is determined by the system’s natural dynamics; for example, in the case of Brownian motion, the transportation cost becomes $L^{2}$ (De Bortoli et al., ).
The SB/EOT can be computed efficiently using the Sinkhorn algorithm, though high-dimensional or large-scale data applications remain challenging (Sinkhorn, ; Cuturi, ).
In recent, many methods have been proposed to approximate SB via distribution learning, utilizing techniques developed in DBM and DDM (De Bortoli et al., ; Liu et al., ; Somnath et al., ; Liu et al., ; Peluchetti, ; Shi et al., ; Lee et al., ).

Despite the progress, most of the methods focus on the continuous spaces, where diffusion processes are represented by Brownian motion, and SB problems in discrete domains are less explored.
It is particularly significant in fields that handle discrete state data, such as graphs or natural language (Austin et al., ; Vignac et al., ).
Directly applying frameworks for approximating SB formulated in continuous spaces to these domains limits its potential since it does not reflect the intrinsic properties of discrete data space.
To bridge the gap, we propose a novel framework called Discrete Diffusion Schrödinger Bridge Matching (DDSBM) utilizing the continuous-time Markov chains (CTMCs) to solve the SB problem in a high-dimensional discrete setting.
Our approach leverages Iterative Markovian Fitting (IMF), which was originally proposed for the SB problem in a continuous domain (Peluchetti, ; Shi et al., ).

We then extend our formulation to the graph domain, where the underlying process can be interpreted as independent modifications to both the nodes and edges (Vignac et al., ).
In this case, the cost function of the corresponding EOT can then be regarded as the graph edit distance (GED), which is especially suited for systems where preserving graph similarity is crucial (Sanfeliu & Fu, ).
The molecular optimization in drug/material discovery is such a case in that molecules are represented as graphs.
In addition, the goal is to obtain the molecules with desired molecular properties while retaining favorable properties in acclaimed molecules.
Since molecular structures are closely related to their properties, it is highly advantageous to minimize structural changes (or graph editing) during the optimization.

Figure: Figure 1: A schematic illustration of DDSBM transforming (a) bridge process to (b) Schrödinger Bridge in discrete spaces.
Refer to caption: x1.png

To validate our framework, we evaluated the performance of DDSBM on molecular optimization tasks, with criteria of demonstrating optimal structural modifications to achieve desired property.
DDSBM shows promising results in molecular distribution shift with minimum structural change compared to the previous graph-to-graph translation models.
As a direct result of this, DDSBM retains multiple properties of the initial molecules, along with desired property.
Lastly, we applied DDSBM for a more challenging task, where proper joint pairing between two molecular spaces can not be defined properly.

Our contributions are as follows:

- •
We propose a novel framework, DDSBM, for the SB problem in discrete state spaces by exploiting the IMF procedure and prove its convergence to the SB solution.
- •
We extend our framework to the graph domain, demonstrating a connection between the objective function and the GED.
- •
By reformulating molecular optimization as the SB problem, we show that our approach successfully addresses molecular optimization tasks.

## 2 Theoretical Background

Notations.
We consider the sample space $\Omega=D([0,\tau],\mathcal{X})$, which is set of all left limited and right continuous (càdlàg) paths over the finite metric space $(\mathcal{X},d_{\mathcal{X}})$.
The space of path measures is denoted by $\mathcal{P}(\Omega)=\mathcal{P}(D([0,\tau],\mathcal{X}))$.
The subset of Markov measures of which the transition probability are continuous and differentiable for time is denoted by $\mathcal{M}$.
The transition probability from state $x$ at time $s$ to state $y$ at time $t$ will be represented by $P_{s:t}(x,y)$.
Similarly, the transition rate or generator of a Markov measure as $A_{s}(x,y)$.
For any $M\in\mathcal{M}$, the reciprocal class of $M$ is denoted by $\mathcal{R}(M)$
(see Definition 2.1).
For any $\mathbb{P}\in\mathcal{P}(\Omega)$, $\mathbb{P}_{t}$ denotes the marginal distribution at time $t$, $\mathbb{P}_{s,t}$ as the joint distribution at time $s$ and $t$, $\mathbb{P}_{t|s}$ the conditional distribution at time $t$ given state at time $s$.
$X=(X_{t})_{t\in[0,\tau]}$ denotes the canonical process given by $X_{t}(\omega)=\omega_{t}$ for $\omega=(\omega_{s})_{s\in[0,\tau]}\in\Omega$.
The reference measure $\mathbb{Q}\in\mathcal{M}$ is an irreducible Markov measure on $\Omega$ with its associated canonical filtration.
We will denote the transition rate (or probability) of any Markov measures $M$ other than $\mathbb{Q}$, using a superscript notation.

### 2.1 Schrödinger Bridge Problem

Schrödinger Bridge (SB) problem is finding a stochastic process that most closely resembles a given reference process, with a condition that initial and final marginal distributions are fixed as $\Gamma$ and $\Xi$.
Specifically, the reference process with the law $\mathbb{Q}$ is given, the optimality is achieved by minimizing Kullback-Leibler (KL) divergence to the reference process, $D_{\text{KL}}(\mathbb{P}\|\mathbb{Q})$.
Thus the definition of the SB solution is as below:

$$ $\mathbb{P}^{\text{SB}}=\arg\min_{\mathbb{P}}\{D_{\text{KL}}(\mathbb{P}\| \mathbb{Q}):\mathbb{P}_{0}=\Gamma,\mathbb{P}_{\tau}=\Xi\}.$ (1) $$

If we additionally fix the initial and terminal coupling $\mathbb{P}_{0,\tau}$, the optimality can be found easily as a mixture of bridges $\mathbb{P}_{0,\tau}\mathbb{Q}_{\cdot|0,\tau}$, which implies that finding the SB solution is equivalent to finding optimal coupling $\mathbb{P}_{0,\tau}^{\text{SB}}$ (Léonard, ).
In particular, such optimal coupling is called static SB solution, which could be defined as follows:

$$ $\mathbb{P}_{0,\tau}^{\text{SB}}=\arg\min_{\mathbb{P}_{0,\tau}}\{D_{\text{KL}}( \mathbb{P}_{0,\tau}\|\mathbb{Q}_{0,\tau}):\mathbb{P}_{0}=\Gamma,\mathbb{P}_{ \tau}=\Xi\}.$ (2) $$

Note that the KL-divergence is decomposed into the entropy term $H(\mathbb{P}_{0,\tau})$ and the cross-entropy term $\mathbb{E}_{\mathbb{P}_{0,\tau}}\left[-\log q(x_{\tau}|x_{0})\right]$, where $q$ denotes the density function of $\mathbb{Q}$.
The cross-entropy term is represented as $L^{2}$ distance when the $\mathbb{Q}$ is associated with the reversible Brownian motion.
In general, the static SB problem is equivalent to the entropy-regularized optimal transport (EOT) problem with the cost function $c(x,y)=-\log q_{\tau|0}(y|x)$ (Léonard, ).
While previous approaches have primarily focused on solving SB problems in continuous spaces, typically involving Brownian motion, our work differs by addressing SB solutions in discrete spaces, specifically dealing with cádlág paths.
Despite this distinction, the core idea remains rooted in the intrinsic properties of the SB solution (see Theorem A.1).

### 2.2 Solution Method

To solve the SB problem, we adopted the Iterative Markovian Fitting (IMF) method, a technique that has been proposed in various studies (Peluchetti, ; Shi et al., ).
These works primarily address diffusion processes in Euclidean space, whereas our approach differs slightly as we focus on a continuous-time Markov chain over finite state space.
Applying IMF to càdlàg paths in discrete spaces requires additional theoretical developments or modifications to ensure proper convergence and stability in this new context, which is described in Appendix A.

The SB solution is a mixture of pinned-down measures of $\mathbb{Q}(\cdot|X_{0}=x_{0},X_{\tau}=x_{\tau})$, where the pair $(x_{0},x_{\tau})$ is drawn from the coupling $\mathbb{P}_{0,\tau}^{\text{SB}}$.
Based on this, the projection method first constructs a reciprocal measure, which is mixture of pinned-down processes from a given initial coupling.
Although each pinned-down process is Markov, the mixture generally loses the Markov property in general as collection of Markov process is non-convex (Léonard et al., ).
Thus, it then identifies the Markov measure that is closest to the mixture.
This yields an improved coupling, and the process is iterated to obtain a measure that is both the mixture of pinned-down measures and Markov– the desired SB solution.

In this context, we define the reciprocal projection to describe the construction of a reciprocal mixture from a given coupling (see Definition 2.1).
Similarly, the term Markov projection is used to describe the approximation of a reciprocal process with a Markov process (see Definition 2.2).

###### Definition 2.1 .

(Reciprocal Projection)
$\Lambda\in\mathcal{P}(\Omega)$ is in the reciprocal class $\mathcal{R}(\mathbb{Q})$ with respect to a Markov measure $\mathbb{Q}$ if $\Lambda=\Lambda_{0,\tau}\mathbb{Q}_{|0,\tau}$.
For a measure $\Lambda\in\mathcal{P}(\Omega)$, its reciprocal projection with respect to $\mathbb{Q}$ is

$$ $\Pi_{\mathcal{R}(\mathbb{Q})}(\Lambda)(\cdot):=\iint_{(\cdot)}\Lambda(dx_{0}, dx_{\tau})\mathbb{Q}(dx_{t}|x_{0},x_{\tau}).$ $$

Among the measures with the coupling $\Lambda_{0,\tau}\neq\mathbb{P}^{\text{SB}}_{0,\tau}$, the minimizer of the KL-divergence to the reference process is the (non-Markov) reciprocal projection $\Pi_{\mathcal{R}}(\Lambda)$.
The reciprocal projection consists of a mixture of bridges, where each bridge is derived from Doob’s h-transform with the realization of an end-point pair $(x,y)\sim\Lambda_{0,\tau}$ (Levin & Peres, ).
Obviously, it preserves the initial-terminal coupling.
Although each pinned-down bridge is Markov (see Lemma A.4), the mixture is generally not Markov.

###### Definition 2.2 .

(Markov Projection)
Given a path measure $\Lambda\in\mathcal{R}(\mathbb{Q})$, a Markov path measure that minimizes the reverse KL-divergence to $\Lambda$ is called as Markov projection of $\Lambda$,

$$ $\Pi_{\mathcal{M}}(\Lambda)=\arg\min_{M}\left\{D_{\text{KL}}(\Lambda\|M):M\in \mathcal{M}\right\}.$ $$

The Markov projection preserves the marginal distribution at all times $t$, but does not preserve the coupling.
Furthermore, the generator and the reverse KL-divergence of the Markov projection over cádlág paths are explicitly derived in the Proposition A.2.

For a given reciprocal process $\Lambda^{(0)}\in\mathcal{R}(\mathbb{Q})$, we consider a sequence $(\Lambda^{(n)})_{n\in\mathbb{N}}$ which is defined by the recurrence relation:

$$ $\displaystyle\Lambda^{(2n+1)}$ $\displaystyle=\Pi_{\mathcal{M}}(\Lambda^{(2n)}),$ (3) $\displaystyle\Lambda^{(2n+2)}$ $\displaystyle=\Pi_{\mathcal{R}(\mathbb{Q})}(\Lambda^{(2n+1)}),$ $$

where $\Lambda_{0}^{(0)}=\Gamma$ and $\Lambda_{\tau}^{(0)}=\Xi$.
Under mild assumptions, the resulting sequence of measures of iterative projection converges to $\mathbb{P}^{\text{SB}}$ in law (see Theorems 2.3 and A.5).

###### Theorem 2.3 .

(Convergence of Iteration)
Assume that $D_{\text{KL}}(\Lambda^{(0)}_{0,\tau}\|\mathbb{P}^{\text{SB}}_{0,\tau})<\infty$, $\Lambda^{(n)}\ll\mathbb{P}^{\text{SB}}$ for all $n\in\mathbb{N}$.
Then the sequence of KL-divergence to $\mathbb{P}^{\text{SB}}$ non-increasing,

$$ $D_{\text{KL}}(\Lambda^{(2n)}\|\mathbb{P}^{\text{SB}})\geq D_{\text{KL}}( \Lambda^{(2n+1)}\|\mathbb{P}^{\text{SB}})\geq D_{\text{KL}}(\Lambda^{(2n+2)}\| \mathbb{P}^{\text{SB}}).$ $$

$D_{\text{KL}}(\Lambda^{(2n)}\|\mathbb{P}^{\text{SB}})=D_{\text{KL}}(\Lambda^{(
2n+1)}\|\mathbb{P}^{\text{SB}})$ if and only if $\Lambda^{(2n)}=\Lambda^{(2n+1)}=\mathbb{P}^{\text{SB}}$.
Moreover, $\Lambda^{(n)}$ converges to $\mathbb{P}^{\text{SB}}$ in law as $n\to\infty$.

###### Definition 2.1 .

###### Definition 2.2 .

###### Theorem 2.3 .

## 3 Methods

Here, we propose the Discrete Diffusion Schrödinger Bridge Matching (DDSBM) framework, which focuses on solving the SB problem in discrete state spaces.
Approaches to the SB problem in continuous spaces model probability flows through stochastic differential equations, while our method uses continuous-time Markov chains (CTMCs) in discrete state spaces, governed by the Kolmogorov equation Equation 9.
We first provide the algorithm of DDSBM (Section 3.1) to apply Iterative Markovian Fitting (IMF) to càdlàg paths in discrete spaces, of which the convergence is ensured theoretically by Theorem 2.3.
We then discuss how DDSBM can be implemented in the graph transformation (Section 3.2) and introduce a graph permutation matching algorithm for efficiency (Section 3.3).

### 3.1 Algorithm for Solving the Schrödinger Bridge Problem

In this section, we present the framework for solving the SB problem in discrete state spaces based on the IMF algorithm proposed by Shi et al. ().
We refer to it as the DDSBM framework.
The framework iteratively applies Markov and reciprocal projections to update the sequence of measures, as described in Equation 3.

It begins with an initial random coupling $\pi$ such that $\pi_{0}=\Gamma$ and $\pi_{\tau}=\Xi$.
Following the definition of the reciprocal class in Definition 2.1, we construct the initial reciprocal bridge to obtain the measure $\Lambda^{(0)}$.

Given a reciprocal measure $\Lambda^{(2n)}\in\mathcal{R}(\mathbb{Q})$, the next step is to compute its Markov projection $M^{(2n+1)}:=\Pi_{\mathcal{M}}(\Lambda^{(2n)})$.
The exact form of the transition rate for $M^{(2n+1)}$ is provided in Proposition A.2.
In practice, the transition rate is approximated by a neural network.

To achieve this, it first samples pairs $(x_{0},x_{\tau})$ from $\Lambda_{0,\tau}^{(2n)}$ and samples intermediate states $x_{t}$ by constructing the bridge $\mathbb{Q}(\cdot|X_{0},X_{\tau})$.
The sampled pairs $(x_{t},x_{\tau})$ are distributed according to $\Lambda_{t,\tau}^{(2n)}$.
Using these realizations, the rate matrix of $M^{(2n+1)}$ is approximated by parameterized Markov measure $M^{\theta}$, by minimizing the following loss function:

$$ $\mathcal{L}(\theta)=\int_{0}^{\tau}\mathbb{E}_{\Lambda_{t,\tau}^{(2n)}}\left[( A_{t}^{\mathbb{Q}_{\cdot|\tau}}-A_{t}^{M^{\theta}})(X_{t},X_{t})+\sum_{y\neq X _{t}}A_{t}^{\mathbb{Q}_{\cdot|\tau}}\log{\frac{A_{t}^{\mathbb{Q}_{\cdot|\tau}} }{A_{t}^{M^{\theta}}}}(X_{t},y)\right]dt,$ (4) $$

where $A_{t}^{\mathbb{Q}_{\cdot|\tau}}$ denotes the generator of pinned down process of $\mathbb{Q}(\cdot|X_{\tau})$.
From the approximated generator $A_{t}^{M^{\theta}}$, we can sample $x_{\tau}$ given $x_{0}$, where $(x_{0},x_{\tau})\sim M^{\theta}_{0,\tau}\approx M^{(2n+1)}_{0,\tau}$.
Note that, until the sequence of measures converges, the new joint coupling $M^{\theta}_{0,\tau}\approx M^{(2n+1)}_{0,\tau}$ will differ from the previous one $\Lambda_{0,\tau}^{(2n)}$.

Once the Markov measure $M^{(2n+1)}$ is obtained, we proceed to compute the corresponding reciprocal measure $\Lambda^{(2n+1)}$ through reciprocal projection, $\Lambda^{(2n+1)}:=\Pi_{\mathcal{R}(\mathbb{Q})}(M^{(2n+1)})$.
In theory, as shown in Proposition A.2, the time marginal distributions are preserved under Markov projection, meaning that $\Lambda^{(2n)}_{t}=M^{(2n+1)}_{t}$ for all $t\in[0,\tau]$.
However, in practice, since the Markov projection is approximated using a neural network, repeatedly applying this approximation can lead to an accumulation of errors in the time marginals.
Such accumulated errors may violate the marginal condition of the SB problem, particularly leading to a potential failure in satisfying the terminal condition $\mathbb{P}_{\tau}=\Xi$.

To compensate these errors, the next Markov measure $M^{(2n+2)}:=\Pi_{\mathcal{M}}(\Lambda^{(2n+1)})$ is approximated in a time-reversal way (see Proposition A.8).
Based on the time-symmetric nature of Markov measures, we can leverage the time-reversed generator $\tilde{A}_{t}^{M^{(2n+2)}}$, which enables the sampling of $x_{0}$ conditioned on $x_{\tau}$, where $(x_{0},x_{\tau})\sim M^{(2n+2)}_{0,\tau}$.
The approximation of $\tilde{A}_{t}^{M^{(2n+2)}}$ is achieved by minimizing the following loss function:

$$ $\mathcal{L}(\phi)=\int_{0}^{\tau}\mathbb{E}_{\Lambda_{0,t}}\left[(\tilde{A}_{t }^{\mathbb{Q}_{\cdot|0}}-\tilde{A}_{t}^{M^{\phi}})(X_{t},X_{t})+\sum_{y\neq X_ {t}}\tilde{A}_{t}^{\mathbb{Q}_{\cdot|0}}\log{\frac{\tilde{A}_{t}^{\mathbb{Q}_{ \cdot|0}}}{\tilde{A}_{t}^{M^{\phi}}}}(X_{t},y)\right]dt,$ (5) $$

where $\tilde{A}_{t}^{\mathbb{Q}_{\cdot|0}}$ denotes the time-reversal generator of pinned down process of $\mathbb{Q}(\cdot|X_{0})$, and $\phi$ represents the parameters of the neural network approximating the time-reversed generator $\tilde{A}_{t}^{M^{\phi}}$.

In this manner, the iterative Markov projection following the reciprocal projection is performed alternately in a forward and backward (time-reversal) fashion.
Finally, this process yields a sequence of measures $(\Lambda^{(n)})_{n\in\mathbb{N}}$ and $(M^{(n)})_{n\in\mathbb{N}^{+}}$, which converges to $\mathbb{P}^{\text{SB}}$ in theory.

### 3.2 Process on the Graph

We present a method for applying the previously described solution to a graph transformation.
A graph is represented as $\mathbf{G}=(\mathbf{V},\mathbf{E})$, where $\mathbf{V}=(v^{(i)})_{i}$ and $\mathbf{E}=(e^{(ij)})_{ij}$ denote node and edge features, respectively.
In a molecular graph, for example, the node and edge features correspond to atomic types and edge features, respectively.
Here, $\mathbf{V}$ and $\mathbf{E}$ are modeled as products of categorical random variables.

As the reference process, we define a jump process in which the nodes and edges of the graph change discretely, assuming that all nodes and edges are independent.
Therefore, the transition probability $P$ and the rate $A$ of the reference process is described as follows:

$$ $\displaystyle P_{s:t}^{\mathbf{G}}(\mathbf{G}_{1},\mathbf{G}_{2})=\prod_{i}P_{ s:t}^{V}(v_{1}^{(i)},v_{2}^{(i)})\prod_{i,j}P_{s:t}^{E}(e_{1}^{(ij)},e_{2}^{( ij)}),$ $\displaystyle\partial_{t}P_{s:t}^{(\cdot)}(x,y)=\sum_{z\in\mathcal{X^{(\cdot)} }}P_{s:t}^{(\cdot)}(x,z)A_{t}^{(\cdot)}(z,y),$ (6) $\displaystyle A^{(\cdot)}_{s}(x,y)=\partial_{t}P_{s:t}^{(\cdot)}(x,y)\big{|}_{ t=s},$ $$

where $\cdot$ denotes $V$ or $E$, $\mathbf{G}_{1}$ denotes $(\mathbf{V}_{1},\mathbf{E}_{1})=\left((v_{1}^{(i)}),(e_{1}^{(ij)})\right)$, $\mathbf{G}_{2}$ denotes $(\mathbf{V}_{2},\mathbf{E}_{2})=\left((v_{2}^{(i)}),(e_{2}^{(ij)})\right)$, and $\mathcal{X}^{V}$ and $\mathcal{X}^{E}$ denote the state space of the nodes and edges, respectively.
More specifically, we use a monotonically decreasing function for signal to noise ratio, $\bar{\alpha}:[0,\tau]\to(0,1]$, in which the transition rate is defined as:

$$ $A_{t}^{(\cdot)}(x,y)=\partial_{t}(\ln{\bar{\alpha}}(t))\left(\delta_{xy}- \mathbf{m}^{(\cdot)}(y)\right),$ (7) $$

where $\delta_{xy}$ denotes the Kronecker delta, and $\mathbf{m}^{(V)}$ and $\mathbf{m}^{(E)}$ denote the prior distribution of nodes and edges as proposed in the previous discrete diffusion work (Vignac et al., ).
According to the Kolmogorov equation, we get the transition probability as,

$$ $P_{s:t}^{(\cdot)}(x,y)=\frac{\bar{\alpha}(t)}{\bar{\alpha}(s)}\delta_{xy}+ \left(1-\frac{\bar{\alpha}(t)}{\bar{\alpha}(s)}\right)\mathbf{m}^{(\cdot)}(y).$ (8) $$

Note that the choice of $\mathbf{m}(\cdot)$ as uniform is associated to the diffusion on the state space $\mathcal{X}$, where the $\mathcal{X}$ is considered fully-connected graph.
Moreover, the stationary distribution of the associated generator always becomes $\mathbf{m}(\cdot)$.
Although non-uniform choice of $\mathbf{m}(\cdot)$ breaks the diffusion formulation on $\mathcal{X}$, it does not harm SB formulation.

### 3.3 Graph Permutation Matching

The transition probability of the reference process depends on graph permutations (see Section 3.2), so graph permutation matching must be considered beforehand.
Although this issue does not affect the sampling phase, it becomes problematic when computing the likelihood of the reference process for two given graphs $\mathcal{G}$ and $\mathcal{G}^{\prime}$, or when constructing a reciprocal process, which is a mixture of Markov bridges between the two graphs (see Appendix B).

One way to handle this is selecting a graph permutation that maximizes the likelihood under the reference process.
Finding the optimal permutation can be formulated as a quadratic assignment problem (QAP), where the objective is to minimize the negative log-likelihood (NLL), consisting of the sum of the NLLs for both the nodes and edges.
While exact solutions are possible through mixed integer programming, the problem is NP-hard, so alternative methods are preferred.
Specifically, we employ a max-pooling algorithm by Cho et al. (), which is an approximation method categorized by continuous relaxation.
After obtaining an approximate solution, we use the Hungarian algorithm to discretize the assignment vector to the final solution.
We observed that the graph matching is highly accurate in molecular graph matching (see Section B.5).
The details of the algorithm are described in Section B.4.
We utilized the graph matching algorithm for every reciprocal bridge construction and likelihood computation.

Recall that the SB problem could be interpreted to the EOT problem, where the cost function corresponds to the NLL.
Thus, defining the reference process can be interpreted as defining a distance (cost) on the set of graphs.
Interestingly, we found that the likelihood of optimal permutation is interpreted as the graph edit distance (see Section B.6).
This leads to the conclusion that the SB problem, with the $\mathbb{Q}$ as Equation 8, is analogous to an OT problem over the metric space of graphs equipped with graph edit distance as metric (GED), where the GED computation is well known to be NP-hard problem.

## 4 Related Works

Diffusion Bridge Models.
Diffusion bridge models (DBMs) have recently shown state-of-the-art results in a variety of continuous domains, such as images, biology, and chemistry (Liu & Wu, ; Liu et al., ; Somnath et al., ; Zhou et al., ; Lee et al., ).
While Igashov et al. (); Yang et al. () have extended these models to discrete domains, they focused on settings where well-defined data pairs exist.
To the best of our knowledge, we are the first to study DBM in discrete domains where proper joint distributions between data points are not provided or well-defined.

Schrödinger Bridge Problem.
The Schrödinger Bridge (SB) problem is an important concept in recent generative modeling (Liu et al., ; Somnath et al., ; Peluchetti, ; Shi et al., ).
In particular, incorporating the SB problem into DBM can address scenarios where no appropriate joint distribution is available, as demonstrated by recent works (Pariset et al., ; Kim et al., ; Dong et al., ).
For example, Somnath et al. () proposed learning an SB based on an assumed (partial) true coupling, while Shi et al. () showed that it is possible to generate high-quality samples from arbitrary couplings that are well-aligned with the initial data.
However, most existing approaches focus on continuous spaces.
To the best of our knowledge, we are the first to propose a framework for solving the SB problem in high-dimensional discrete spaces.

Molecular Optimization.
Molecular optimization is a promising strategy in drug/material discovery that aims to improve acclaimed molecules to satisfy multiple domain-specific properties.
One major approach to the molecular optimization problem is to formulate it as a graph transformation problem, which can be categorized into latent-based and graph-editing approach.
The latent based approaches such as JT-VAE by Jin et al. () and HierG2G by Jin et al. () encode an input graph into a single latent vector and then decode it into a whole graph that follows a certain data distribution.
The graph-editing approaches such as Modof by Chen et al. () and DeepBioisostere by Kim et al. () learn a graph editing procedure to transport a given molecule to another.
Despite their promising results in optimizing molecular graphs, their training schemes rely on paired data created by predefined rules, which limits not only the general applicability but also the optimality of the transformations.
In this work, we propose a more flexible framework for molecular optimization by formulating it as a graph SB problem, leading to more optimal graph transformation accompanying less structural changes.

## 5 Results and Discussions

Here, we demonstrate the effectiveness of the Discrete Diffusion Schrödinger Bridge Matching (DDSBM) framework on graph transformation tasks.
Specifically, we apply DDSBM to a chemical domain, where the graph transformation task is nontrivial due to additional constraints imposed by molecular graphs and their associated properties.
We conducted experiments on two different molecule datasets: ZINC250K (Kusner et al., ) and Polymer (St. John et al., ).
Throughout this section, we first elaborate on the common experimental setups and metrics for evaluation.
The second and third sections provide a detailed analysis of ZINC250K and Polymer experiments.
Finally, we empirically show utilizing an initial coupling based on the Tanimoto similarity (Willett et al., ), the widely-used metric for assessing chemical similarity, is effective for training DDSBM.

### 5.1 Experimental setup and metrics

Experimental Setup. To train the models on graph transformation problems, an initial coupling between two distributions is necessary.
We randomly coupled the data of initial and terminal distributions to obtain paired data.
The molecular pairs are divided into training and test datasets in a ratio of 8:2.
The effect of different initial couplings is discussed in Section 5.4.
Detailed explanations about model architectures and hyperparameters can be found in Appendix C.

Metrics.
By definition of SB (see Equation 1), both joint distribution and marginal distributions at each side must be examined to assess the degree to which the SB has been successfully achieved.
Specifically, we evaluate these distributions from two perspectives: graph structure and molecular properties.
While graph structure provides valuable insights, it alone may not be sufficient for molecular optimization tasks.
Therefore, the assessment of molecular properties is essential for a more comprehensive validation of the models, as these properties might exhibit weak correlations with the graph structures.

In the context of graph structure, we calculate Fréchet ChemNet Distance (FCD) score (Preuer et al., ) between the target molecule distribution and the generated molecules from each model, which focuses on the marginal distribution of the molecules.
As the metric for evaluating the quality of the joint distribution, we report the negative log-likelihood (NLL), which measures the cost of transforming one graph into another based on the reference process, aligning with the graph edit distance (GED) (see Section B.6).

For molecular properties, we measure the Wasserstein-1 distance ($W_{1}$) between the property distributions of the target and generated molecules representing the properties to be modified, which corresponds to the marginal distribution.
In terms of the joint distribution, we calculate the mean absolute difference (MAD) of other key properties between the initial and generated molecules that should be retained, where drug-likeness (QED) (Bickerton et al., ) and synthetic accessibility score (SAscore) (Ertl & Schuffenhauer, ) are typical examples in optimizing drug candidates.

### 5.2 Small molecule transformation

First, we validate our proposed methods on the SB problem between two molecule distributions constructed from the standard ZINC250K dataset.
We constructed two sets of molecules whose $\log{P}$ values are largely different.
Molecules from the ZINC250K dataset were randomly selected and divided into two sets whose $\log{P}$ values follow the Gaussian distributions centered at 2 and 4 with variance of 0.5, respectively (For more details, see Section D.1).
We compared our methods with three baseline models that perform graph-to-graph translation.
AtomG2G and HierG2G are latent-based models that encode an input graph into a latent vector and decode it into another graph (Jin et al., ).
Diffusion Bridge Model (DBM) is a bridge model trained with the same reference process as DDSBM, which is equivalent to the first Markov projection in DDSBM.
We refer readers to Section C.2 for more details about the baselines and implementation of our models.
We note that, for these three baseline models, the initial coupling is utilized during the whole training procedure without change, while DDSBM dynamically alters the training data pair by Iterative Markovian Fitting (IMF).

**Table 1: Distribution shift performance on ZINC. Reference refers to metrics from the initial coupling, used as a standard to evaluate each model’s graph translation. For both AtomG2G and HierG2G, we’ve excluded the generated molecules that are too large with more atoms than the maximum number of atoms in our dataset for computing metrics other than validity. $\uparrow$ and $\downarrow$ denote higher and lower values are better, respectively. The best performance is highlighted in bold.**
| Model | Type | Validity$(\uparrow)$ | FCD$(\downarrow)$ | NLL$(\downarrow)$ | LogP $W_{1}(\downarrow)$ | QED MAD$(\downarrow)$ | SAscore MAD$(\downarrow)$ |
| --- | --- | --- | --- | --- | --- | --- | --- |
| Reference^1 | - | - | 4.807 / 0.279 | 360.862 | 2.007 | 0.153 | 0.595 |
| AtomG2G | Latent | 100.0 | 5.480 | 347.280 | 0.222 | 0.135 | 0.742 |
| HierG2G | Latent | 100.0 | 4.413 | 344.208 | 0.053 | 0.135 | 0.657 |
| DBM | Bridge | 89.5 | 1.038 | 285.576 | 0.161 | 0.145 | 0.611 |
| DDSBM | Schrödinger Bridge | 95.3 | 0.896 | 162.228 | 0.134 | 0.118 | 0.398 |

Table 1 shows overall results of our method compared to the three baseline models.
From Table 1, we find that DDSBM constantly outperforms the other models in terms of both FCD and NLL.
The lower FCD suggests that DDSBM-generated molecules are more similar to those in the target dataset, while the minimal NLL indicates that DDSBM applies minimal structural change on initial graphs.
This result demonstrates that DDSBM achieves more optimal graph transformation between fixed initial and terminal distributions.
When it comes to the model type, bridge models achieved lower FCD and NLL compared to the latent-based models.
The bridge models, DBM and DDSBM, dynamically transform a graph based on the reference dynamics, favoring the retention of the given structure, whereas whole-graph reconstruction using a latent vector does not.
This leads to lower NLL values for the bridge models, which is ensured by the Definition 2.2.
Additionally, the lower FCD values of the bridge models highlight that constructing a dynamic bridge within the graph domain enhances the performance of distribution learning for the target distribution.
Furthermore, we attribute the best performance of DDSBM to the gradual updating of training pairs, which become more similar than the random initial pairs, simplifying the graph transportation process.

Next, we analyze molecular properties of the source and generated molecules.
DDSBM resulted in the lowest MAD in QED and SAscore, meaning deviation of molecular properties other than $\log{P}$ is the smallest for DDSBM.
It is noteworthy that DDSBM achieves minimal changes in various molecular properties despite being trained solely to optimize graph transformations with minimal structural alterations, without explicit knowledge on molecular properties.
Meanwhile, the latent-based models exhibited larger MAD values in QED and SAscore, indicating they are vulnerable to losing other properties of the initial molecules.
This can be inferred from the result of HierG2G; although HierG2G achieves the lowest $W_{1}$ value in $\log{P}$ so that it modulates $\log{P}$ the best to the desired degree, its much larger NLL value suggests that it might generate a graph with less consideration of the initial graph constraint, as illustrated in Figure 4.
For a better understanding of overall results, we provide the corresponding distributions of properties of the models in Figure 6.

Despite the promising results, we observe that all models except DDSBM have a high reliance on a pre-defined initial coupling.
Thus, we conducted additional experiments using pseudo-optimal initial coupling based on Tanimoto similarity, which is detailed in Section 5.4.
Interestingly, we found that introducing the pseudo-optimal initial coupling not only accelerated the training of DDSBM in practice but also allowed it to achieve performance on par with the previous results.

### 5.3 Polymer graph transformation

The Polymer dataset (St. John et al., ) consists of 91,000 monomer molecules with their optical excitation energies (GAPs) obtained by time-dependent density functional theory calculations.
We reconstructed the Polymer dataset for a transport problem between two sets of molecules with distinct GAPs, corresponding to green and blue optical colors, respectively.
The reconstructed dataset contains 7,603 molecular pairs.
This task is considered as a more challenging application because the relationship between the graph structure and the target GAP property is highly non-linear, making it hard to predict the effect of specific structural changes on the GAP.
In this context, we apply our DDSBM model to find the optimal transformation between the two sets of molecules.

The performance of DDSBM is compared to DBM, which serves as a baseline.
The GAPs of the molecules generated by the models were obtained using the pre-trained MolCLR model (Wang et al., ), which has a mean absolute error (MAE) of 0.027 eV for the GAP prediction.
Table 2 shows the overall results of our method on the Polymer dataset.
Apparently, DDSBM outperforms DBM on all evaluated metrics, which is consistent with the results from the experiments on ZINC250K.
In terms of validity, DBM shows significantly lower scores, which contrasts with the results observed on the ZINC250K dataset.
This can be attributed to the characteristics of the molecules in the Polymer dataset, which have relatively large sizes and multiple ring structures.
In this context, minimal transformations are advantageous for achieving high validity, and DBM may have struggled to learn these changes from the randomly paired data.
Examples of the generated molecular graphs are visualized in Figure 5.

**Table 2: Distribution shift performance on Polymer. Reference refers to metrics from the initial coupling, used as a standard to evaluate each model’s graph translation. $\uparrow$ and $\downarrow$ denote higher and lower values are better, respectively. The best performance is highlighted in bold.**
| Model | Type | Validity$(\uparrow)$ | FCD$(\downarrow)$ | NLL$(\downarrow)$ | GAP $W_{1}(\downarrow)$ |
| --- | --- | --- | --- | --- | --- |
| Reference^1 | - | - | 1.469 / 0.384 | 749.800 | 0.312 |
| DBM | Bridge | 35.7 | 2.803 | 585.043 | 0.258 |
| DDSBM | Schrödinger Bridge | 96.6 | 1.101 | 215.590 | 0.136 |

### 5.4 The effect of initial coupling

We further analyze the molecular optimization task in Section 5.2 with different initial couplings.
In that section, we used randomly coupled molecules as the initial coupling for all the models.
However, except for DDSBM, all of them assume that suitable molecule pairs have already been identified to be provided.
To meet this, we adopted Tanimoto similarity (Willett et al., ) as a pseudo-metric to find similar molecule pairs between the two molecule distributions, denoted as Tanimoto similarity-based coupling.
We note that the similarity-based coupling is another optimal transport problem of maximizing the sum of pair-wise molecular similarities, where we employed the Hungarian method to obtain a sub-optimal solution.

Using the Tanimoto similarity-based coupling, we retrained all models discussed in Section 5.2 and compared their performance.
Here, we denote a newly trained DDSBM model as DDSBM-T to deviate it from the model trained on the randomly coupled data.
From Table 3, we observe that all models achieved lower NLL values compared to when they were trained with random coupling.
The HierG2G model exhibits much lower FCD and NLL values compared to those of random coupling, indicating that previous graph transformation methods can be improved if more optimal pairs are provided as training data.

The distinct feature of DDSBM-T is that all baseline models learn graph transformations between the molecule pairs with high Tanimoto similarity, while DDSBM-T learns graph transformations with minimal cost defined by its reference process $\mathbb{Q}$.

Figure: Figure 2: Comparison of $D_{\textrm{KL}}(\Lambda|M^{\theta})$ across IMF iterations for two types of initial couplings: Random and Tanimoto similarity.
Refer to caption: extracted/5895698/Figures/source/KL_Loss_Convergence_0.999_uniform.png

It cannot be guaranteed that the distance defined by our reference process is better than the Tanimoto similarity from the perspective of molecular optimization.
Essentially, the success of molecular optimization should be measured by how well the target property is adjusted while preserving other key properties.
Apparently, Table 3 shows that DDSBM-T outperforms the other baselines in terms of molecular property metrics.
This suggests that the graph transformation from DDSBM retains other molecular properties, attaining the goal of the molecule optimization task.

Finally, we analyze the effect of initial coupling on the training process, especially focusing on the approach to convergence.
We illustrate the training losses of DDSBM models from two different initial couplings, random and Tanimoto similarity-based, in Figure 2, respectively.
DDSBM-T shows consistently lower loss values up to the sixth IMF iteration and reaches convergence at the third IMF iteration.

**Table 3: Distribution shift performance on ZINC with initial coupling based on the Tanimoto similarity. As in Table 1, reference refers to metrics from the initial coupling, used as a standard to evaluate each model’s graph translation. The experimental setting is the same as described Section 5.2, except for the initial coupling.**
| Model | Type | Validity$(\uparrow)$ | FCD$(\downarrow)$ | NLL$(\downarrow)$ | LogP $W_{1}(\downarrow)$ | QED MAD$(\downarrow)$ | SAscore MAD$(\downarrow)$ |
| --- | --- | --- | --- | --- | --- | --- | --- |
| Reference^1 | - | - | 4.811 / 0.315 | 245.765 | 2.011 | 0.126 | 0.367 |
| AtomG2G | Latent | 100.0 | 5.676 | 346.499 | 0.220 | 0.136 | 0.738 |
| HierG2G | Latent | 100.0 | 1.171 | 264.072 | 0.189 | 0.127 | 0.446 |
| DBM | Bridge | 90.2 | 0.749 | 220.594 | 0.141 | 0.127 | 0.508 |
| DDSBM-T | Schrödinger Bridge | 95.6 | 0.911 | 152.856 | 0.103 | 0.110 | 0.393 |

## 6 Discussion

In this paper, we propose Discrete Diffusion Schrödinger Bridge Matching (DDSBM), a novel framework utilizing continuous-time Markov chains to solve the SB problem in a high-dimensional discrete state space.
To this end, we extend Iterative Markovian Fitting (IMF), proving its convergence to SB.
We successfully apply our framework to graph transformation, specifically for molecular optimization.
Experimental results demonstrate that DDSBM effectively transforms molecules to achieve the desired property with minimal graph transformation, while retaining other features.

However, our method has several limitations.
The IMF requires iterative sampling from the learned Markov process, which can be more computationally intensive than simple bridge matching.
Additionally, the graph permutation matching may introduce challenges for wider applicability in general graph tasks.

## Appendix A Propositions and proof

### A.1 Notations

In this section, we introduce the notations that will be used throughout the proofs of the propositions.
$\Omega=D([0,\tau],\mathcal{X})$ denotes the space of all left-limited and right-continuous (càdlàg) paths over a finite space $\mathcal{X}$.
We assume the state space $\mathcal{X}$ has connected finite graph structure (fully-connected graph), which imply it becomes metric space with graph distance metric $d_{\mathcal{X}}(\cdot,\cdot)$.
Accordingly, the sample space $\Omega$ is equipped with Skorokhod topology with Skorokhod metric $d_{\Omega}(\cdot,\cdot)$, and associated Borel $\sigma$-algebra.
$X=(X_{t})_{t\in[0,\tau]}$ denotes the canonical process given by:

$$ $X_{t}(\omega)=\omega_{t},\quad t\in[0,\tau],\quad\omega=(\omega_{s})_{s\in[0, \tau]}.$ $$

The reference measure $\mathbb{Q}$ is an irreducible Markov measure on $\Omega$ with its associated canonical filtration.
Assume that the transition probability of $\mathbb{Q}$, denoted by $P_{s:t}(x,y)$, from $(s,x)\in[0,\tau]\times\mathcal{X}$ to $(t,y)\in[0,\tau]\times\mathcal{X}$, is continuous and differentiable over time.
The measure is generated by the transition rate function $A_{s}(x,y)$, which gives the rate of transition from $x\in\mathcal{X}$ to $y\in\mathcal{X}$ at time $s\in(0,\tau)$, and satisfies the Kolmogorov forward equation:

$$ $\displaystyle\frac{\partial P_{s:t}(x,y)}{\partial t}$ $\displaystyle=\sum_{z\in\mathcal{X}}P_{s:t}(x,z)A_{t}(z,y),$ (9) $\displaystyle A_{s}(x,y)$ $\displaystyle=\left[\frac{\partial P_{s:t}(x,y)}{\partial t}\right]_{t=s}.$ $$

We also assume that $\mathbb{Q}$ can construct a bridge $\mathbb{Q}(\cdot|X_{0}=x,X_{\tau}=y)$ for all $x,y\in\mathcal{X}$.
For any Markov measure $M$, the corresponding generator is denoted as $A^{(M)}$.

### A.2 Theorem A.1

###### Theorem A.1 .

(Uniqueness of the Schrödinger Bridge Solution)
If the reference process $\mathbb{Q}$ is Markov, then under mild conditions the Schrödinger Bridge solution $\mathbb{P}^{\text{SB}}$ exists and is unique.
Furthermore, the solution is mixture with static Schrödinger Bridge solution represented as $\mathbb{P}^{\text{SB}}=\mathbb{P}^{\text{SB}}_{0,\tau}\mathbb{Q}_{\cdot|0,\tau
}\in\mathcal{R}(\mathbb{Q})$, and also it is in $\mathcal{M}$.
Conversely, a process $\mathbb{P}=\mathbb{P}_{0,\tau}\mathbb{Q}_{\cdot|0,\tau}$ is Markov if and only if $\mathbb{P}=\mathbb{P}^{\text{SB}}$.

###### Proof.

This is direct consequence of Theorem 2.12 of (Léonard, ).
∎

###### Theorem A.1 .

###### Proof.

### A.3 Markov projection

###### Proposition A.2 .

(solution of Markov projection)
Let $M^{*}=\Pi_{\mathcal{M}}(\Lambda)$, $\Lambda\in\mathcal{R}(\mathbb{Q})$ and the generator of $\mathbb{Q}$ be $A_{t}(x,y)$ with transition probability $P_{s:t}(x,y)$.
Under mild assumptions, the generator of the Markov measure $M^{*}$ becomes

$$ $\displaystyle A^{(M^{*})}_{t}(X_{t},y)$ $\displaystyle=\mathbb{E}_{\Lambda_{\tau|t}}\left[A_{t}(X_{t},y;X_{\tau})\bigg{ |}X_{t}\right]$ $\displaystyle A_{t}(x,y;z)$ $\displaystyle=A_{t}(x,y)\frac{P_{s:\tau}(y,z)}{P_{s:\tau}(x,z)}-\delta_{xy} \sum_{u}A_{t}(y,u)\frac{P_{t:\tau}(u,z)}{P_{t:\tau}(x,z)}$ $$

, where $z\in\mathcal{X}$

The reverse KL-divergence is

$$ $\displaystyle D_{\text{KL}}(\Lambda\|M^{*})$ $\displaystyle=\int_{0}^{\tau}\mathbb{E}_{\Lambda_{0,t}}\left[(A^{(\Lambda_{|0} )}_{t}-A^{(M^{*})}_{t})(X_{t},X_{t})+\sum_{y\neq X_{t}}A^{(\Lambda_{|0})}_{t} \log\frac{A^{(\Lambda_{|0})}_{t}}{A^{(M^{*})}_{t}}(X_{t},y)\right]dt,$ $$

where the $A^{\Lambda_{|0}}$ is the generator for the conditioned measure $\Lambda_{|0}$ which is defined as

$$ $A^{\Lambda_{|0}}_{t}(X_{t},y)=\mathbb{E}_{\Lambda_{\tau|0,t}}\left[A_{t}(X_{t} ,y;X_{\tau})\bigg{|}X_{t},X_{0}\right].$ $$

Moreover, for any $t\in[0,\tau],\Lambda_{t}=M^{*}_{t}$.

###### Proposition A.2 .

#### A.3.1 KL-divergence of Markov measure

Consider two Markov path measure $\tilde{M}\ll M$.
Based on the Girsanov’s formula, we can express the Radon-Nikodym derivative as follow:

$$ $\frac{dM}{d\tilde{M}}(\omega)=\frac{dM_{0}}{d\tilde{M}_{0}}(\omega_{0})\exp{ \left(\int_{0}^{\tau}\log\frac{A_{t}^{M}}{A_{t}^{\tilde{M}}}(\omega_{t-}, \omega_{t})dN_{t}(\omega)+\int_{0}^{\tau}(A_{t}^{M}-A_{t}^{\tilde{M}})(\omega_ {t},\omega_{t})dt\right)},$ $$

where $N_{t}(\omega)$ is the number of jumps of the path $\omega$ up to time $t$ and $\omega_{t-}$ is the left limit of the path at time $t$ (see Chazottes et al. () or Appendix 1 of Kipnis & Landim ()).
Note that, due to the compactness of time interval, the number of jumps of each path in $\Omega$ is at most finite, and thus $N_{t}(\omega)$ is bounded.
Also, the escape rate of the state $x$ associated to $M$ is $-A_{t}^{M}(x,x)$.
Thus, we can construct a natural martingale (Lemma 5.1 of Kipnis & Landim ())

$$ $N_{t}+\int_{0}^{t}A_{s}(\omega_{s},\omega_{s})ds,$ $$

which is zero-mean process.

For any continuous and bounded function $\phi:\mathcal{X}\to\mathbb{R}$, we can change the integrator $dN_{t}$ as follow:

$$ $\mathbb{E}_{M}\left[\int_{0}^{t}\phi(\omega_{s})dN_{s}\right]=\mathbb{E}_{M} \left[\int_{0}^{t}-\phi(\omega_{s})A_{s}^{M}(\omega_{s},\omega_{s})ds\right]= \int_{0}^{t}\mathbb{E}_{x\sim M_{s}}\left[-\phi(x)A_{s}^{M}(x,x)\right]ds.$ $$

The KL-divergence is expectation of logarithm of Radon-Nikodym derivative, which leads to:

$$ $\displaystyle D_{\text{KL}}(M\|\tilde{M})$ $\displaystyle=D_{\text{KL}}(M_{0}\|\tilde{M}_{0})+\mathbb{E}_{M}\left[\int_{0} ^{\tau}\log\frac{A_{t}^{M}}{A_{t}^{\tilde{M}}}(\omega_{t-},\omega_{t})dN_{t}( \omega)+\int_{0}^{\tau}(A_{t}^{M}-A_{t}^{\tilde{M}})(\omega_{t},\omega_{t})dt \right],$ $\displaystyle=D_{\text{KL}}(M_{0}\|\tilde{M}_{0})+\int_{0}^{\tau}\mathbb{E}_{x \sim M_{s}}\left[-A_{s}^{M}(x,x)\sum_{y\neq x}p_{s}(x,y)\log\frac{A_{s}^{M}}{A _{s}^{\tilde{M}}}(x,y)\right]ds$ $\displaystyle+\int_{0}^{\tau}\mathbb{E}_{x\sim M_{s}}\left[(A_{t}^{M}-\tilde{A }_{t}^{\tilde{M}})(x,x)\right]ds,$ $$

where $p_{s}(x,y)$ is the probability of jump from $x$ to $y$ given that a jump occurs, which is $\frac{A_{s}^{M}(x,y)}{-A_{s}^{M}(x,x)}$
By applying this, we obtain KL-divergence represented solely in terms of transition rate $A^{M}$ and $A^{\tilde{M}}$:

$$ $D_{\text{KL}}(M\|\tilde{M})=D_{\text{KL}}(M_{0}\|\tilde{M}_{0})+\int_{0}^{\tau }\mathbb{E}_{x\sim M_{s}}\left[\sum_{y\neq x}A_{s}^{M}(x,y)\log\frac{A_{s}^{M} }{A_{s}^{\tilde{M}}}(x,y)+(A_{t}^{M}-A_{t}^{\tilde{M}})(x,x)\right]ds.$ (10) $$

#### A.3.2 KL-divergence of reciprocal measure to Markov measure

###### Lemma A.3 .

If a reciprocal measure $\Lambda\in\mathcal{R}(\mathbb{Q})$ is conditioned on $X_{0}$ being a.s. constant, then the corresponding measure $\Lambda_{\cdot|0}$ is Markov.
For any $M\in\mathcal{M}$ such that $M_{0}=\Lambda_{0}$ and $\Lambda\ll M$, the KL-divergence $D_{\text{KL}}(\Lambda\|M)$ disintegrates as follow:

$$ $D_{\text{KL}}(\Lambda\|M)=\mathbb{E}_{\Lambda_{0}}\left[D_{\text{KL}}(\Lambda_ {\cdot|0}\|M_{\cdot|0})\right].$ $$

###### Proof.

According to proposition 2.5 of Léonard et al. () and lemma 1.4 of Jamison (), conditioning $\Lambda$ on $X_{0}$ not only preserves its reciprocal property, but also transforms it into a Markov process.
Due to the absolute continuity together with $|\mathcal{X}|<\infty$, the KL-divergence is finite.
Then, the KL-divergence can be reformulated as follows:

$$ $\displaystyle D_{\text{KL}}(\Lambda\|M)$ $\displaystyle=\mathbb{E}_{\Lambda}\left[\frac{d\Lambda}{dM}\right],$ $\displaystyle=\mathbb{E}_{\Lambda}\left[\frac{d\Lambda_{\cdot|0}}{dM_{\cdot|0} }\right],$ $\displaystyle=\mathbb{E}_{\Lambda_{0}}\mathbb{E}_{\Lambda_{\cdot|0}}\left[ \frac{d\Lambda_{\cdot|0}}{dM_{\cdot|0}}\right],$ $\displaystyle=\mathbb{E}_{\Lambda_{0}}\left[D_{\text{KL}}(\Lambda_{\cdot|0}\|M _{\cdot|0})\right].$ $$

∎

According to Lemma A.3, while a reciprocal measure $\Lambda\in\mathcal{R}(\mathbb{Q})$ is in general non-Markov, the conditional measure $\Lambda_{\cdot|0}$ is Markov.
Note that we can compute the KL-divergence between two Markov measure based on Equation 10.

###### Lemma A.3 .

###### Proof.

#### A.3.3 Generator of conditioned process

To compute $D_{\text{KL}}(\Lambda\|M)$ based on Lemma A.3, Equation 10, we need the generator of the conditioned process $\Lambda_{\cdot|0}$.
Before deriving the generator of $\Lambda_{\cdot|0}$, we first consider the pinned process of $\mathbb{Q}$ conditioned on $X_{\tau}=z$ in prior.

###### Lemma A.4 .

Let ${(X_{t})_{0\leq t\leq\tau}}$ be a Markov process under the reference measure $\mathbb{Q}$ with transition probability $P_{s:t}(\cdot,\cdot),(s\leq t)$, and generator $A_{s}(\cdot,\cdot)$.
Consider the process conditioned on $X_{\tau}=z$ with corresponding measure denoted by $\mathbb{Q}^{(z)}$.
Then, the conditioned process is also Markov, and its generator is given by:

$$ $A_{s}(x,y;z)=A_{s}(x,y)\frac{P_{s:\tau}(y,z)}{P_{s:\tau}(x,z)}-\delta_{xy} \left[\sum_{u}A_{s}(x,u)\frac{P_{s:\tau}(u,z)}{P_{s:\tau}(x,z)}\right].$ $$

###### Proof.

The conditional probability of $X_{t}$ given the natural filtration $\mathcal{F}_{u}$ and $X_{s}$ with $u\leq s\leq t\leq\tau$ under the measure $\mathbb{Q}^{(z)}$ is as follows:

$$ $\displaystyle\mathbb{Q}^{(z)}(X_{t}=y|X_{s}=x,\mathcal{F}_{u})$ $\displaystyle=\mathbb{Q}(X_{t}=y|X_{s}=x,X_{\tau}=z,\mathcal{F}_{u}),$ $\displaystyle=\mathbb{Q}(X_{t}=y|X_{s}=x,X_{\tau}=z),$ $\displaystyle(\because\mathbb{Q}\in\mathcal{M})$ $\displaystyle=\mathbb{Q}^{(z)}(X_{t}=y|X_{s}=x),$ $$

which confirms $\mathbb{Q}^{(z)}$ is Markov.

Next, the transition probability of $\mathbb{Q}^{(z)}$, denoted by $P_{s:t}(x,y;z)$, is derived as:

$$ $\displaystyle P_{s:t}(x,y;z)$ $\displaystyle=\mathbb{Q}^{(z)}(X_{t}=y|X_{s}=x),$ $\displaystyle=\mathbb{Q}(X_{t}=y|X_{s}=x)\frac{\mathbb{Q}(X_{\tau}=z|X_{t}=y)} {\mathbb{Q}(X_{\tau}=z|X_{s}=x)},$ $\displaystyle=P_{s:t}(x,y)\frac{P_{t:\tau}(y,z)}{P_{s:\tau}(x,z)}.$ $$

Note that, we assumed the measure $\mathbb{Q}$ construct bridge everywhere, the probability ratio has finite value.

Finally, the corresponding generator $A_{s}(x,y;z)$ is obtained using the Kolmogorov forward equation:

$$ $\displaystyle A_{s}(x,y;z)$ $\displaystyle=\frac{\partial}{\partial t}P_{s:t}(x,y;z)\bigg{|}_{t=s},$ $\displaystyle=A_{s}(x,y)\frac{P_{s:\tau}(y,z)}{P_{s:\tau}(x,z)}+\delta_{xy} \left[\frac{\partial_{s}P_{s:\tau}(y,z)}{P_{s:\tau}(x,z)}\right],$ $\displaystyle=A_{s}(x,y)\frac{P_{s:\tau}(y,z)}{P_{s:\tau}(x,z)}-\delta_{xy} \left[\sum_{u}A_{s}(x,u)\frac{P_{s:\tau}(u,z)}{P_{s:\tau}(x,z)}\right].$ $$

∎

We now consider the transition probability and generator of conditioned measure $\Lambda_{\cdot|0}$ of $\Lambda\in\mathcal{R}(\mathbb{Q})$.

###### Lemma A.5 .

For a reciprocal measure $\Lambda\in\mathcal{R}(\mathbb{Q})$, the conditioned process with $X_{0}=x_{0}$ is denoted by $\Lambda_{\cdot|0=x_{0}}$.
The generator of $\Lambda_{\cdot|0=x_{0}}$ is given by the conditional expectation:

$$ $A_{s}^{\Lambda_{|0=x_{0}}}(x,y)=\mathbb{E}_{z\sim\Lambda_{\tau|0,s}}\left[A_{s }(x,y;z)|X_{0}=x_{0},X_{s}=x\right],$ $$

where $A_{s}(x,y;z)$ is the generator of the conditioned process of $\mathbb{Q}$ with $X_{\tau}=z$.

###### Proof.

We denote the transition probability conditioned on $X_{0}=x_{0}$ by $P^{\Lambda_{\cdot|0=x_{0}}}_{s:t}$, which is computed as follows:

$$ $P^{\Lambda_{\cdot|0=x_{0}}}_{s:t}(x,y)=\int_{\mathcal{X}}\Lambda_{\cdot|0=x_{0 }}(X_{\tau}=z|X_{s}=x)\Lambda_{\cdot|0=x_{0}}(X_{t}=y|X_{s}=x,X_{\tau}=z)dz.$ $$

The first term in the integrand is

$$ $\displaystyle\Lambda_{\cdot|0=x_{0}}(X_{\tau}=z|X_{s}=x)$ $\displaystyle=\Lambda(X_{\tau}=z|X_{0}=x_{0},X_{s}=x),$ $\displaystyle=\frac{\Lambda(X_{\tau}=z|X_{0}=x_{0})\Lambda(X_{s}=x|X_{0}=x_{0} ,X_{\tau}=z)}{\Lambda(X_{s}=x|X_{0}=x_{0})},$ $\displaystyle=\frac{\nu_{\tau}(z;x_{0})}{\nu_{s}(x;x_{0})}P_{0:s}(x_{0},x;z),$ $\displaystyle=\frac{\nu_{\tau}(z;x_{0})}{\nu_{s}(x;x_{0})}P_{0:s}(x_{0},x) \frac{P_{s:\tau}(x,z)}{P_{0:\tau}(x_{0},z)},$ $\displaystyle=\frac{\nu_{\tau}(z;x_{0})}{\nu_{s}(x;x_{0})}\frac{\mu_{s}(x;x_{0 })}{\mu_{\tau}(z;x_{0})}P_{s:\tau}(x,z),$ $$

where $\nu_{s}(\cdot;x_{0}),\mu_{s}(\cdot;x_{0})$ are the probability mass functions of $\Lambda_{s|0},\mathbb{Q}_{s|0}$, respectively, with $X_{0}=x_{0}$.
Given the initial-terminal condition, $\Lambda$ is equivalent to the reference $\mathbb{Q}$ based on the definition of reciprocal class Definition 2.1.
Similarly, the second term is same as $\mathbb{Q}(X_{t}=y|X_{0}=x_{0},X_{s}=x,X_{\tau}=z)$, which can be expressed as $P_{s:t}(x,y;z)$ due to the Markov property of $\mathbb{Q}$.
Therefore, the transition probability is:

$$ $\displaystyle P^{\Lambda_{\cdot|0=x_{0}}}_{s:t}(x,y)$ $\displaystyle=\frac{\mu_{s}(x;x_{0})}{\nu_{s}(x;x_{0})}P_{s:t}(x,y)\int_{ \mathcal{X}}P_{t:\tau}(y,z)\frac{\nu_{\tau}(z;x_{0})}{\mu_{\tau}(z;x_{0})}dz.$ $$

Accordingly, the generator is derived as:

$$ $\displaystyle A^{\Lambda_{\cdot|0=x_{0}}}_{s}(x,y)$ $\displaystyle=\partial_{t}P^{\Lambda_{\cdot|0=x_{0}}}_{s:t}(x,y)\bigg{|}_{t=s},$ $\displaystyle=\frac{\mu_{s}(x;x_{0})}{\nu_{s}(x;x_{0})}\partial_{t}P_{s:t}(x,y )\bigg{|}_{t=s}\int_{\mathcal{X}}P_{t:\tau}(y,z)\frac{\nu_{\tau}(z;x_{0})}{\mu _{\tau}(z;x_{0})}dz$ $\displaystyle\quad+\frac{\mu_{s}(x;x_{0})}{\nu_{s}(x;x_{0})}P_{s:t}(x,y)\int_{ \mathcal{X}}\partial_{t}P_{t:\tau}(y,z)\bigg{|}_{t=s}\frac{\nu_{\tau}(z;x_{0}) }{\mu_{\tau}(z;x_{0})}dz,$ $\displaystyle=\frac{\mu_{s}(x;x_{0})}{\nu_{s}(x;x_{0})}A_{s}(x,y)\int_{ \mathcal{X}}P_{s:\tau}(y,z)\frac{\nu_{\tau}(z;x_{0})}{\mu_{\tau}(z;x_{0})}dz$ $\displaystyle\quad-\delta_{xy}\frac{\mu_{s}(x;x_{0})}{\nu_{s}(x;x_{0})}\int_{ \mathcal{X}}\sum_{u}\left[A_{s}(y,u)P_{s:\tau}(u,z)\right]\frac{\nu_{\tau}(z;x _{0})}{\mu_{\tau}(z;x_{0})}dz,$ $\displaystyle=\int_{\mathcal{X}}\frac{\mu_{s}(x;x_{0})}{\nu_{s}(x;x_{0})}\left [A_{s}(x,y)P_{s:\tau}(y,z)-\delta_{xy}\sum_{u}A_{s}(y,u)P_{s:\tau}(u,z)\right] \frac{\nu_{\tau}(z;x_{0})}{\mu_{\tau}(z;x_{0})}dz,$ $\displaystyle=\int_{\mathcal{X}}\frac{\nu_{\tau}(z;x_{0})}{\mu_{\tau}(z;x_{0}) }\frac{\mu_{s}(x;x_{0})}{\nu_{s}(x;x_{0})}P_{s:\tau}(x,z)A_{s}(x,y;z)dz,$ $\displaystyle=\int_{\mathcal{X}}\Lambda_{|0}(X_{\tau}=z|X_{s}=x)A_{s}(x,y;z)dz,$ $\displaystyle=\mathbb{E}_{z\sim\Lambda_{\tau|0,s}}[A_{s}(x,y;z)|X_{0}=x_{0},X_ {s}=x].$ $$

In conclusion, the generator of $\Lambda_{\cdot|0}$ is given as the conditional expectation,

$$ $A_{s}^{\Lambda_{\cdot|0}}(x,y)=\mathbb{E}_{\Lambda_{\tau|0,s}}[A_{s}(x,y;X_{ \tau})|X_{0},X_{s}=x].$ $$

∎

###### Lemma A.4 .

###### Proof.

###### Lemma A.5 .

###### Proof.

#### A.3.4 Marginal distribution of mixture of Markov chains

A reciprocal process $\Lambda\in\mathcal{R}(\mathbb{Q})$ can be represented as a mixture of Markov chains as:

$$ $\Lambda(\cdot)=\sum_{x,y}\mathbb{Q}(\cdot|X_{0}=x,X_{\tau}=y)\Lambda_{0,\tau}( x,y).$ $$

We here propose the mixture representation of Markov chain similar to theorem 2 of Peluchetti () which describing mixture of diffusion process over Euclidean space.

###### Proposition A.6 .

Consider a family of Markov measures on $\Omega$ indexed by $u\in\mathcal{I}$

$$ $\displaystyle\partial_{t}P_{s:t}^{u}(x,y)$ $\displaystyle=\sum_{z}P_{s:t}^{u}(x,z)A_{t}^{u}(z,y),$ $\displaystyle A_{s}^{u}(x,y)$ $\displaystyle=\partial_{t}P_{s,t}^{u}(x,y)\bigg{|}_{t=s},$ $$

where $P^{u}$ and $A^{u}$ is the transition probability and generator for each measure associated to the process $X^{u}=(X_{t}^{u})_{t\in[0,\tau]}$.
We here assume that $A_{s}^{u}$ is finite for every $u$.
Let the corresponding Markov measure be $M^{u}$, and $\mu^{u}_{t}$ be the density of $M^{u}_{t}$, the marginal distribution of $X_{t}^{u}$.
Let a mixture of $M^{u}$ with the index distribution $\mathcal{U}$ over $\mathcal{I}$ be $\Lambda$,

$$ $\Lambda(\cdot)=\int_{\mathcal{I}}M^{u}(\cdot)\mathcal{U}(du).$ $$

We denote mixture marginal density $\mu_{t}$ and mixture initial distribution $\Lambda_{0}$ as:

$$ $\displaystyle\mu_{t}(x)=\int_{\mathcal{I}}\mu_{t}^{u}(x)\mathcal{U}(du),$ $\displaystyle\Lambda_{0}(\cdot)=\int_{\mathcal{I}}M^{u}_{0}(\cdot)\mathcal{U}( du).$ $$

Let $X=(X_{t})_{t\in[0,\tau]}$ be another Markov chain generated by:

$$ $\displaystyle\partial_{t}P_{s:t}(x,y)=\sum_{z}P_{s:t}(x,z)A_{t}(z,y),$ $\displaystyle A_{t}(x,y)=\frac{1}{\mu_{t}(x)}\int_{\mathcal{I}}A_{t}^{u}(x,y) \mu_{t}^{u}(x)\mathcal{U}(du),$ $\displaystyle X_{0}\sim\Lambda_{0}.$ $$

Then, the marginal density of $X_{t}$ is $\mu_{t}$.
It is assumed that exchange of $\partial_{t}$ and $\int_{\mathcal{I}}\mathcal{U}(du)$ is justified.

###### Proof.

We start from verifying $A_{t}(x,y)$ admits conditions of transition rate function.
For $x\neq y$ it is trivial that $A_{t}(x,y)$ is finite and non-negative.
Also, $\sum_{y\in\mathcal{X}}A_{t}(x,y)$ becomes zeros for all $x$ and $t$ as:

$$ $\displaystyle\sum_{y\in\mathcal{X}}A_{t}(x,y)$ $\displaystyle=\sum_{y\in\mathcal{X}}\frac{1}{\mu_{t}(y)}\int_{\mathcal{I}}A_{t }^{u}(x,y)\mu_{t}^{u}(x)\mathcal{U}(du),$ $\displaystyle=\frac{1}{\mu_{t}(y)}\int_{\mathcal{I}}\sum_{y\in\mathcal{X}}A_{t }^{u}(x,y)\mu_{t}^{u}(x)\mathcal{U}(du),$ $\displaystyle=0.$ $$

Thus, mixture of generators $A_{t}(x,y)$ holds the condition for a generator of Markov measures.

Next, for $t\in(0,\tau)$,

$$ $\displaystyle\partial_{t}\mu_{t}(x)$ $\displaystyle=\partial_{t}\int_{\mathcal{I}}\mu_{t}^{u}(x)\mathcal{U}(du),$ $\displaystyle=\int_{\mathcal{I}}\partial_{t}\mu_{t}^{u}(x)\mathcal{U}(du),$ $\displaystyle(\text{assumption)}$ $\displaystyle=\int_{\mathcal{I}}\sum_{y\in\mathcal{X}}A^{u}_{t}(y,x)\mu_{t}^{u }(y)\mathcal{U}(du),$ $\displaystyle(\because\text{Kolmogorov equation})$ $\displaystyle=\sum_{y\in\mathcal{X}}\left(\int_{\mathcal{I}}A^{u}_{t}(y,x)\mu_ {t}^{u}(y)\mathcal{U}(du)\right),$ $\displaystyle=\sum_{y\in\mathcal{X}}A_{t}(y,x)\mu_{t}(y).$ $$

The equality of the left-hand side and the final line corresponds to the Kolmogorov equation for the process $X$ generated by $A_{t}$.
This shows that $\mu_{t}$ is the marginal distribution of $X_{t}$.
∎

###### Proposition A.6 .

###### Proof.

#### A.3.5 Minimizer of the KL-divergence

The Markov projection of a reciprocal process $\Lambda$, denoted $\Pi_{\mathcal{M}}(\Lambda)$, is a Markov process $M^{*}$ which minimizes the reverse KL-divergence $D_{\text{KL}}(\Lambda,M)$.
We here characterize $M^{*}$ by specifying its generator.
While the generator $A_{t}^{\Lambda_{\cdot|0}}$ of $\Lambda_{\cdot|0}$ represented as a conditional expectation given $X_{0},X_{t}$ as noted in Lemma A.5, that of the $M^{*}$ is supposed to be represented as the conditional expectation without $X_{0}$ condition.
The following lemma says that $\mathbb{E}_{\Lambda_{0|t}}\left[A_{t}^{\Lambda_{\cdot|0}}(X_{t},y)\right]$ is the generator.

###### Lemma A.7 .

Let the Markov projection of a reciprocal process $\Lambda\in\mathcal{R}(\mathbb{Q})$ be $M^{*}=\Pi_{\mathcal{M}}(\Lambda)$.
Then the generator of $M^{*}$ is $A_{t}^{M^{*}}(X_{t},y)=\mathbb{E}_{\Lambda_{\tau|t}}\left[A_{t}(X_{t},y;X_{
\tau})\bigg{|}X_{t}\right]$.

###### Proof.

Because $\mathbb{Q}$ is assumed to be able to construct bridge everywhere, $A_{t}(x,y;z)=0\iff A_{t}(x,y)=0\quad\forall x\neq y,z\in\mathcal{X},t\in[0,\tau)$.

We claim that the Markov measure $M$ generated by $A_{t}^{M}(x,y)=\mathbb{E}_{\Lambda_{\tau|t}}\big{[}A_{t}(X_{t},y;X_{\tau})\big
{|}X_{t}=x\big{]}$ with $M_{0}=\Lambda_{0}$ is a.s. $M^{*}$.
We can re-formulate the generator as:

$$ $\displaystyle A_{t}^{M}(X_{t},y)$ $\displaystyle=\mathbb{E}_{\Lambda_{\tau|t}}\left[A_{t}(X_{t},y;X_{\tau})\big{| }X_{t}\right],$ $\displaystyle=\sum_{x_{\tau}}\Lambda(X_{\tau}=x_{\tau}|X_{t})A_{t}(X_{t},y;x_{ \tau}),$ $\displaystyle=\sum_{x_{\tau}}\sum_{x_{0}}\Lambda(X_{0}=x_{0},X_{\tau}=x_{\tau} |X_{t})A_{t}(X_{t},y;x_{\tau}),$ $\displaystyle=\sum_{x_{0}}\Lambda(X_{0}=x_{0}|X_{t})\sum_{x_{\tau}}\Lambda(X_{ \tau}=x_{\tau}|X_{0}=x_{0},X_{t})A_{t}(X_{t},y;x_{\tau}),$ $\displaystyle=\sum_{x_{0}}\Lambda(X_{0}=x_{0}|X_{t})A_{t}^{\Lambda_{\cdot|0}}( X_{t},y),$ $\displaystyle=\mathbb{E}_{\Lambda_{0|t}}\left[A_{t}^{\Lambda_{\cdot|0}}(X_{t}, y)\big{|}X_{t}\right],$ (11) $$

which conclude that it becomes conditional expectation of the generator of $\Lambda_{\cdot|0}$.
Also, based on the Proposition A.6, the $M_{t}=\Lambda_{t}$ for all $t$.

Note that $A_{t}^{\Lambda_{\cdot|0}}(x,y)=\mathbb{E}_{\Lambda_{\tau|0,t}}\big{[}A_{t}(x,y
;X_{\tau})\big{|}X_{0},X_{t}=x\big{]}$, which implies $A_{t}^{M}(x,y)=0\Longrightarrow A_{t}^{\Lambda_{\cdot|0}}(x,y)=0\quad\forall x
\neq y,z\in\mathcal{X},t\in[0,\tau)$.
Thus, $\Lambda\ll M$.
Because $M^{*}$ is the minimizer of KL-divergence $D_{\text{KL}}(\Lambda\|\cdot)$, we can assume $M^{*}_{0}=\Lambda_{0}$.
Based on this, we want to show that $D_{\text{KL}}(\Lambda\|M^{*})-D_{\text{KL}}(\Lambda\|M)=0$.
Recall that Equation 10, where the KL-divergence is formulated as

$$ $\displaystyle D_{\text{KL}}(\Lambda\|M)$ $\displaystyle=\int_{0}^{\tau}\mathbb{E}_{\Lambda_{0,t}}\left[\left(A_{t}^{ \Lambda_{\cdot|0}}-A_{t}^{M}\right)(X_{t},X_{t})+\sum_{y\neq X_{t}}A_{t}^{ \Lambda_{\cdot|0}}\log\frac{A_{t}^{\Lambda_{\cdot|0}}}{A_{t}^{M}}(X_{t},y) \right]dt,$ $\displaystyle=\int_{0}^{\tau}\underbrace{\mathbb{E}_{\Lambda_{0,t}}\left[\sum_ {y\neq X_{t}}\left(A_{t}^{\Lambda_{\cdot|0}}\log\frac{A_{t}^{\Lambda_{\cdot|0} }}{A_{t}^{M}}-A_{t}^{\Lambda_{\cdot|0}}+A_{t}^{M}\right)(X_{t},y)\right]}_{f(t ,\Lambda,M)}dt.$ $$

Then, $D_{\text{KL}}(\Lambda\|M^{*})-D_{\text{KL}}(\Lambda\|M)=\int_{0}^{\tau}\Delta
dt\leq
0$ by definition, where $\Delta:=f(t,\Lambda,M^{*})-f(t,\Lambda,M)$.

$$ $\displaystyle\Delta$ $\displaystyle=\mathbb{E}_{\Lambda_{0,t}}\left[\sum_{y\neq X_{t}}\left(A_{t}^{M ^{*}}-A_{t}^{M}+A_{t}^{\Lambda_{\cdot|0}}\log\frac{A_{t}^{M}}{A_{t}^{M^{*}}} \right)(X_{t},y)\right],$ $\displaystyle=\mathbb{E}_{\Lambda_{t}}\mathbb{E}_{\Lambda_{0|t}}\left[\sum_{y \neq X_{t}}\left(A_{t}^{M^{*}}-A_{t}^{M}+A_{t}^{\Lambda_{\cdot|0}}\log\frac{A_ {t}^{M}}{A_{t}^{M^{*}}}\right)(X_{t},y)\right],$ $\displaystyle=\mathbb{E}_{\Lambda_{t}}\sum_{y\neq X_{t}}\mathbb{E}_{\Lambda_{0 |t}}\left[\left(A_{t}^{M^{*}}-A_{t}^{M}+A_{t}^{\Lambda_{\cdot|0}}\log\frac{A_{ t}^{M}}{A_{t}^{M^{*}}}\right)(X_{t},y)\right],$ $\displaystyle=\mathbb{E}_{\Lambda_{t}}\sum_{y\neq X_{t}}\left(A_{t}^{M^{*}}-A_ {t}^{M}+\mathbb{E}_{\Lambda_{0|t}}\left[A_{t}^{\Lambda_{\cdot|0}}|X_{t}\right] \log\frac{A_{t}^{M}}{A_{t}^{M^{*}}}\right)(X_{t},y),$ $\displaystyle(\because M,M^{*}\in\mathcal{M})$ $\displaystyle=\mathbb{E}_{\Lambda_{t}}\sum_{y\neq X_{t}}\left(A_{t}^{M^{*}}-A_ {t}^{M}+A_{t}^{M}\log\frac{A_{t}^{M}}{A_{t}^{M^{*}}}\right)(X_{t},y),$ $\displaystyle(\because\lx@cref{creftype~refnum}{eq:tmp})$ $\displaystyle=\mathbb{E}_{M_{t}}\sum_{y\neq X_{t}}\left(A_{t}^{M^{*}}-A_{t}^{M }+A_{t}^{M}\log\frac{A_{t}^{M}}{A_{t}^{M^{*}}}\right)(X_{t},y).$ $\displaystyle(\because M_{t}=\Lambda_{t})$ $$

We can deduce that $\int_{0}^{\tau}\Delta dt=D_{\text{KL}}(M\|M^{*})\geq 0$, which conclude that $M=M^{*}$.

∎

###### Lemma A.7 .

###### Proof.

#### A.3.6 Proof of Proposition A.2

Now we prove the proposition Proposition A.2 using above lemmas.

###### Proof.

Proposition A.2

By the Lemma A.3 and Equation 10, the KL-divergence of $\Lambda$ to any Markov measure $M\in\mathcal{M}$ disintegrates as:

$$ $D_{\text{KL}}(\Lambda\|M)=\int_{0}^{\tau}\mathbb{E}_{\Lambda_{0,t}}\left[\left (A_{t}^{\Lambda_{\cdot|0}}-A_{t}^{M}\right)(X_{t},X_{t})+\sum_{y\neq X_{t}}A_{ t}^{\Lambda_{\cdot|0}}\log\frac{A_{t}^{\Lambda_{\cdot|0}}}{A_{t}^{M}}(X_{t},y) \right]dt,$ $$

where the $A_{t}^{\Lambda_{\cdot|0}}$ is the generator of pinned process of $\Lambda_{\cdot|0}$ stated by Lemma A.5.
According to Lemma A.7, the Markov measure $M^{*}$ generated by $A_{t}(x,y)=\mathbb{E}_{\Lambda_{\tau|t}}\left[A_{t}(x,y;X_{\tau})|X_{t}=x\right]$
is the minimizer of $D_{\text{KL}}(\Lambda\|M)$ for $M\in\mathcal{M}$.
In last, Proposition A.6 ensure the time marginals of $\Lambda,M^{*}$ for all $t\in[0,\tau]$ are equivalent.
∎

###### Proof.

### A.4 Time-reversal Markov projection

###### Proposition A.8 .

Let $M^{*}=\Pi_{\mathcal{M}}(\Lambda),\Lambda\in\mathcal{R}(\mathbb{Q})$ and the forward generator of $\mathbb{Q}$ be $A_{t}(x,y)$ with the transition probability $P_{s:t}(x,y)$.
Let $X=(X_{t})_{t\in[0,\tau]}$ be the canonical process and $Y=(Y_{t})_{t\in[0,\tau]}$ be the time reversal of $X$, where $Y_{t}=X_{\tau-t}$.
Under mild assumptions, the time-reversal generator of $M^{*}$ becomes

$$ $\displaystyle\tilde{A}_{t}^{M^{*}}(Y_{t},x)$ $\displaystyle=\mathbb{E}_{0|t}\left[\tilde{A}_{t}(Y_{t},x;X_{0})\bigg{|}Y_{t}\right]$ $\displaystyle\tilde{A}_{s}(y,x;z)$ $\displaystyle=\partial_{t}\tilde{P}_{s:t}(y,x;z)\bigg{|}_{t=s}$ $\displaystyle=A_{\tau-s}(x,y)\frac{P_{0:\tau-s}(z,x)}{P_{0:\tau-s}(z,y)}- \delta_{xy}\sum_{u}A_{\tau-s}(u,x)\frac{P_{0:\tau-s}(z,u)}{P_{0:\tau-s}(z,x)},$ $$

where $z\in\mathcal{X}$, and $\tilde{A}_{s}(\cdot,\cdot;z)$ is the generator for the conditioned process of $Y$ with $Y_{\tau}=z$.
The reverse KL-divergence is

$$ $D_{\text{KL}}(\Lambda\|M^{*})=\int_{0}^{\tau}\mathbb{E}_{\Lambda_{\tau,t}} \left[(\tilde{A}^{\Lambda_{\cdot|\tau}}-\tilde{A}^{M^{*}})(Y_{t},Y_{t})+\sum_{ x\neq Y_{t}}\tilde{A}^{\Lambda_{\cdot|\tau}}\log{\frac{\tilde{A}^{\Lambda_{ \cdot|\tau}}}{\tilde{A}^{M^{*}}}}(Y_{t},x)\right]dt,$ $$

where the $\tilde{A}^{\Lambda_{\cdot|\tau}}$ is the time-reversal generator for the conditioned measure $\Lambda_{\cdot|\tau}$ which is defined as

$$ $\tilde{A}_{t}^{\Lambda_{\cdot|\tau}}(Y_{t},x)=\mathbb{E}_{\Lambda_{0|t,\tau}} \left[\tilde{A}(Y_{t},x;Y_{\tau})\bigg{|}Y_{t},Y_{0}\right].$ $$

###### Proof.

Follow the proof of Proposition A.2. Note that the Markov measure is time-symmetric.
∎

###### Proposition A.8 .

###### Proof.

### A.5 Proof of Theorem 2.3

###### Proof.

With iterative projection, we have a sequence of measure $\Lambda^{(n)}$ such that:

$$ $\displaystyle\Lambda^{(2n+1)}=\Pi_{\mathcal{M}}(\Lambda^{(2n)}),$ $\displaystyle\Lambda^{(2n+2)}=\Pi_{\mathcal{R}(\mathbb{Q})}(\Lambda^{(2n+1)}),$ $$

where $\Lambda^{(0)}\in\mathcal{R}(\mathbb{Q})$.
Let $\Pi^{(n)}=\Lambda^{(2n)}$ and $M^{(n)}=\Lambda^{(2n+1)}$.
We will omit the superscription $\cdot^{(n)}$ if there is no confusion.

Let $\mathbb{P}\in\mathcal{M}$ be a Markov process generated by a generator $A^{\mathbb{P}}_{t}(x,y)<\infty$ with initial distribution $M_{0}$ with the Kolmogorov forward equation.
Assuming that $D_{\text{KL}}(\Pi\|\mathbb{P})<\infty$, then

$$ $D_{\text{KL}}(\Pi\|\mathbb{P})-D_{\text{KL}}(\Pi\|M)=\mathbb{E}_{\Pi}\left[ \log\frac{dM}{d\mathbb{P}}\right]<\infty,$ $$

which implying the Radon-Nikodym derivative $\frac{dM}{d\mathbb{P}}<\infty$ over the support of $\Pi$.

Then we will first show the equality:

$$ $\mathbb{E}_{\Pi}\left[\log\frac{dM}{d\mathbb{P}}\right]=\mathbb{E}_{M}\left[ \log\frac{dM}{d\mathbb{P}}\right].$ $$

As $\Pi_{0}=M_{0}$, it is sufficient to show

$$ $\mathbb{E}_{\Pi}\left[\log\frac{dM_{\cdot|0}}{d\mathbb{P}_{\cdot|0}}\right]= \mathbb{E}_{M}\left[\log\frac{dM_{\cdot|0}}{d\mathbb{P}_{\cdot|0}}\right],$ $$

where

$$ $\log\frac{dM_{\cdot|0}}{d\mathbb{P}_{\cdot|0}}(\omega)=\int_{0}^{\tau}\log \frac{A_{t}^{M}}{A_{t}^{\mathbb{P}}}(\omega_{t-},\omega_{t})dN_{t}(\omega)+ \int_{0}^{\tau}(A_{t}^{M}-A_{t}^{\mathbb{P}})(\omega_{t},\omega_{t})dt$ $$

is given from Section A.3.1.
Since $\Pi_{|0}$ is Markov from Lemma A.3 with generator

$$ $A_{t}^{\Pi_{|0}}(x,y)=\mathbb{E}_{\Pi_{\tau|0,t}}\left[A_{t}(x,y;X_{\tau})|X_{ 0},X_{t}=x\right],$ $$

which is derived at the Lemma A.5.
Then,

$$ $\displaystyle\mathbb{E}_{\Pi}\left[\int_{0}^{\tau}\log\frac{A_{t}^{M}}{A_{t}^{ \mathbb{P}}}(X_{t-},X_{t})dN_{t}\right]$ $\displaystyle=\mathbb{E}_{\Pi}\left[\int_{0}^{\tau}\sum_{y\neq X_{t}}A_{t}^{ \Pi_{\cdot|0}}\log\frac{A_{t}^{M}}{A_{t}^{\mathbb{P}}}(X_{t},y)dt\right],$ $\displaystyle=\mathbb{E}_{\Pi}\left[\int_{0}^{\tau}\sum_{y\neq X_{t}}\mathbb{E }_{\Pi}\left[A_{t}(X_{t},y;X_{\tau})|X_{0},X_{t}\right]\log\frac{A_{t}^{M}}{A_ {t}^{\mathbb{P}}}(X_{t},y)dt\right],$ $\displaystyle=\mathbb{E}_{\Pi}\left[\int_{0}^{\tau}\sum_{y\neq X_{t}}\mathbb{E }_{\Pi}\left[A_{t}(X_{t},y;X_{\tau})|X_{t}\right]\log\frac{A_{t}^{M}}{A_{t}^{ \mathbb{P}}}(X_{t},y)dt\right],$ $$

by the tower property of conditional expectation.
From the Lemma A.7 and Proposition A.6, we can replace the $\mathbb{E}_{\Pi}\left[A_{t}(X_{t},y;X_{\tau})|X_{t}\right]$ as $A_{t}^{M}(X_{t},y)$,

$$ $\displaystyle\mathbb{E}_{\Pi}\left[\int_{0}^{\tau}\sum_{y\neq X_{t}}\mathbb{E} _{\Pi}\left[A_{t}(X_{t},y;X_{\tau})|X_{t}\right]\log\frac{A_{t}^{M}}{A_{t}^{ \mathbb{P}}}(X_{t},y)dt\right],$ $\displaystyle=\mathbb{E}_{\Pi}\left[\int_{0}^{\tau}\sum_{y\neq X_{t}}A_{t}^{M} \log\frac{A_{t}^{M}}{A_{t}^{\mathbb{P}}}(X_{t},y)dt\right],$ $\displaystyle=\mathbb{E}_{M}\left[\int_{0}^{\tau}\sum_{y\neq X_{t}}A_{t}^{M} \log\frac{A_{t}^{M}}{A_{t}^{\mathbb{P}}}(X_{t},y)dt\right],$ $$

where the last equation justified since $\Pi_{t}=M_{t}$ for all $t\in[0,\tau]$.
The other term, $\int_{0}^{\tau}(A_{t}^{M}-A_{t}^{\mathbb{P}})(X_{t},X_{t})dt$, can be treated as the same way, which establishes the equality.
From the result, we get $D_{\text{KL}}(\Pi\|\mathbb{P})=D_{\text{KL}}(\Pi\|M)+D_{\text{KL}}(M\|\mathbb{
P})$.
The equality is derived for other diffusion processes Peluchetti (); Shi et al. (); Liu & Wu (), but in this paper we extended it to continuous Markov chain with discrete state space case.

By letting $\mathbb{P}=\mathbb{P}^{\text{SB}}$, we get

$$ $D_{\text{KL}}(\Pi\|\mathbb{P}^{\text{SB}})\geq D_{\text{KL}}(M\|\mathbb{P}^{ \text{SB}}),$ $$

the equality holds if and only if $\Pi=M$.
By the assumption $D_{\text{KL}}(\Lambda^{(0)}\|\mathbb{P}^{\text{SB}})<\infty$, the sequence of KL-divergence $\{D_{\text{KL}}(\Lambda^{(n)}\|\mathbb{P}^{\text{SB}})\}_{n\in\mathbb{N}}$ is non-increasing and bounded, which implies the KL-divergence converges.
Thus,

$$ $\lim_{n\to\infty}D_{\text{KL}}(\Pi^{(n)}\|\mathbb{P}^{\text{SB}})-D_{\text{KL} }(M^{(n)}\|\mathbb{P}^{\text{SB}})=\lim_{n\to\infty}D_{\text{KL}}(\Pi^{(n)}\|M ^{(n)})=0.$ $$

We here utilize Aldous’ tightness criteria for showing the tightness of $\mathbb{P}^{\text{SB}}$.
Since the state space is finite it is sufficient to show the following:
For any $\varepsilon>0$, there exists $\delta>0$ such that for all stopping time $t$,

$$ $\mathbb{P}^{\text{SB}}(d_{\mathcal{X}}(X_{t+\theta},X_{t})\geq\eta)\leq\varepsilon,$ $$

for arbitrary small $\eta>0$ and $0<\theta<\delta$ (see Section 16 of Billingsley () Section 3 of Ethier & Kurtz ()).
To check the criteria, we will show that the leaving rate of $\mathbb{P}^{\text{SB}}$ is bounded.
From the Lemma A.5 and Lemma A.7, we know the analytic form of the transition rate $A_{t}$.
Since Theorem A.1 ensure that the $\mathbb{P}^{\text{SB}}$ is Markov as well as is in reciprocal class,

$$ $\displaystyle A_{t}^{\mathbb{P}^{\text{SB}}}(x,x)$ $\displaystyle=\mathbb{E}_{\mathbb{P}^{\text{SB}}}\left[A_{t}(x,x;X_{\tau})|X_{ 0}=x_{0},X_{t}=x\right],$ $\displaystyle=\sum_{X_{\tau}\in\mathcal{X}}\bigg{[}\mathbb{P}^{\text{SB}}(X_{ \tau}|X_{t}=x,X_{0}=x_{0})A_{t}(x,x)$ $\displaystyle\quad-\frac{\mathbb{P}^{\text{SB}}(X_{\tau}|X_{t}=x,X_{0})}{P_{t: \tau}(x,X_{\tau})}\sum_{u\in\mathcal{X}}A(x,u)P_{t:\tau}(u,X_{\tau})\bigg{]},$ $\displaystyle=\sum_{X_{\tau}\in\mathcal{X}}\bigg{[}\mathbb{P}^{\text{SB}}(X_{ \tau}|X_{t}=x)A_{t}(x,x)$ $\displaystyle\quad-\frac{\mathbb{P}^{\text{SB}}(X_{\tau}|X_{t}=x,X_{0})}{P_{t: \tau}(x,X_{\tau})}\sum_{u\in\mathcal{X}}A(x,u)P_{t:\tau}(u,X_{\tau})\bigg{]},$ $\displaystyle=A_{t}(x,x)-\sum_{X_{\tau}\in\mathcal{X}}\bigg{[}\frac{\mathbb{P} ^{\text{SB}}(X_{\tau}|X_{t}=x,X_{0})}{P_{t:\tau}(x,X_{\tau})}\sum_{u\in \mathcal{X}}A(x,u)P_{t:\tau}(u,X_{\tau})\bigg{]},$ $\displaystyle=A_{t}(x,x)-\sum_{X_{\tau}\in\mathcal{X}}\bigg{[}\frac{d\mathbb{P }^{\text{SB}}_{t,\tau}}{d\mathbb{Q}_{t,\tau}}(x,X_{\tau})\sum_{u\in\mathcal{X} }A(x,u)P_{t:\tau}(u,X_{\tau})\bigg{]}.$ $$

The Radon-Nikodym derivative is finite and positive for $x$, such that $\mathbb{P}^{\text{SB}}(X_{t}=x)>0$.
Because the state space $\mathcal{X}$ is finite, there is finite upper bound $u_{t}$ for every possible pair $(x,x_{\tau})$.
Thus,

$$ $\displaystyle A_{t}^{\mathbb{P}^{\text{SB}}}(x,x)$ $\displaystyle\geq A_{t}(x,x)-\sum_{X_{\tau}\in\mathcal{X}}\bigg{[}u_{t}\sum_{u \in\mathcal{X}}A(x,u)P_{t:\tau}(u,X_{\tau})\bigg{]},$ $\displaystyle=A_{t}(x,x)-u_{t}\sum_{u\in\mathcal{X}}A(x,u)\sum_{X_{\tau}\in \mathcal{X}}P_{t:\tau}(u,X_{\tau}),$ $\displaystyle=A_{t}(x,x)-u_{t}\sum_{u\in\mathcal{X}}A(x,u),$ $\displaystyle=A_{t}(x,x).$ $$

The leaving rate $c_{t}^{\text{SB}}(x)$ of $\mathbb{P}^{\text{SB}}$ at $(t,x)$ is bounded by that of $Q$, $c_{t}(x)$.
Assuming the leaving rate of $Q$ is uniformly bounded by $c$, $c_{t}^{\text{SB}}$ is uniformly bounded by $c$.
Back to the Aldous’ tightness criteria, by choosing $\delta=\varepsilon/c$ we can see that $\mathbb{P}^{\text{SB}}$ is tight.

The sequence $\{\Pi^{(n)}\}_{n\in\mathbb{N}}$ is tight.
Since $\mathbb{P}^{\text{SB}}$ is tight, for any $\varepsilon>0$ we can choose a compact and measurable $K$ (under Skorokhod topology and associated Borel $\sigma$ algebra).
For any measurable $K$,

$$ $\displaystyle D_{\text{KL}}(\Pi^{(n)}\|\mathbb{P}^{\text{SB}})$ $\displaystyle=\mathbb{E}_{\Pi^{(n)}}\left[-\log\frac{d\mathbb{P}^{\text{SB}}}{ d\Pi^{(n)}}|K^{c}\right]\Pi^{(n)}(K^{c})+\mathbb{E}_{\Pi^{(n)}}\left[-\log \frac{d\mathbb{P}^{\text{SB}}}{d\Pi^{(n)}}|K\right]\Pi^{(n)}(K),$ $\displaystyle\geq-\log\frac{\mathbb{P}^{\text{SB}}(K^{c})}{\Pi^{(n)}(K^{c})} \Pi^{(n)}(K^{c})-\log\frac{\mathbb{P}^{\text{SB}}(K)}{\Pi^{(n)}(K)}\Pi^{(n)}(K),$ $$

by the Jensen inequality.
If $\{\Pi^{n}\}_{n\in\mathbb{N}}$ is not tight, for each compact $K$ and $\lambda>0$, there is at least one $n\in\mathbb{N}$ where $\Pi^{(n)}(K^{c})\geq\lambda,\Pi^{(n)}(K)<1-\lambda$.
Thus, for $\varepsilon>0$ there exists $n\in\mathbb{N}$, such that

$$ $-\log\frac{\mathbb{P}^{\text{SB}}(K^{c})}{\Pi^{(n)}(K^{c})}\Pi^{(n)}(K^{c}) \geq-\log(\varepsilon/\lambda)\lambda,$ $$

which implies that the lower bound of $D_{\text{KL}}(\Pi^{(n)}\|\mathbb{P}^{\text{SB}})$ can be arbitrary large.
However, the KL-divergence is non-increasing and upper bounded by $D_{\text{KL}}(\Pi^{(0)}\|\mathbb{P}^{\text{SB}})<\infty$, contradiction occurs.
Thus, the sequence of measure should be tight.
Similarly, we can show the $\{M^{(n)}\}_{n\in\mathbb{N}}$ is tight.

The space $(\Omega,d_{\Omega})$ is a Polish space (Section 12 of Billingsley () or Section 3.5 of Ethier & Kurtz ()), and therefore, by Prokhorov’s theorem, the collections of measures ${\Pi^{(n)}}_{n\in\mathbb{N}}$ and ${M^{(n)}}_{n\in\mathbb{N}}$ are relatively compact.
Each subsequence of $\{M^{(n)}\}_{n\in\mathbb{N}}$ has a sub-subsequence $\{\Pi^{(i)}\}_{i\geq 0}$ weakly converges to $\Pi^{(\infty)}$ as $i\to\infty$.
Similarly, there is a sub-subsequence $\{M^{(i)}\}_{i\geq 0}$ that converges in law to $M^{(\infty)}$.
By the lower semi-continuity of the KL-divergence, $D_{\text{KL}}(\Pi^{(\infty)}\|M^{(\infty)})\leq\liminf_{i\to\infty}(\Pi^{(i)}
\|M^{(i)})=0$, which implies the two convergence point is equal to $\mathbb{P}^{(\infty)}$.
The resulting measure is Markov and is in $\mathcal{R}(\mathbb{Q})$, because the state space is finite, which deduce $\mathbb{P}^{\text{SB}}=\mathbb{P}^{(\infty)}$.
We choose an arbitrary convergent point of sub-subsequence resulting in $\mathbb{P}^{\text{SB}}$, the convergence is ensured (Theorem 2.6 Billingsley ()).
∎

###### Proof.

## Appendix B Graph permutation matching

### B.1 Introduction

A graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$ is defined as the set of nodes $\mathcal{V}=\{v_{i}\}$ and the set of edges $\{e_{ij}\}$, where $i$ and $j$ denote the node indices.
Under a permutation $\sigma\in S_{n}$, the structure of the graph $\mathcal{G}$ remains unchanged, as the node and edge sets are invariant to indexing:

$$ $\displaystyle\sigma(\mathcal{V})=\{v_{\sigma(i)}\}=\mathcal{V},$ $\displaystyle\sigma(\mathcal{E})=\{e_{\sigma(i)\sigma(j)}\}=\mathcal{E},$ $$

However, the vectorized representation of a graph, $\mathbf{G}=(\mathbf{V},\mathbf{E})$ is affected by the permutation because $\mathbf{V}$ and $\mathbf{E}$ are treated as ordered sets.
Thus, a graph $\mathcal{G}$ can be considered as a set of all permuted version of its vectorized representation:

$$ $\mathcal{G}=\{\sigma(\mathbf{G}):\sigma\in S_{n}\},$ $$

where the $\mathbf{G}$ is an arbitrarily indexed vectorization of $\mathcal{G}$.

The likelihood of transitioning from an initial graph $\mathbf{G}=(\mathbf{V},\mathbf{E})$ to a terminal graph $\mathbf{G}^{\prime}=(\mathbf{V}^{\prime},\mathbf{E}^{\prime})$, denoted $p(\mathbf{G}^{\prime}|\mathbf{G})=\mathbb{Q}(X_{\tau}=\mathbf{G}^{\prime}|X_{0
}=\mathbf{G})$, depends on both node and edge correspondences, whereas the likelihood $p(\mathcal{G}^{\prime}|\mathbf{G})$, which is permutation-invariant, does not.
More specifically, for any permutation $\sigma$, the likelihood $p(\sigma(\mathbf{G}^{\prime})|\mathbf{G})$ changes, even though the underlying graph $\mathcal{G}^{\prime}$ remains unchanged.
Calculating the likelihood of the graph itself (as opposed to the graph vector) would require summing the likelihoods of all possible permutations of $\mathbf{G}^{\prime}$, but this approach is computationally prohibitive due to the factorial number of permutations.

This dependency on graph permutation introduces challenges not only in computing likelihoods but also in constructing reciprocal measures.
Ideally, a reciprocal measure over graph domain would account for all permutations of graph vectors.
However, as with likelihood computations, the construction of reciprocal measures is highly sensitive on the alignment of graph vectors, making proper handling of these permutations essential for the iterative Markovian fitting (IMF) algorithm.

### B.2 Sub-optimal vs. Optimal Permutation

In practice, the key insight is that the likelihood difference of graph vectors between the optimal permutation $\sigma^{*}$ and sub-optimal permutations $\sigma$ is substantial, allowing us to neglect sub-optimal permutations.
The design choice of the reference process $\mathbb{Q}$, particularly its signal to noise ratio $\frac{\bar{\alpha}(\tau)}{\bar{\alpha}(0)}\gg 0$ (see Equation 8), ensures the transitions preserving the initial states are far more probable.
As a result, the likelihood for sub-optimal permutations, where node and edge correspondences are mismatched, is expected to be significantly lower than for the optimal permutation.
The reason is that even a small mismatch in node alignment between $\sigma^{*}$ and $\sigma$ can cause a large number of edge state mismatches, leading to a dramatic decrease in the total likelihood.

Thus, the likelihood $p(\sigma^{*}\mathbf{G}^{\prime}|\mathbf{G})\gg p(\sigma\mathbf{G}^{\prime}|
\mathbf{G})$ and the contribution of sub-optimal permutation in calculating $p(\mathcal{G}^{\prime}|\mathcal{G})$ is negligible.
It is sufficient that focus solely on the optimal permutation $\sigma^{*}$ for the likelihood computation, as the probability of sub-optimal arrangements is effectively zero in practical terms.

The prioritization of optimal permutation could be also utilized in the construction of reciprocal measures.
Suppose we have only one graph for initial and terminal, where the one is $\mathcal{G}$ and $\mathcal{G}^{\prime}$, respectively.
Without loss of generality, we assume an even distribution over graph vectors $\mathbf{G}\in\mathcal{G}$, and then consider permutations over $\sigma\mathbf{G}^{\prime}\in\mathcal{G}^{\prime}$.
Note that $p(\sigma\mathbf{G}^{\prime}|\mathbf{G})=p(\sigma^{\prime}\sigma\mathbf{G}^{
\prime}|\sigma^{\prime}\mathbf{G})$.
According to the static SB solution, the graph vectors $\sigma\mathbf{G}^{\prime}\in\mathcal{G}^{\prime}$ are distributed according to the likelihood $p(\sigma\mathbf{G}^{\prime}|\mathbf{G})$ for each graph vector $\mathbf{G}\in\mathcal{G}$, where the optimal permutation is selected most frequently.
In this context, neglecting sub-optimal permutation can be viewed as the static OT solution.

As the original SB problem is formulated with the dynamics over graph ”vector” domain, arranging the graph vectors is a part of the SB problem.
However, by relying on the optimal permutation, we can reduce complexity of the SB problem aroused by the graph vectors alignment.
This strategy ensures the transition from a graph vector to the optimally permuted graph vector, resulting in more reliable convergence of the iterative algorithm.

### B.3 Standard formulation of Graph matching problem

In this section we will discuss about the standard formulation of the graph matching problem and how the task of finding the optimal permutation $\sigma^{*}$, which minimizes the negative log-likelihood (NLL) $-\log{p(\sigma\mathbf{G}^{\prime}|\mathbf{G})}$, can be framed as a graph matching problem.

Consider two graph vectors $\mathbf{G}=(\mathbf{V},\mathbf{E})$ and $\mathbf{G}^{\prime}=(\mathbf{V}^{\prime},\mathbf{E}^{\prime})$ along with the cost function $\mathbf{c}^{V}$ and $\mathbf{c}^{E}$, which define the cost of mapping nodes and edges, respectively.
The binary assignment matrix $\mathbf{X}\in\{0,1\}^{n\times n^{\prime}}$ represents the matching between nodes, where $n$ and $n^{\prime}$ denote the number of nodes of $\mathbf{G}$ and $\mathbf{G}^{\prime}$, respectively.
If $v_{i}\in\mathbf{V}$ matches $v_{a}^{\prime}\in\mathbf{V}^{\prime}$, then $\mathbf{X}_{i,a}=1$, while all other entries for node $v_{i}$ are zero.

The total node matching cost is given by $\sum_{\mathbf{X}_{i,a}=1}\mathbf{c}^{V}(v_{i},v_{a}^{\prime})$, where $\mathbf{c}^{V}(v_{i},v_{a}^{\prime})$ denotes the cost of matching from $v_{i}\in\mathbf{V}$ to $v_{a}^{\prime}\in\mathbf{V}^{\prime}$.
Similarly, the total edge matching cost is $\sum_{\mathbf{X}_{i,a}=1,\mathbf{X}_{j,b}=1}\mathbf{c}^{E}(e_{ij},e_{ab}^{
\prime})$, where $\mathbf{c}^{E}(e_{ij},e_{ab}^{\prime})$ is the cost of matching from $e_{ij}\in\mathbf{E}$ to $e_{ab}^{\prime}\in\mathbf{E}^{\prime}$.

We define the cost matrix $\mathbf{A}\in\mathbb{R}^{nn^{\prime}\times nn^{\prime}}$, where the diagonal components $\mathbf{A}_{ia,ia}=\mathbf{c}^{V}(v_{i},v_{a}^{\prime})$ represent node costs and the off-diagonal components $\mathbf{A}_{ia,jb}=\mathbf{c}^{E}(e_{ij},e_{ab}^{\prime})$ represent edge costs.
By flattening the assignment matrix to be $\mathbf{x}\in\{0,1\}^{nn^{\prime}}$, we can write the total cost function $f$ as

$$ $\displaystyle f(\mathbf{x})$ $\displaystyle=\sum_{\mathbf{x}_{ia}=1,\mathbf{x}_{jb}=1}\mathbf{c}^{E}(e_{ij}, e_{ab}^{\prime})+\sum_{\mathbf{x}_{ia}=1}\mathbf{c}^{V}(v_{i},v_{a}^{\prime}),$ $\displaystyle=\mathbf{x}^{\top}\mathbf{A}\mathbf{x}.$ $$

Thus, the graph matching problem is find the optimal assignment $\mathbf{x}$ that minimizing cost $f(\mathbf{x})$, which is written in formal:

$$ $\displaystyle\mathbf{x}^{*}=\arg\min_{\mathbf{x}}\mathbf{x}^{\top}\mathbf{A} \mathbf{x},$ (12) $\displaystyle s.t.\sum_{i=1}^{n}\mathbf{x}_{ia}\leq 1,\sum_{a=1}^{n^{\prime}}x _{ia}\leq 1,$ $$

where the inequalities become equalities if $n=n^{\prime}$.
This formulation corresponds to the well-known quadratic assignment problem (QAP), which is NP-hard.

The NLL $-\log{p(\mathbf{G}^{\prime}|\mathbf{G})}$ can be decomposed by the sum of the NLLs of nodes and edges:

$$ $-\log{p(\mathbf{G}^{\prime}|\mathbf{G})}=\sum_{\delta_{ia}=1,\delta_{jb}=1}- \log{P_{0:\tau}^{E}(e_{ij},e_{ab}^{\prime})}+\sum_{\delta_{ia}=1}-\log{P_{0: \tau}^{V}(v_{i},v_{a}^{\prime})}.$ $$

By interpreting the NLLs of nodes and edges as cost function $\mathbf{c}^{V}$ and $\mathbf{c}^{E}$, respectively, the problem of finding the optimal permutation $\sigma^{*}$ can be formulated as a QAP, as in Equation 12.

### B.4 Solution method

Exact solution methods for the QAP, such as mixed integer programming (MIP), requires combinatorial optimization, which incurs prohibitive computational costs Sahni & Gonzalez ().
Many accelerated algorithms adopt branch-and-bound strategies that utilize bounds of objective functions Anstreicher (); Gilmore (); Loiola et al. (); Xia (), reducing the exploration space.
However, the combinatorial optimization cannot be avoidable in these strategy.
Alternatively, continuous relaxation methods Cho et al. (); Leordeanu & Hebert () solve the QAP in Equation 12 using a continuous vector $\mathbf{x}$, bypassing combinatorial explorations.
However, the continuous relaxation yields non-binary vectors, which require a discretization step, potentially introducing errors.

In our work, we compared two continuous relaxation methods: spectral method (SM) and maximum polling method (MPM), followed by the Hungarian algorithm for post-discretization.
Both methods iteratively update the continuous version of assignment vector $\mathbf{x}$ as:

$$ $\mathbf{x}_{k+1}=\frac{\mathbf{x}_{k}-\mathbf{A}\odot\mathbf{x}_{k}\varepsilon }{v},$ (13) $$

where $\odot$ denotes matrix multiplication for SM and max-pooled matrix multiplication for MPM, and $\varepsilon$ denotes step size, and $v$ denotes the normalizer.
Specifically, the $\mathbf{A}\odot\mathbf{x}$ for SM is described as

$$ $(\mathbf{A}\odot\mathbf{x})_{ia}=\mathbf{x}_{ia}\mathbf{A}_{ia:ia}+\sum_{j\in \mathcal{N}_{i}}\sum_{b\in\mathcal{N}_{a}}\mathbf{x}_{jb}\mathbf{A}_{ia;jb}.$ $$

In the MPM case, the operation is defined as

$$ $(\mathbf{A}\odot\mathbf{x})_{ia}=\mathbf{x}_{ia}\mathbf{A}_{ia:ia}+\sum_{j\in \mathcal{N}_{i}}\max_{b\in\mathcal{N}_{a}}\mathbf{x}_{jb}\mathbf{A}_{ia;jb}.$ $$

In practice, we modified the cost matrix $\mathbf{A}$ slightly for efficiency.
To accelerate the process, we neglect all edges corresponding to dummies in $\mathbf{G}^{\prime}$, which significantly reduces computational cost, making $|\mathcal{N}_{a}|$ scale linearly with $n^{\prime}$.
Additionally, a small Gaussian perturbation was applied to $\mathbf{A}$, slightly altering the minima in the continuous vector space.
We solved each QAP ten times to compensate the effects of randomness, improving performance with reasonable computational costs (see Section B.5).

### B.5 Performance of graph matching algorithm

In this section, we compare the performance of the following algorithms under different conditions: (1) SM algorithm, (2) MPM algorithm, and (3) MPM algorithm with randomness.
The hyper-parameters $\bar{\alpha}(\tau)/\bar{\alpha}(0)=0.3$ for the reference process $\mathbb{Q}$.

To evaluate the effectiveness of QAP solvers on molecular graphs, we selected 100 molecules from the ZINC test set.
In all experiments, the source molecule was treated as fully connected (with dummy types), while the target molecule retained only its original edges.
We performed up to 1,000 iterations of updates according to Equation 13 with a specified tolerance.

To assess the performance of the algorithms, we permuted the molecular graphs and then tested whether the original indices could be recovered, where optimality is achieved with the inverse permutation.
For this task, we evaluated the exact matching ratio, which indicates the proportion of cases where the original indices were successfully recovered.
Due to the inherent symmetry in molecular graphs, multiple optimal permutations may exist. Therefore, instead of relying solely on exact index recovery, we base the success criterion on the objective function value.
If the Negative Log-Likelihood (NLL) error falls below 1e-2, the solution is considered successful.

For different molecule pairs where no optimal permutation is available, which is common in practical scenarios, exact solvers become computationally prohibitive even with relatively small molecular graphs.
As such, we evaluated the performance of the QAP solvers by measuring the reduction in NLL between the initial and optimized permutations.

We conducted experiments with various hyperparameters, and the results are illustrated in Table 4.
All algorithms successfully found exact matches for the same molecules, but the MPM algorithm performed better than the SM algorithm when applied to different molecule pairs.
For all subsequent experiments, we adopted the MPM algorithm with randomness in the cost matrix as the QAP solver.

**Table 4: Comparison of different algorithms based on configuration parameters. NLL drop represents the improvement in likelihood between randomly paired molecules, while exact match reflects the percentage of permutations successfully recovered after applying a random permutation to the molecules. $\uparrow$ and $\downarrow$ denote higher and lower values are better, respectively.**
| Algorithm | Configuration | NLL drop($\uparrow$) | Exact matching($\uparrow$) |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Pooling | Tolerance | Precision | Max # iterations | Noise coefficient | # trials |  |  |  |
| (1) SM | Sum | 1e-9 | FP64 | 1000 | 0 | 1 | 2.251 | 100 % |
| (2) MPM | Max | 1e-5 | FP32 | 1000 | 0 | 1 | 13.256 | 100 % |
| (3) MPM + Randomness | Max | 1e-5 | FP32 | 1000 | 1e-6 | 10 | 13.256 | 100 % |
| (4) Our setting | Max | 1e-4 | FP32 | 2500 | 1e-6 | 10 | 14.324 | 100 % |

### B.6 Relation to graph edit distance

The graph edit distance (GED) is a widely used and flexible metric for measuring the dissimilarity between two graphs.
It is defined as the minimum cost required to transform one graph into another through a sequence of unit operations.
Each unit operation can be a removal, substitution, or insertion, and can be applied to either node or edges.
The total cost of the transformation is the sum of the costs assigned to these components.

According to Bougleux et al. (), finding the minimal-cost transformation between two graphs is equivalent to the QAP problem with an associated cost function.
Though both problems are NP-hard, the equivalence is meaningful in that there are numerous approximation algorithms for the QAP problem that operate within polynomial time.

The basic idea for re-formulation into QAP problem involves the introduction of dummy nodes and dummy edges, where removal (insertion) operations could be replaced by substitution into (from) dummies.
The cost function of unit operation is defined as $\mathbf{c}^{V}(v_{i},v_{a}^{\prime})$ and $\mathbf{e}^{E}(e_{ij},e_{ab}^{\prime})$ for node replacement and edge replacement, respectively.
Let $\alpha=\bar{\alpha}(\tau)/\bar{\alpha}(0)$, the cost replacement is defined as:

$$ $\displaystyle\mathbf{c}^{V}(v_{i},v_{a}^{\prime})=\begin{cases}-\log{\frac{(d^ {V}-1)\alpha+1}{d^{V}}}&\text{if }v_{i}=v_{a}^{\prime}\\ -\log{\frac{-\alpha+1}{d^{V}}}&\text{otherwise,}\end{cases}$ $\displaystyle\mathbf{c}^{E}(e_{ij},e_{ab}^{\prime})=\begin{cases}-\log{\frac{( d^{E}-1)\alpha+1}{d^{E}}}&\text{if }e_{ij}=e_{ab}^{\prime}\\ -\log{\frac{-\alpha+1}{d^{E}}}&\text{otherwise,}\end{cases}$ $$

where $d^{V}$ and $d^{E}$ denote the cardinality of $\mathcal{X}^{V}$ and $\mathcal{X}^{E}$, respectively.
Similar to Equation 12, we can formulate it to a quadratic problem with the objective function $f$ as:

$$ $\displaystyle f(\mathbf{x})$ $\displaystyle=\sum_{\mathbf{x}_{ia}=1,\mathbf{x}_{jb}=1}\mathbf{c}^{E}(e_{ij}, e_{ab}^{\prime})+\sum_{\mathbf{x}_{ia}=1}\mathbf{c}^{V}(v_{i},v_{a}^{\prime}),$ $\displaystyle=\sum_{\mathbf{x}_{ia}=1,\mathbf{x}_{jb}=1}-\log{P_{0:\tau}^{E}(e _{ij},e_{ab}^{\prime})}+\sum_{\mathbf{x}_{ia}=1}-\log{P_{0:\tau}^{V}(v_{i},v_{ a}^{\prime})},$ $\displaystyle=p(\sigma\mathbf{G}^{\prime}|\mathbf{G}),$ $$

where, the $\sigma$ is the graph permutation associated to the assignment vector $\mathbf{x}$.

However, the GED is not same to the NLL.
Note that the edit cost functions penalize every operations with same cost $-\log{\frac{-\alpha+1}{d}}$ except for the identity operation.
Though the cost of identity operation is lesser than the others, the optimal transformation would not contains any identity operation.
In real, the cost of identity operation does not affect the optimal edit path, implying that GED is not equal but proportion to the NLL, where the difference proportional to the number of the nodes and edges that are equal under the optimal graph matching $\sigma^{*}$.
Though the GED is not exactly same to the NLL, the problem is equivalent to the graph matching problem.

This observation provides a clear interpretation of the underlying dynamics of the SB problem, revealing that the associated OT cost is effectively the GED.
Therefore, solving the SB problem can be understood as finding the OT plan between graph distributions, where the transport cost is defined by the GED.

## Appendix C Implementation details

### C.1 Parameterization

In this section, we briefly describe the neural network parameterization and the practical training loss.
According to Equations 4 and 5, the neural network approximates the generator $A_{t}^{\Lambda_{\cdot|0}}(x,y)$ and $\tilde{A}_{t}^{\Lambda_{\cdot|\tau}}(y,x)$, which are formulated as conditional expectation of $A(x,y;z)$ and $\tilde{A}(y,x;z)$, respectively (see Lemma A.5).

The target generator takes the form:

$$ $A_{s}(x,y;z)=\begin{cases}A_{s}(x,y)\frac{P_{s:\tau}(y,z)}{P_{s:\tau}(x,z)},& \text{if }x\neq y\\ -\sum_{u\neq x}A_{s}(x,u)\frac{P_{s:\tau}(u,z)}{P_{s:\tau}(x,z)},&\text{if }x= y\end{cases}$ (14) $$

, where $A_{s}$ and $P_{s:\tau}$ are tractable attributes of $\mathbb{Q}$.
Let the the transition probability matrix be $P(s,\tau):[0,\tau]^{2}\to\mathbb{R}^{|\mathcal{X}|\times|\mathcal{X}|}$, with $P_{s:\tau}(x,y)=\mathbf{e}_{x}^{\top}P(s,\tau)\mathbf{e}_{y}$.
Thus, the neural network predicts the distribution of $z$ given $X_{t}$, denoted as $z_{\theta}(t,X_{t})\in\mathbb{R}^{|\mathcal{X}|}$, and our parameterization choice is as follows:

$$ $A_{s}^{M^{\theta}}(X_{t},y;t,X_{t})=A_{s}(x,y)\frac{\mathbf{e}_{y}P(s,\tau)z_{ \theta}(t,X_{t})}{\mathbf{e}_{x}P(s,\tau)z_{\theta}(t,X_{t})},$ $$

for $X_{t}\neq y$.
The time reverse generator is also defined similarly.

Though the loss formulation was defined as continuous manner, we approximate it with the discretizations.
Firstly, we will replace $A_{t}(X_{t},X_{t};z)-A_{t}^{M^{\theta}}(X_{t},X_{t})$ with as follows:

$$ $\displaystyle A_{t}(X_{t},X_{t};z)-A_{t}^{M^{\theta}}(X_{t},X_{t})$ $\displaystyle\approx\frac{1}{\Delta t}(1+A_{t}(X_{t},X_{t};z)\Delta t)\log{ \frac{1+A_{t}(X_{t},X_{t};z)\Delta t}{1+A_{t}^{M^{\theta}}(X_{t},X_{t})\Delta t }},$ $\displaystyle\approx\frac{1}{\Delta t}P_{t:t+\Delta t}^{\mathbb{Q}_{\cdot|\tau =z}}(X_{t},X_{t})\log\frac{P_{t:t+\Delta t}^{\mathbb{Q}_{\cdot|\tau=z}}(X_{t}, X_{t})}{P_{t:t+\Delta t}^{M^{\theta}}(X_{t},X_{t})}.$ $$

Similarly,

$$ $\displaystyle A_{t}(x,y;z)\log{\frac{A_{t}(x,y;z)}{A_{t}^{M^{\theta}}(x,y)}} \approx\frac{1}{\Delta t}P_{t:t+\Delta t}^{\mathbb{Q}_{\cdot|\tau=z}}(x,y)\log \frac{P_{t:t+\Delta t}^{\mathbb{Q}_{\cdot|\tau=z}}(x,y)}{P_{t:t+\Delta t}^{M^{ \theta}}(x,y)}.$ $$

Thus the discretized loss function corresponding to Equation 4 is:

$$ $L(\theta)=\sum_{t_{i}}\frac{1}{\Delta t}\mathbb{E}_{\Lambda_{t_{i},\tau}}\left [\sum_{y}P_{t:t+\Delta t}^{\mathbb{Q}_{\cdot|\tau=z}}(x,y)\log\frac{P_{t:t+ \Delta t}^{\mathbb{Q}_{\cdot|\tau=z}}(x,y)}{P_{t:t+\Delta t}^{M^{\theta}}(x,y) }\right].$ $$

The time reversed loss can be discretized as same way.

### C.2 Neural Network Parameterization and Hyperparameter settings

DDSBM and DBM.
Our neural network parameterization is based on Vignac et al. (), which uses a graph transformer network (Dwivedi & Bresson, ).
In short, it takes as input a noisy graph $\mathbf{G}_{t}=(\mathbf{V}_{t},\mathbf{E}_{t})$ and predicts a distribution over the target graphs.
Structural and spectral features are also used as inputs to improve the expressivity of neural networks.
We refer the reader to Vignac et al. () for more details.

For noise scheduling, we employ a slightly different strategy than Vignac et al. ().
While Vignac et al. () uses a cosine schedule for $\bar{\alpha}_{t}$, we implement a symmetric scheduling of $\alpha_{t}$ by incorporating $\alpha_{\text{min}}$, as defined below:

$$ $\alpha(t)=\frac{\partial_{t}\bar{\alpha}(t)}{\bar{\alpha}(t)}=\cos\left(\frac{ t/\tau+s}{1+s}\cdot\frac{\pi}{2}\right)^{2}\cdot(1-\alpha_{\text{min}})+\alpha _{\text{min}},$ (15) $$

for $0\leq t\leq\tau/2$, with $\alpha(t)=\alpha(\tau-t)$ for the remaining half of the schedule.
For instance, with 100 diffusion steps, $\bar{\alpha}(\tau)\approx 0.95$ when $\alpha_{\text{min}}=0.999$, and $\bar{\alpha}(\tau)\approx 0.90$ when $\alpha_{\text{min}}=0.99795$.

Figure: Figure 3: Plot of $\alpha(t)$ and $\bar{\alpha}(t)$ as functions of timestep $t$ where $\alpha_{\text{min}}=0.999$
Refer to caption: extracted/5895698/Figures/source/alpha_schedule.png

We trained DDSBM models with IMF iterations, and DBM models, an one-directional variant of DDSBM with fixed joint molecular pairs, were trained with the same number of gradient updates to ensure consistency.
Both DBM and DDSBM reported in this work were trained using four RTX A4000 GPUs.
The detailed hyperparameters for training are shown in Table 5.

**Table 5: Training hyperparmeters of DBM and DDSBM.**
| Task | Model | diffusion steps | $\alpha_{\text{min}}$ | epoch | SB iterations |
| --- | --- | --- | --- | --- | --- |
| ZINC250k | DBM | 100 | 0.999 | 1800 | – |
| ZINC250k | DDSBM | 100 | 0.999 | 300 | 6 |
| Polymer | DBM | 100 | 0.999 | 1250 | – |
| Polymer | DDSBM | 100 | 0.999 | 250 | 5 |

Graph-to-Graph Translation.
For HierG2G and AtomG2G, we used the default settings provided on the official GitHub repository.
Both models were trained until maximum epochs by default with a single RTX A4000 GPU.

## Appendix D Supplementary results

### D.1 Data processing for ZINC250K dataset experiment

In Section 5.2, we experimented with the standard ZINC250K dataset.
In our molecular graph representation, $\mathbf{G}=(\mathbf{V},\mathbf{E})$, the node vectors $\mathbf{V}=(v^{(i)})_{i}$ represent the atomic types, and the edge vectors $\mathbf{E}=(e^{(ij)})_{ij}$ represent bond orders.
Due to this implementation choice, node features other than the atomic type cannot be represented, causing occasional failures in decoding into molecules.
Still, to focus on the SB problem itself, we filtered out molecules whose graph representations are not directly converted into molecules by the RDKit package.
Among the node features, formal charge and explicit hydrogen information is crucial in decoding process since, without them, each atom’s valency cannot be inferred so that corresponding molecule cannot be uniquely determined.
Thus, the following criteria were applied to filter molecules from the ZINC250K dataset: (1) all atoms do not possess a formal charge and (2) all aromatic atoms do not have an explicit hydrogen.

### D.2 Examples for generated molecules

To show the difference between molecules generated from DDSBM and others, we visualized some selected examples for the ZINC250K and polymer datasets, respectively (see Figures 4 and 5).

Figure: Figure 4: Visualization of molecules generated by DDSBM, DBM, and HierG2G compared to the source molecule.
Refer to caption: extracted/5895698/Figures/source/ZINC_gen.png

Figure: Figure 5: Visualization of molecules generated by DDSBM and DBM with the source molecule. The samples generated by DBM and DDSBM were selected from the molecules predicted to have a blue color with GAP values in the range of 2.56–2.75 eV.
Refer to caption: extracted/5895698/Figures/source/Figure_POLYMER.png

Figure: Figure 6: Density plots comparing the distributions of three molecular properties ($\log{P}$, $\Delta$QED, and $\Delta$SAscore) across four different molecular generation methods: AtomG2G, HierG2G, DBM, and DDSBM.
Refer to caption: extracted/5895698/Figures/source/Figure_Distributions.png

## Peer Reviews

confidence: 2
contribution: 3
initial_rating: 6
presentation: 3
questions: - What is the Hungarian algorithm ? You mention it in section 3.3. I'd appreciate if you could drop a reference or give some details in appendix.
- The baselines outperform your approach in terms of Validity. Even though the results on FCD and NLL indicate that your method is stronger, I'd appreciate if you could comment on this metric.
soundness: 3
strengths: - The Schrödinger bridge problem is well studied in the continuous-state diffusion literature, but has, to the best of my knowledge, not yet been applied to the discrete setting. The paper nicely bridges this gap, hence proposing a valuable contribution to the discrete diffusion literature.
- The paper looks sound and technically strong.
- Empirical results seem to indicate that this approach outperforms previous baselines by a large margin.
summary: This work presents Discrete Diffusion Schrödinger Bridge Matching (DDSBM), a framework that adapts continuous-time Markov chains (CTMCs) to address the Schrödinger Bridge (SB) problem within high-dimensional discrete spaces. Using the Iterative Markovian Fitting (IMF) technique, DDSBM enables optimal graph modifications by minimizing structural changes, which is particularly valuable for molecular optimization tasks in drug and material discovery. In this context, the framework aligns with graph edit distance (GED), allowing efficient property-driven modifications while preserving the molecule's structure. Experiments show DDSBM achieves minimal structural shifts and successfully maintains desirable molecular properties, outperforming traditional graph translation methods.
weaknesses: - The choice of the FCD metric the evaluate the graph structure distributions is surprising, since it is typically used to assess distribution learning from a chemical perspective. Why not use NSPDK for example ?
- I’d like to see some error bars on your results.
- The method description is quite technical, and might be hard to understand for the average graph machine learning practioner. It could be worth it to include a short tutorial on the Schrödinder Bridge problem to ease the reader understanding.

---

confidence: 3
contribution: 2
initial_rating: 3
presentation: 2
questions: 1. What is the relevance and difference between the proposed method and existing discrete SB models in Weakness?

2. How does the performance of molecular optimization compare with discrete flow matching models?
soundness: 3
strengths: Strong performance compared to diffusion bridge model.
summary: This paper introduces discrete diffusion Schrödinger bridge model using continuous-time Markov chains. The proposed model was validated on molecular optimization problem.
weaknesses: Overall, this paper failed to distinguish its own contributions from existing work and to compared itself with relevant works. First of all, the proposed discrete diffusion Schrödinger bridge matching (DDSBM) highly relies on continuous-time Markov chains (CMTC). However, this paper failed to refer flow matching models [1, 2] which introduce CMTC for discrete data domains. 

[1] Campbell, Andrew, et al. "Generative Flows on Discrete State-Spaces: Enabling Multimodal Flows with Applications to Protein Co-Design." Forty-first International Conference on Machine Learning.

[2] Gat, Itai, et al. "Discrete flow matching." arXiv preprint arXiv:2407.15595 (2024).

And, this is not the first paper that solving Schrödinger bridge (SB) problem in discrete state spaces. There are works for discrete SB [3, 4]. However, the relevance and difference between the proposed method were not discussed.

[3] Harchaoui, Zaid, Lang Liu, and Soumik Pal. "Asymptotics of discrete Schrödinger bridges via chaos decomposition." Bernoulli 30.3 (2024): 1945-1970.

[4] Chow, Shui-Nee, et al. "A discrete Schrodinger bridge problem via optimal transport on graphs." calculus of variations 20.33 (2021): 34.

---

confidence: 4
contribution: 3
initial_rating: 5
presentation: 2
questions: How significantly does the choice of graph-matching algorithm impact the overall performance of the method?
If it has a substantial influence, there is another reason that experiments for some simpler modalities without such complex matching may help to clarify the contribution.
soundness: 3
strengths: 1. This method adapts diffusion Schrodinger bridge matching which opens a new avenue for this task.
2. The experimental setting for molecular optimization is a very relevant downstream task for the discrete domain, and it builds up some baselines for that task.
3. The work provides corresponding theoretical guarantees.
summary: Adapt diffusion Schrodinger bridge matching to the discrete domain and instantiate it through molecular optimization
weaknesses: 1. Globally, the structure is complete but the writing is not easy to follow. The writing can be enhanced at least in the following way:
* The square brackets around the citation should not appear twice.
* The name of section 2.2 'SOLUTION METHOD' needs to be optimized.
* The supplementary material also contains the body part of the paper.
* The structure section 2/3 is a bit confusing:
  * There should be an introduction to the continuous methods (IPF, IMF) as prior knowledge since the method highly dependent on them, then a transition to the method for the discrete domain. Here, section 2 starts by discussing the discrete domain and introducing the definitions together, which makes it a bit confusing which parts come from previous work and which part comes from this work.
  *  Consequently, the insights about the transition from continuous domain to discrete domain are missing.

2. Under the claim of proposing a discrete model for DSBM, the results only contain 2 molecule datasets with few supportive ablation experiments being given. It will support better the contribution of this work to give toy results for examples except for molecule optimization such as toy datasets with discrete features, or results on other types of graphs, or more ablations, or pure generation results (as in continuous SB papers), or more visualizations of the resulting optimization chain. No need to have them all, but similar supportive evidence would definitely make the experiments more complete.

Meta-review: The authors introduce Discrete Diffusion Schrödinger Bridge Matching (DDSBM), a new method for transforming data in discrete spaces like graphs using continuous-time Markov chains. DDSBM improves key properties of molecules with minimal changes to their structure, preserving other important features.

Strengths: see comments by reviewers.

Weaknesses: The main major weakness after rebuttal was raised by reviewer 7pXr, who requested comparisons with the latest discrete flow matching methods. 

Overall the paper was well received the rest of the reviewers.

Decision: Accept (Poster)