Feynman-Kac Operator Expectation Estimator

## Abstract

Abstract The Feynman-Kac Operator Expectation Estimator (FKEE) is an innovative method for estimating the target Mathematical Expectation 𝔼 X ∼ P ⁢ [ f ⁢ ( X ) ] subscript 𝔼 similar-to 𝑋 𝑃 delimited-[] 𝑓 𝑋 \mathbb{E}_{X\sim P}[f(X)] blackboard_E start_POSTSUBSCRIPT italic_X ∼ italic_P end_POSTSUBSCRIPT [ italic_f ( italic_X ) ] without relying on a large number of samples, in contrast to the commonly used Markov Chain Monte Carlo (MCMC) Expectation Estimator. FKEE comprises diffusion bridge models and approximation of the Feynman-Kac operator. The key idea is to use the solution to the Feynmann-Kac equation at the initial time u ⁢ ( x 0 , 0 ) = 𝔼 ⁢ [ f ⁢ ( X T ) | X 0 = x 0 ] 𝑢 subscript 𝑥 0 0 𝔼 delimited-[] conditional 𝑓 subscript 𝑋 𝑇 subscript 𝑋 0 subscript 𝑥 0 u(x_{0},0)=\mathbb{E}[f(X_{T})|X_{0}=x_{0}] italic_u ( italic_x start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT , 0 ) = blackboard_E [ italic_f ( italic_X start_POSTSUBSCRIPT italic_T end_POSTSUBSCRIPT ) | italic_X start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = italic_x start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ] . We use Physically Informed Neural Networks (PINN) to approximate the Feynman-Kac operator, which enables the incorporation of diffusion bridge models into the expectation estimator and significantly improves the efficiency of using data while substantially reducing the variance. Diffusion Bridge Model is a more general MCMC method. In order to incorporate extensive MCMC algorithms, we propose a new diffusion bridge model based on the Minimum Wasserstein distance. This diffusion bridge model is universal and reduces the training time of the PINN. FKEE also reduces the adverse impact of the curse of dimensionality and weakens the assumptions on the distribution of X 𝑋 X italic_X and performance function f 𝑓 f italic_f in the general MCMC expectation estimator. The theoretical properties of this universal diffusion bridge model are also shown. Finally, we demonstrate the advantages and potential applications of this method through various concrete experiments, including the challenging task of approximating the partition function in the random graph model such as the Ising model.

## 1 Introduction

### 1.1 Background

Markov Chain Monte Carlo (MCMC) is a prevalent computational method used in fields such as statistics, machine learning, and computational science. It is primarily applied for sampling from complex distributions, Bayesian inference, and optimization[, ]. MCMC algorithms are typically divided into two categories: those that sample from the target distribution and those that estimate the statistical characteristics of the target distribution, such as the expectation. While alternative methods to traditional MCMC samplers, like generative models [] and diffusion models [], have been explored, MCMC remains the standard for expectation estimation. However, traditional MCMC estimators, based on the law of large numbers (LLN) and the ergodic theorem of Markov chains (ETMC), face limitations in data efficiency and impose complex constraints on the distribution $P$ and performance function $f$. Thus, developing algorithms that overcome these limitations by integrating modern sampling methods with deep learning is of great importance.

### 1.2 Motivation

Advantages of MCMC algorithm:
MCMC algorithms are effective for sampling from target distributions and are accompanied by two types of expectation estimators. The first type, based on LLN, uses averages from multiple samples at the terminal time of the Markov chain. The second type, based on ETMC, averages values along the path of the Markov chain. These methods utilize statistical principles effectively, particularly for high-dimensional distributions, mitigating the curse of dimensionality in integral approximations.

Disadvantages of MCMC algorithm: Despite their advantages, MCMC algorithms are not optimal for expectation estimation. The efficiency of MCMC estimators depends on the distribution $P$ and the function $f$. Different MCMC algorithms are required for different $P$, and due to burn-in periods, MCMC algorithms often waste many points. Additionally, the sample size $N$ shoule be large enough to achieve accurate estimates, leading to variances on the order of $\mathcal{O}(\sqrt{N})$. Quasi-Monte Carlo methods offer variances on the order of $\mathcal{O}(N^{\frac{1}{2}+\delta})$ with $\delta\leq\frac{1}{2}$ [], but this diminishes efficiency and introduces bias. Error probabilities can be estimated through concentration inequalities [], but these depend on the Lipschitzian norm of $f$.

Estimating mathematical expectations is crucial in both machine learning and statistics. A unified expectation estimator is theoretically and practically significant. Therefore, we focus on two essential questions:

(A) Is it possible to unify most existing MCMC algorithms into a cohesive framework to create a universal sampler for expectation estimation?

(B) How can we develop a universal expectation estimator that leverages samples from universal samplers for accurate expectation estimates without relying on post-processing or specialized methods?

For the first question: We propose the following solutions. we know the Markov model is determined by the transfer density of one step. The transition density function associated with the discrete Markov chain generated by the MCMC algorithm can be interpreted as the transition density function of a specific stochastic differential equation (SDE) of Markov properties. In this study, we refer to this SDE as the diffusion bridge model. This encompasses a broad class of SDEs that share identical transition densities with the Markov chains in the MCMC algorithm. The distribution of the terminals in such SDEs aligns with a predefined target distribution, which can take the form of discrete points or a probability density function. Moreover, the starting point of this SDE can be either arbitrary or fixed.

For the second question: In the context of the diffusion bridge model, we can view the expectation estimation problem as a decoding problem. It is easy to observe that MCMC is not the optimal decoding method because a substantial number of burn-in samples go to waste when estimating Mathematical expectations by using the MCMC algorithm. However, these samples harbor valuable information, specifically pertaining to the gradient information of the drift and diffusion coefficients along the paths derived from the SDE. We capitalize on this information by integrating it through the Physics-Informed Neural Network (PINN) approach [, , ]. This process, akin to approximating the Feynman-Kac Operator, is referred to as solving the Feynman-Kac model. Notably, this approximation is meshless and effectively overcomes the curse of dimensionality. By amalgamating different combinations of the aforementioned models, we derive the Feynman-Kac Operator Expectation Estimator (FKEE). In the framework of FKEE, the diffusion bridge model can extend the range of $P$, while the Feynman-Kac model can extend the range of $f$.

Our contributions can be summarized as follows:

- •
Establishing a Link Between Sampling Methods and High-Dimensional Partial Differential Equations: We establish a connection between widely used sampling techniques and the intricate realm of high-dimensional partial differential equations. Our innovative approach not only introduces a fundamental algorithm but also explores the synergy of merging two prominent algorithmic classes. We substantiate this through experimental analysis, demonstrating the potential of our approach as a bridge between these two disparate domains.
- •
Introducing a more versatilediffusion bridge model: We introduce a highly adaptable diffusion bridge model. This model not only allows for the specification of target distributions at terminal moments but also facilitates the reconstruction of the entire Markov chain. It can be employed in conjunction with the Feynman-Kac model for expectation estimation, as well as independently for resampling target distributions to estimate expectations.
- •
Expanding the Scope of Expectation Estimators: Our approach enhances the efficiency of Markov chains, requiring fewer assumptions and not relying on LLN or the Markov ergodic theorem. This method holds promise for approximating mathematical expectations with non-independent samples in real-world scenarios.

## 2 Main methodologies

### 2.1 Notations

Let $\mu_{t}$ represent the distribution of $X_{t}$ and $\hat{\mu}_{t}=\frac{1}{N}\sum^{N}_{i=1}\delta_{X^{i}_{t}}$ denote the empirical distribution of $X_{t}^{i},1\leq i\leq N$, where $\delta_{X}$ is the Dirac measure at $X$. $A_{i,j}$ denotes the elements of row $i$ and column $j$ of the matrix $A$, $diag(A)$ represents the diagonal matrix of matrix $A$ and $diag(A)_{i}$ denotes the i-th element on the main diagonal of the diagonal matrix of $A$. $\mathcal{C}^{2}$ denotes the space of continuous functions with second order derivatives. Denote $Y\sim P:=\mu^{*}$ by the target distribution.

### 2.2 Diffusion bridge model

The sampling methods mentioned in related work can be generalized into a common framework: for most MCMC sampling methods, we can consider using a Markov-type SDE as follows:

$$ $dX_{t}=b(X_{t},t)dt+\sigma(X_{t},t)dW_{t},\quad\quad X_{0}=x_{0},$ (1) $$

where $b$ : $\mathbb{R}^{d}\times[0,T]\to\mathbb{R}^{d}$ is a vector-valued function , $\sigma$ : $\mathbb{R}^{d}\times[0,T]\to\mathbb{R}^{d\times d}$ is a matrix-valued function, and $\{W_{t}\}_{t\geq 0}$ is a Brownian motion taking values on $\mathbb{R}^{d}$.

For a non-stationary diffusion model, these coefficients $b$ and $\sigma$ must satisfy certain regularity conditions to ensure the existence and uniqueness of a strong solution. For diffusion models with a stationary distribution, the uniqueness of the stationary distribution must hold. For a given distribution $P$ or a set of discrete points from $P$, we encode the information of $P$ into $(X_{0},b,\sigma,T)$. This encoding algorithm should ensure that the distance between $\mu_{T}$ and $P$ is sufficiently small. The encoding loss to be minimized is:

$$ $\mathcal{L}_{e}=\mathcal{W}_{2}(\mu_{T},\mu^{*}),$ (2) $$

where $\mu_{T}$ is the measure of solution of SDE (1) and $\mathcal{W}_{2}$ is 2-Wasserstein distance. The encoding loss $\mathcal{L}_{e}$ consists of two main components: structural loss and discretization loss. Structural loss is typically induced by the accuracy of $(b,\sigma)$, and discretization loss is usually due to the need for a sufficiently large $T$ and the numerical discretization of the SDE. According to the triangle inequality, the error can be decomposed as:

$$ $\mathcal{L}_{e}\leq\underbrace{\mathcal{W}_{2}(\mu_{T},\bar{\mu}_{T})}_{\text{ discretization loss}}+\underbrace{\mathcal{W}_{2}(\bar{\mu}_{T},\mu^{*})}_{ \text{structural loss}},$ (3) $$

where $\bar{\mu}_{T}$ is the distribution of numerical solution. The encoding loss depends on whether the distribution $P$ has an explicit density.

Specifically, if the density of $P$ is known, the structural loss has two components: the approximation error of $(X_{0},b,\sigma)$ itself and the asymptotic error. The asymptotic error exists only in certain methods when the target distribution is the stationary distribution, such as in Langevin MCMC, where the error of $(X_{0},b,\sigma)$ is zero and the asymptotic error decreases exponentially. However, discretization loss arises from the SDE solver.

If the density of $P$ is unknown, the constructed SDE will exhibit both types of losses. The core focus of the diffusion model is to minimize these two losses. The error in $(b,\sigma)$ is controlled by a specific loss function, while discretization loss is controlled by minimizing $T$ as much as possible and using a high-precision SDE solver.

In this paper we propose a diffusion bridge model that minimizes the encoding loss through the use of parameterized tuples $(X_{0},b,\sigma)$. This method is similar to the Neural SDE [, ]. Specifically, we use the following Neural SDE:

$$ $dX_{t}=b_{\theta_{1}}(X_{t},t)dt+\sigma_{\theta_{2}}(X_{t},t)dW_{t},\quad\quad X _{0,\theta_{3}}=x_{0,\theta_{3}},$ (4) $$

where $\mathcal{P}_{\theta}=(X_{0,\theta_{3}},b_{\theta_{1}},\sigma_{\theta_{2}})$ is a neural network, typically a multi-layer perceptron (MLP) with the tanh activation function. Here, the time $T$ and time step $h$ are given in advance for the SDE solver. Diffusion bridge model matching means that we use neural network methods to find the appropriate $(X_{0}^{*},b^{*},\sigma^{*})$ such that the distribution of $X_{T}$ at the moment $T$ is just the given target distribution $P$. We need to categorise the target distribution to determine the matching method. This depends on whether the target distribution has an explicit probability density function.

Encoding loss for diffusion bridge models

We examine the error of the diffusion bridge model. Unlike other design loss functions, we aim to control both errors in loss (3) simultaneously. Different loss functions correspond to different problems, necessitating the classification of the target distribution.

Only a few discrete observations: We propose a matching algorithm that deals with only a subset of discrete points from the target distribution $P$. Specifically, we employ a diffusion bridge model to parameterize $(X_{0,\theta},b_{\theta},\sigma_{\theta})$ using a neural network. Given the empirical distribution of the target $\hat{\mu}^{*}$, we simulate $N$ trajectories of Brownian motion and use the Euler-Maruyama method [] to obtain the solution $\bar{X}_{T}$. Subsequently, we match the obtained solutions to the given points and utilize the Wasserstein distance loss function:

$$ $\mathcal{P_{\theta}}^{*}=\mathop{\rm argmin}_{\mathcal{P_{\theta}}\in\Theta} \mathcal{W}_{2}(\hat{\bar{\mu}}_{T},\hat{\mu}),$ (5) $$

where $\hat{\bar{\mu}}_{T}$ (and $\hat{\mu}$ respectively) is the empirical distribution of independent identical copies of $\bar{X}_{T}$ (and $\mu^{*}$ respectively).
Given $T$ and $h$, we can estimate the discretization loss and control the structural loss through the Wasserstein distance loss. This method has two additional applications:

(i) Resample (Generate) samples: For a set of high-quality samples (not within the burn-in period of MCMC), this method can be used for resampling. By matching a diffusion bridge model to the given points, we can simulate the SDE to obtain more samples. High-quality samples can also be obtained through other methods, such as Perfect Sampling [].

(ii) Matching Markov chains and generating more Markov chains quickly: For trajectories $Y_{i}^{N}$ of $N$ independently run Markov chains obtained from MCMC algorithms $Y_{i}\sim\mu_{t_{i}}^{*}$ where $i\leq M$, we aim to find a set of $(X_{0},b,\sigma)$ such that $\bar{X}_{i}$ and $Y_{i}$ are close at $M$ moments in the sense of the Wasserstein distance. Here, $\bar{X}_{i}$ is the solution to the SDE defined by $(X_{0},b,\sigma)$. This can be achieved by optimizing the Wasserstein distance loss:

$$ $\mathcal{P_{\theta}}^{*}=\mathop{\rm argmin}_{\mathcal{P_{\theta}}\in\Theta} \sum_{i=1}^{M}\mathcal{W}_{2}(\hat{\bar{\mu}}_{t_{i}},\hat{\mu}_{t_{i}}^{*}),$ (6) $$

where $\hat{\bar{\mu}}_{t_{i}}$ and $\hat{\mu}_{t_{i}}^{*}$ are empirical distributions corresponding to $\bar{X}_{i}$ and $Y_{i}$ respectively. This process allows for matching either a segment or the entire Markov chain, potentially starting the matching process from a later moment to minimize reliance on points within the burn-in period.

Figure: Algorithm 1 Diffusion bridge model (DBM)

This algorithm is a simplified version of a more detailed one available in the Appendix 10.

In the following we present some theoretical results with the proofs given in the Appendix. We first provide an estimate for the discrete loss.

Theoretical results

###### Theorem 2.1 .

Assuming that $b$ and $\sigma$ are L-lipschitz functions and Linear growth. SDE solver is the Euler-Maruyama method ,we can obtain the following estimate:

$$ $\mathcal{W}_{2}(\mu_{T},\bar{\mu}_{T})\leq Ch^{\frac{1}{2}}\exp{(4L^{2}T)},$ (7) $$

where $C$ depends on $X_{0}$, but it is independent of $h$. We can pre-select suitable $T$ and $h$ to control this error.

###### Proof.

In the Appendix.
∎

For convenience we use the following notations to indicate that the measures depend on the parameter. $\mu_{T}:=\mu_{T}^{\mathcal{P_{\theta}}}$, $\bar{\mu}_{T}:=\bar{\mu}_{T}^{\mathcal{P_{\theta}}}$, $\hat{\mu}_{T}:=\hat{\mu}_{T}^{N,\mathcal{P_{\theta}}}$, $\hat{\bar{\mu}}_{T}:=\hat{{\bar{\mu}}}^{N,\mathcal{P_{\theta}}}_{T}$ and $\hat{\mu}_{t_{i}}:=\hat{\mu}_{t_{i}}^{N,\mathcal{P_{\theta}}}$, $\hat{\mu}_{t_{i}}:=\hat{\mu}_{t_{i}}^{N}$.$\hat{\mu}:=\hat{\mu}^{N}.$

In fact, we use the Minimal Wasserstein distance estimator, as the properties of this estimator have been outlined in []. We apply it to our specific problem to control the structural loss. We first introduce the following assumptions:

###### Assumption 2.2 .

The model is identifiable: there exists a unique $\mathcal{P_{\theta}}^{*}\in\Theta$ such that $\mathcal{W}_{2}(\mu^{\mathcal{P^{*}_{\theta}}},\nu)=\mathcal{W}_{2}(\mu^{*},\nu)$ for every $\nu$ and

$$ $\mathcal{P}^{*}_{\theta}={\mathop{\rm argmin}_{\mathcal{P_{\theta}}\in\Theta}} \mathcal{W}_{2}(\mu^{*},\bar{\mu}_{T}^{\mathcal{P_{\theta}}}).$ $$

This assumption ensures the existence of a deterministic parameter in the SDE.

###### Assumption 2.3 .

Data processes are sufficient: The data generation process error is satisfied $\mathcal{W}_{2}(\hat{\mu}^{N},\mu^{*})\to 0$, $\mathbb{P}$-almost surely, as $N\to\infty$.

###### Assumption 2.4 .

Continuity: The map $\mathcal{P_{\theta}}\mapsto\bar{\mu}_{T}^{\mathcal{P_{\theta}}}$ is continuous in the sense that $D(\mathcal{P}_{\theta}^{N},\mathcal{P_{\theta}})\to 0$ implies ${{\bar{\mu}}}^{\mathcal{P}_{\theta}^{N}}_{T}\stackrel{{\scriptstyle w}}{{\to}}
{{\bar{\mu}}}^{\mathcal{P_{\theta}}}_{T}$ as $N\to\infty$. $D$ is the metric of the parameter.

###### Assumption 2.5 .

Level boundedness: For some $\epsilon>0$
The set $B(\epsilon)=\left\{\mathcal{P_{\theta}}\in\Theta:\mathcal{W}_{2}(\mu^{*},\bar{
\mu}_{T}^{\mathcal{P_{\theta}}})\leq\inf_{\theta\in\Theta}\mathcal{W}_{2}(\mu^
{*},\bar{\mu}_{T}^{\mathcal{P_{\theta}}})+\epsilon\right\}$ is bounded.

###### Theorem 2.6 .

Consistency of the structural loss.
Assuming that 2.2,2.3,2.4 and 2.5 hold, the loss function in (5) $\mathcal{W}_{2}(\hat{{\bar{\mu}}}^{N,\mathcal{P_{\theta}}}_{T},\hat{\mu}^{N})
\leq\epsilon_{l}$ where $\epsilon_{l}\to 0$, $\mathbb{P}$-almost surely. Then we exist $E\subset\Omega$ with $P(E)=1$ such that for all $\omega\in E$ :

$$ $\inf\limits_{\mathcal{P}_{\theta}\in\Theta}\mathcal{W}_{2}(\hat{\mu}^{N}( \omega),\bar{\mu}_{T}^{\mathcal{P_{\theta}}})\to\inf\limits_{\mathcal{P}_{ \theta}\in\Theta}\mathcal{W}_{2}(\mu^{*},\bar{\mu}_{T}^{\mathcal{P_{\theta}}}),$ (8) $$

and there exists $n(\omega)$ such that,for all $N\geq n(\omega)$

$$ $\mathcal{P}_{\theta}^{N}\to\mathcal{P}^{*}_{\theta}\text{ as }N\to\infty, \epsilon_{l}\to 0,\mathbb{P}\text{-almost surely}.$ (9) $$

###### Proof.

The proof is based on []. However, the key difference is that we introduced a loss function control term, which enhances the result. We discuss details in the appendix and give a proof.
∎

After introducing the theorem above, we provide an error estimation for the diffusion bridge model.

###### Theorem 2.7 .

Consistency of the diffusion bridge model: Assuming that 2.2,2.3,2.4 and 2.5 hold, the loss function in (5) satisfied $\mathcal{W}_{2}(\hat{{\bar{\mu}}}^{N,\mathcal{P_{\theta}}}_{T},\hat{\mu}^{N})
\leq\epsilon_{l}$ where $\epsilon_{l}\to 0$,, $\mathbb{P}$-almost surely for ${\mathcal{P}_{\theta}}$.

$$ $\mathcal{W}_{2}(\mu^{\mathcal{P}_{\theta}}_{T},\mu^{*})\to 0\text{ as }N\to \infty,\epsilon_{l}\to 0,\mathbb{P}\text{-almost surely}.$ $$

###### Proof.

This result can be proved by combining Theorem 2.1 with Theorem 2.6.
∎

Know the target distribution:
This scenario has been extensively studied using MCMC algorithms and SDE-type samplers. Our method can still be applied to match a diffusion bridge, utilizing two primary matching methods. The first method involves specifying a density function, using existing MCMC algorithms to obtain $N$ discrete points at each position $X_{t}$, and then employing the Wasserstein distance loss as described above for matching.

Note: The design of the loss function is not unique. The diffusion bridge matching method presented here serves as a baseline algorithm that can be replaced by many other algorithms. We employ a Neural SDE bridge for the following reasons:
(1) We aim to minimize the number of steps to reach the target distribution within the smallest possible time interval to reduce the amount of PINN training.
(2) In cases with only partially observed samples, where density information is absent, the matching process is not unique and relies on the chosen model.
(3) We simulate an equal number of Brownian motion paths and use a fully trainable initial value for drift and diffusion coefficients to ensure maximum flexibility. The Wasserstein distance guarantees the stability of training and overall match between the generated samples and the target value.

In practical scenarios, the maximum number of training points for the PINN is $MN$, where $M$ is the number of iterations of the SDE solver, satisfying $(M-1)h=T$. Given $\varepsilon$ and $M_{0}$, we aim to achieve $M\leq M_{0}$ by choosing appropriate $h$ and
$T$ such that the following conditions hold simultaneously:

$$ ${[\frac{T}{h}]\leq M_{0}}\quad\text{ and }\quad Ch^{\frac{1}{2}}\exp{(4L^{2}T) }\leq\varepsilon.$ (10) $$

This is relatively easy to achieve because we have parameterized the initial values, allowing us to control them and, consequently, control $C$. This approach is distinctive and innovative compared to other diffusion bridge models. Additionally, this method can also match the Markov chain in cases where the density is known. The first one involves specifying a density function and then using existing MCMC algorithms to obtain $N$ discrete points at each position $X_{t}$. We then employ the same Wasserstein distance loss (6) as mentioned above for matching. The second method involves using transition density matching in []. Specifically, given a density function $f$, we can determine a transition density function $p(y|x,h)$ in MCMC algorithms. Then, by discretizing the SDE using the Euler-Maruyama method, we obtain the following transition density:

$$ $\hat{p}(y|x,h)=\mathcal{N}(y;x+b_{\theta_{1}}(x,t)h,h\sigma_{\theta_{2}}(x,t) \sigma_{\theta_{2}}^{T}(x,t))\text{.}$ (11) $$

We can consider the following loss function:

$$ $\mathcal{P_{\theta}}^{*}=\mathop{\rm argmin}_{\mathcal{P_{\theta}}\in\Theta} \iint[\hat{p}(y|x,h)-p(y|x,h)]^{2}dydx+[X_{0}-X_{0,\theta_{3}}]^{2}\text{.}$ (12) $$

Using an SDE-type sampler directly as the diffusion bridge is also feasible, eliminating the need for a matching process. A straightforward method for this purpose is the Langevin diffusion. Our experiments demonstrate the improved estimates provided by FKEE for the Langevin diffusion equation. The parameter pairs determining the diffusion bridge are $(X_{0},b,\sigma,T)$. Many MCMC algorithms can be reduced to a diffusion bridge model, as shown in Table 1 in Appendix 8.1.

###### Theorem 2.1 .

###### Proof.

###### Assumption 2.2 .

###### Assumption 2.3 .

###### Assumption 2.4 .

###### Assumption 2.5 .

###### Theorem 2.6 .

###### Proof.

###### Theorem 2.7 .

###### Proof.

### 2.3 Feynman-Kac model

This section presents our main contribution: a novel approach to expectation estimation. We aim to estimate $\mathbb{E}_{X\sim P}[f(X)]$ by decoding $P$. All relevant information about $P$ is encapsulated in $(X_{0},b,\sigma,T)$. The decoding loss measures the accuracy of our estimate, while the decoding speed impacts the algorithm’s efficiency. Our key innovation is the direct utilization of information within $(X_{0},b,\sigma,T)$, as it contains all the necessary information about $P$. This is the core of our algorithm. The decoding process can be viewed as an approximation of the Feynman-Kac operator, formally obtained by solving the Feynman-Kac equation.

The Feynman-Kac operator [] is crucial in translating between deterministic PDEs and stochastic processes through the Feynman-Kac formula (Feynman-Kac equation). The Feynman-Kac equation [] is a powerful method for solving PDEs by linking them to stochastic processes. The basic idea is to represent the solution of a PDE as the expectation of a function of a stochastic process and use Monte Carlo methods to approximate this expectation.

In our approach, we can reverse the process to obtain new methods for deriving MCMC results. Specifically, we can use the solution of a PDE to accurately express the corresponding MCMC results. We consider the following simplified version of the Feynman-Kac formula, which is commonly encountered.

###### Theorem 2.8 .

Feynman-Kac formula: Assuming that SDE (1) has strong solutions and $f$ is a function in ${C}^{2}$, we have the following Feynman-Kac equation, which has unique solutions on the interval $[0,T]$.

$$ $\frac{\partial u(x,t)}{\partial t}+\sum_{i=1}^{d}\frac{\partial u(x,t)}{ \partial x_{i}}b_{i}(x,t)+\frac{1}{2}\sum_{i=1}^{d}\sum_{j=1}^{d}\frac{ \partial^{2}u(x,t)}{\partial x_{i}\partial x_{j}}(\sigma(x,t)\sigma(x,t)^{T})_ {i,j}=0,$ $$

with the boundary conditions

$$ $u(x,T)=f(x).$ $$

The solution to the Feynman-Kac equation at the initial time is $u(x_{0},0)=\mathbb{E}[f(X_{T})|X_{0}=x_{0}]$.

###### Proof.

The proof of this theorem involves using the Itô formula and properties of martingales with stochastic integrals. A brief explanation is provided in the appendix. For more details, please refer to [].
∎

Fast calculation method. Calculating this equation involves computing the Hessian matrix of a function and some partial derivatives, which can be obtained using any library with automatic differentiation, such as Pytorch. If we consider only the diagonal diffusion coefficients $\sigma$, the algorithm can be accelerated. For instance, in Langevin diffusion where $\sigma=I_{d\times d}$, we need to calculate the second-order partial derivatives of the main diagonal.

For the Neural SDE, to reduce computation, we can consider a diagonal diffusion matrix function $\sigma$ : $\mathbb{R}^{d}\times[0,T]\to\Lambda(\mathbb{R}^{d})$, where $\Lambda(\mathbb{R}^{d})$ is the set of real-valued diagonal matrices. We only calculate the second-order derivatives of the diagonal elements to avoid the entire diffusion matrix function. To achieve this, we design the following loss functions:

$$ $\mathcal{L}_{1}=\iint_{\mathcal{D}\times[0,T]}\left[\frac{\partial u_{\theta}( x,t)}{\partial t}+\sum_{i=1}^{d}\frac{\partial u_{\theta}(x,t)}{\partial x_{i} }b_{i}(x_{t},t)\right.\left.+\frac{1}{2}\sum_{i=1}^{d}\frac{\partial^{2}u_{ \theta}(x,t)}{\partial x^{2}_{i}}\text{diag}(\sigma^{2}(x,t))_{i}\right]^{2}\, dx\,dt$ $$

and

$$ $\mathcal{L}_{2}=\int_{\mathcal{D}}\left[u_{\theta}(x,T)-f(x)\right]^{2}dx.$ $$

Finally, we obtain the solution by optimizing these two loss functions.

$$ $u^{*}(x,t)=\mathop{\rm argmin}_{u_{\theta}(x,t)}\left[\lambda_{1}\mathcal{L}_{ 1}+\lambda_{2}\mathcal{L}_{2}\right],$ $$

where $\lambda_{1}$ and $\lambda_{2}$ are the weights of the two loss functions. $u_{\theta}(x,t)$ is the neural network with a tanh activation function. Ultimately, we can obtain the expectation $u^{*}(x_{0},0)=\mathbb{E}[f(X_{T})|X_{0}=x_{0}]$. Detailed implementation of the algorithms is provided in the Appendix.

Figure: Algorithm 2 Feynman-Kac model (FCM)

The term $\mathcal{L}_{d}=\lambda_{1}\mathcal{L}_{1}+\lambda_{2}\mathcal{L}_{2}$ represents the empirical decoding loss, which incurs statistical and optimization errors compared to the true decoding loss. Error analysis for this equation can be found in many works related to PINN, for example, in [].

Viewing MCMC expectation estimators from a decoding perspective The decoding loss of traditional MCMC expectation estimators might be suboptimal because these estimators often do not fully utilize the information about $(X_{0},b,\sigma,T)$. Typically, these estimators rely on simulating a subset of samples for averaging, which can introduce local bias and fail to provide a comprehensive estimate of the entire distribution. The method of control functions [, ] attempts to mitigate this by reusing information from $P$, but it is not universally applicable to any $P$ and $f$. Additionally, these methods require stronger assumptions to guarantee accurate estimates, influenced by the Law of Large Numbers (LLN) and the Ergodic Theorem for Markov Chains (ETMC), which can further reduce the efficiency of MCMC algorithms. Therefore, traditional MCMC expectation estimators can be seen as incomplete decoding.

Discussion of the choice of the Feynman-Kac model Our approach fundamentally changes how expectations are calculated by utilizing the full distributional information in $P$ to approximate $(X_{0}^{*},b^{*},\sigma^{*})$. However, when approximating $(X_{0}^{*},b^{*},\sigma^{*})$ in many diffusion bridge models, it is often necessary to simulate part of the Brownian motion trajectory to estimate the loss function. This results in some positions $(x,t)$ corresponding to $(b,\sigma)$ being accurate, while others depend on the network’s generalization ability. Consequently, the appearance of $x$ in our position $(x,t)$ occurs randomly, necessitating a meshless PDE solver.

For certain $(b,\sigma)$ with exact analytical forms and diffusion bridges that exhibit better generalization, a non-meshless PDE solver may suffice. The second critical issue is the change in how expectations are computed, introducing the dimension $d$ with respect to the MCMC expectation estimator. To overcome the curse of dimensionality, we need a PDE solver capable of handling this problem. For low-dimensional, non-meshless scenarios, finite element methods [] are viable. However, in more general cases, we require meshless PDE solvers that can address the curse of dimensionality. We have chosen a classical PDE solver called PINN, but other PDE solvers meeting these conditions are also feasible.

###### Theorem 2.8 .

###### Proof.

### 2.4 Feynman-Kac Operator Expectation Estimator

FKEE consists of two parts: the Diffusion Bridge Model and the Feynman-Kac Model. The Diffusion Bridge Model provides the coefficients and initial values of the SDE for the Feynman-Kac Model. For a target distribution, we first use the Diffusion Bridge Model to obtain the corresponding coefficients and save them. The Feynman-Kac Model then uses these coefficients to directly approximate $\mathbb{E}_{X\sim P}[f(X)]$. Since the Feynman-Kac Model is trained using PINN, we can leverage GPU arithmetic acceleration or parallelism to efficiently handle high-dimensional distributions and obtain the corresponding results.

## 3 Discussion

Discussions on $P$: In conventional MCMC algorithms like Langevin diffusion, extensive analysis is conducted on the properties of potential energy functions, particularly the requirements for Lipschitz continuous gradients and strong convexity, as detailed in [, ]. However, our approach diverges by not depending on these specific properties of energy functions for convergence and speed. Instead, we require the corresponding SDE to have strong solutions and the Feynman-Kac equation to be well defined. To better understand the applicability of our method, consider the Itô-type SDE (1), which corresponds to the Fokker–Planck–Kolmogorov (FPK) equation [, ]:

$$ $\frac{\partial p(x,t)}{\partial t}=-\sum_{i}\frac{\partial}{\partial x_{i}}[b_ {i}(x,t)p(x,t)]+\frac{1}{2}\sum_{i,j}\frac{\partial^{2}}{\partial x_{i} \partial x_{j}}\left\{\left[\sigma(x,t)\sigma^{\top}(x,t)\right]_{ij}p(x,t) \right\},$ $$

where $p(x,t)$ is the probability density function of $X_{t}$. For the stationary distribution, we set $\frac{\partial p(x,t)}{\partial t}=0$. There can be multiple pairs $(b,\sigma)$ that satisfy this stationary FPK equation, with Langevin diffusion being a special case where $\sigma=I_{d}$. Our method can handle various other cases as well, such as those discussed in []. Additionally, our method applies to more general pairs $(b,\sigma)$ that satisfy the FPK equation, even under finite time and non-stationary conditions.

Discussion on $f$:
In classical MCMC expectation estimators, the following computation forms are used:

$$ $\mathbb{E}\left[f(X)\right]=\frac{1}{N}\sum_{i=1}^{N}f(X_{T}^{i})\text{,}$ (13) $$

where $X_{T}^{i}$ is the value at moment $T$ of the $i$th Markov chain, with different Markov chains being independent. This represents the classical Monte Carlo integral calculation, where error is based on the Law of Large Numbers (LLN). For simple distributions, such as the uniform distribution in $\left[0,1\right]$, bias occurs when $f(x)=x^{n}$ for large $n$.

Another estimator is applicable only when $P$ is the stationary distribution:

$$ $\mathbb{E}\left[f(X)\right]=\frac{1}{N-M}\sum_{t=M}^{N}f(X_{t})\text{.}$ (14) $$

Here, averaging is done over the time span of a Markov chain, with $M$ denoting the number of samples discarded during the burn-in period, characterized by correlated samples. Error in this case is influenced by the Ergodic Theorem for Markov Chains (ETMC), making optimal $M$ selection challenging for complex problems. Incorporating relevant samples can reduce the impact on MCMC expectation estimator efficiency. In difficult scenarios, properties of $f$ can often lead to larger biases. Our approach enhances MCMC efficiency by utilizing points within the burn-in period for PDE loss computation and refining assumptions on $f$, offering a novel bias reduction method.

In summary, one method relies on the LLN, often requiring Lipschitz continuity of $f$, with the estimator’s variance related to the Lipschitz coefficient of $f$. The other is based on the ETMC, also imposing requirements on the Lipschitz coefficient and the density function of $P$. Our method, however, only requires that the boundary conditions of the PDE satisfy a specific smoothness, namely $f\in C^{2}$. This significantly broadens the scope of this approach.

## 4 Experiments

Partition Function Computation for Random Graph Models.
In our first example, we focus on computing the partition function for random graph models, simplifying the setup from [] to estimate the mathematical expectation and the corresponding partition function. Background and details can be found in Appendix 9.1. Our method shows remarkable efficiency in handling distributions on random graphs, particularly with complexities introduced by the target expectation index and the function $H(x)$. For example, in a high-dimensional scenario with $n=15$, involving graphical models with partition functions summing over $2^{225}$ discrete points, our method proves effective. This improvement is evident when comparing our method with the sample points $w_{i}$, $v_{i}$ and the MCMC estimator from []. Further validation can be found at [https://github.com/zysophia/Doubly_Adaptive_MCMC](https://github.com/zysophia/Doubly_Adaptive_MCMC). We achieve comparable accuracy with only 2000 points in the Markov chain, while reducing computation time compared to current state-of-the-art MCMC expectation estimators. The results of this experiment are presented in Appendix 9.1. This experiment highlights the effect of $f$ on the target distribution’s expectation and the algorithm’s efficiency, defined here as using fewer points on the Markov chain to achieve higher accuracy in approximating expectations.

Evaluating the Diffusion Bridge Model.
The second example evaluates the impact of our proposed diffusion bridge model, detailed in Appendix 9.2. We generate a set of samples, build a diffusion bridge model based on these samples, and use it to generate additional samples. We then compare the distributions of the initial and subsequent sample sets, showcasing the effectiveness of the diffusion bridge model. The resampled data’s impact can also be observed in the first experiment.

Baseline Experiments on Properties of $P$ and Low Variance.
Additional baseline experiments on the properties of $P$ and low variance are presented in Appendix 9.3. We simulate the same trajectory using the Langevin diffusion model but employ different methods for expectation computation, resulting in distinct outcomes visualized as distributions. We also explore high-dimensional examples to further validate our approach.

## 5 Conclusion

We have introduced a heuristic method for estimating mathematical expectations by bridging the gap between deep learning PDE solvers and sampling methods. This approach reduces reliance on traditional assumptions (LLN and ETMC) and expands the applicability to a broader range of $P$ and $f$. We presented a versatile diffusion bridge model to extend the range of $P$ and utilized PDE methods to broaden the scope of $f$. Our method demonstrates potential significance across multiple domains.

## 6 Appendix / supplemental material

## 7 Related work

The diffusion model belongs to a class of stochastic differential equations, which are used to approximate the target distribution. It has been widely used for generative models [], variational inference [, ], etc.
The diffusion bridge model is a variant of the diffusion model. Early development of diffusion bridge models involved simulating processes originating from two endpoints []. Alternative approaches for constructing diffusion bridge models are outlined by [, ]. Since the diffusion bridge model essentially functions as a sampling algorithm, it plays a pivotal role in addressing the crucial task of high-dimensional distribution sampling. Sampling high-dimensional distributions is a fundamental task with applications across various fields. Common methods include MCMC, random flow, and generative models. Recent work includes stream-based methods [, , , , ], MCMC-based methods [, , ] and generative models [], score-based models []. Normalizing Flows []. These models can be broadly categorized into two groups: those based on given discrete points and those relying on a given density function. The former primarily serves for learning and generating real world data such as text and images, while the latter is used for sampling, statistical estimation, and similar purposes. Notably, Langevin diffusion [, , ] is a classical model within the latter category.

The Feynman-Kac model is a technique employed to solve partial differential equations (PDEs) by using deep learning. Deep learning has found application in solving PDEs of the Feynman-Kac equation type, as demonstrated by [, ]. [] highlights the approximation of Feynman-Kac type expectations through the approximation of discrete Markov chains, thereby enhancing the order of convergence. When employing PINN to solve Feynman-Kac type PDEs, the sampling algorithm can be linked to the path of the SDE. This approach enables the acquisition of adaptive sampling points from the paths of SDE, which proves more efficient than uniform point selection []. Further analysis of the approximation error for this class of equations is presented in [].

## 8 Definitions and related theory

### 8.1 Comparison of Diffusion bridge model.

**Table 1: Comparison of Diffusion bridge model**
| Method | $X_{0}$ | $b$ | $\sigma$ | $T$ | Descriptions |
| --- | --- | --- | --- | --- | --- |
| Classical MCMC | $\forall\mathbf{x}_{0}\in\mathbb{R}^{d}$ | $p(y|x)$ | $p(y|x)$ | $\infty$ | $p(y|x)$ is the transfer probability density function in the MCMC algorithm. The meaning of $p(y|x)$ is that the corresponding coefficients can be obtained by a SDE. |
| Langevin MCMC | $\forall\mathbf{x}_{0}\in\mathbb{R}^{d}$ or<br>$\forall X_{0}\sim P_{0}$ | $\frac{1}{2}\nabla_{x}\log p(x)$ | $I_{d\times d}$ | $\infty$ | $p(x)$ is target density function and density function of a stationary distribution. $P_{0}$ is the initial distribution. |
| Score-based SDE and diffusion models (DDPM) | $\forall X_{0}\sim P_{0}$ | $f(\mathbf{x},t)-g^{2}(t)$<br>$\nabla_{\mathbf{x}}\log p_{t}(\mathbf{x})$ | $g(t)$ | $\infty$ | $\nabla_{\mathbf{x}}\log p_{t}(\mathbf{x})$ is obtained from the data and $f(x,t)$ and $g(t)$ are known. $P_{0}$ is the prior distribution. |
| Flow match ODE | $\forall X_{0}\sim P_{0}$ | $v(x,t)$ | 0 | 1 | $v(x,t)$ is obtained by matching the data. $P_{0}$ is the initial distribution. |
| Neural SDE bridge (taken in this paper) | $x_{0}=x_{0,\theta{{}_{3}}}$ | $b_{\theta{{}_{1}}}(x,t)$ | $\sigma_{\theta{{}_{2}}}(x,t)$ | $<\infty$ | $b_{\theta{{}_{1}}}(x,t)$ and $\sigma_{\theta{{}_{2}}}(x,t)$ is obtained from the data or match method. |

### 8.2 Definitions

Wasserstein distance: The most commonly used measure of distance between probability distributions is the Wasserstein distance. It calculates the minimum cost of transporting mass from one distribution to another, based on the distance between the points being transported and the amount of mass being moved. The Wasserstein distance is especially beneficial for comparing distributions with different shapes since it considers the structure of distributions instead of only their statistical moments. This distance metric is widely applied in fields like image processing, computer vision, and machine learning. The definition is

$$ $\mathcal{W}_{p}(\mu,\nu):=\left(\inf_{\pi\in\Pi(\mu,\nu)}\int_{\mathbb{R}^{d} \times\mathbb{R}^{d}}|x-y|^{p}\pi(dx,dy)\right)^{\frac{1}{p}}=\inf\left\{\left [\mathbb{E}|X-Y|^{p}\right]^{\frac{1}{p}},\mathbb{P}_{x}=\mu,\mathbb{P}_{Y}= \nu\right\}.$ $$

$\Pi(\mu,\nu)$ denotes the class of measures on $\mathbb{R}^{d}\times\mathbb{R}^{d}$ with marginal distributions $\mu$ and $\nu$.

Euler-Maruyama method: [] is a frequently used approach for solving SDE through an iterative format. This method has been shown to converge to a strong order of $\mathcal{O}(h^{\frac{1}{2}})$, where the error is dependent on the Lipschitz coefficients of the drift and diffusion coefficients. When generating paths using this method, it is recommended to use smaller step sizes, to minimize the errors associated with the method.

$$ $X_{t+h}=X_{t}+b(X_{t},t)h+\sigma(X_{t},t)(W_{t+h}-W_{t}),X_{0}=x_{0}.$ $$

Numerical solvers for stochastic differential equations of any accuracy are allowed when constructing sample paths for diffusion.

Physics-informed neural networks: PINN [] is a deep learning method for solving partial differential equations. The main idea is to use neural networks for fitting solutions to PDE problems, PINN incorporates the residuals of the PDE (the difference between the left-hand side and the right-hand side of the PDE equation) into the loss function, and then updates the weights and parameters of the neural network through a backpropagation algorithm. Specifically, we consider follow PDE:

$$ $F(u_{t},u_{x},u_{xx})=g(u,x,t)$ $$

and the boundary condition is

$$ $G(u_{t},u_{x},u_{xx})=0$ $$

We choose a neural network $u^{\theta}(x,t)$ to approximate the solution $u(x,t)$. By automatic differentiation, we can easily obtain the term $u^{\theta}_{t}$,$u^{\theta}_{x}$ and $u^{\theta}_{xx}$.
We then need to sample the region of the target and calculate the value of the empirical loss function for these points. Finally the solution $u^{\theta}_{t}$ is obtained by optimising the combination of the two loss functions.

$$ $\textit{Loss PDE}=F(u^{\theta}_{t},u^{\theta}_{x},u^{\theta}_{xx})-g(u^{\theta },x,t)\quad\textit{Loss boundary}=G(u^{\theta}_{t},u^{\theta}_{x},u^{\theta}_{ xx})$ $$

$$ $\textit{Loss}=\lambda_{1}\textit{Loss boundary}+\lambda_{2}\textit{Loss PDE}$ $$

$\lambda_{1}$ and $\lambda_{2}$ are the weights of the two loss functions.

### 8.3 Proof of theoretical results

###### Theorem 8.1 .

If the drift and diffusion coefficients satisfy the conditions [] in SDE (1):

- •
Lipschitz condition
$|b(x,t)-b(y,t)|\leq K|x-y|\text{ and }||\sigma(x,t)-\sigma(y,t)||_{F}\leq K|x-y|$
for all $t\in\left[{0},T\right]$ and $x,y\in\mathbb{R}^{d}$
- •
Linear growth bound
There exists a constant $C$ such that
$|b(x,t)|^{2}\leq C^{2}(1+|x|^{2})\quad\text{and}\quad||\sigma(x,t)||_{F}^{2}
\leq C^{2}(1+|x|^{2})$
for all $t\in\left[{0},T\right]$ and $x,y\in\mathbb{R}^{d}\text{. }$
- •
Measurability
$b(x,t)\text{ and }\sigma(x,t)\text{ is jointly measurable.}$
- •
Initial value
$X_{0}\text{ is }\mathcal{F}_{{0}}\text{-measurable with }\mathbb{E}(|X_{{0}}|^
{2})<\infty\text{. }$

where $||\cdot||_{F}$ denotes the Frobenius norm of a matrix.
Then the SDE has a unique strong solution.
The solution can be controlled by the initial value, i.e., $\mathbb{E}X^{2}_{T}\leq C\mathbb{E}X^{2}_{0}$.

Proof for Theorem 2.1.

###### Proof.

Base on the Lipschitz condition, Linear growth bound condition in Theorem 8.1, and the Itô isometry, we can derive the following:

$$ $\displaystyle\mathcal{W}^{2}_{2}(\mu_{T},\bar{\mu}_{T})\leq\mathbb{E}|X_{T}- \bar{X}_{T}|^{2}$ $\displaystyle\leq 4\mathbb{E}\left|\int_{0}^{T}b(X_{t},t)-b(\bar{X}_{t},t)dt \right|^{2}+4\mathbb{E}\left|\int_{0}^{T}\sigma(X_{t},t)-\sigma(\bar{X}_{t},t) dW_{t}\right|^{2}.$ $\displaystyle+4\mathbb{E}\left|\int_{[T]}^{T}\sigma(\bar{X}_{t},t)dW_{t}\right |^{2}+4\mathbb{E}\left|\int_{[T]}^{T}b(\bar{X}_{t},t)dt\right|^{2}.$ $\displaystyle\leq 4L^{2}\mathbb{E}\int_{0}^{T}|X_{t}-\bar{X}_{t}|^{2}dt+4L^{2} \mathbb{E}\int_{0}^{T}|X_{t}-\bar{X}_{t}|^{2}dt+4\mathbb{E}\left|\int_{[T]}^{T }\sigma(\bar{X}_{t},t)dW_{t}\right|^{2}+4\mathbb{E}\left|\int_{[T]}^{T}b(\bar{ X}_{t},t)dt\right|^{2}.$ $\displaystyle\leq 8L^{2}\mathbb{E}\int_{0}^{T}|X_{t}-\bar{X}_{t}|^{2}dt+4C^{2} \mathbb{E}\left|\int_{[T]}^{T}\bar{X}_{t}^{2}dt\right|+4C^{2}\mathbb{E}\left| \int_{[T]}^{T}\bar{X}_{t}dt\right|^{2}+4C^{2}h^{2}.$ $\displaystyle\leq 8L^{2}\mathbb{E}\int_{0}^{T}|X_{t}-\bar{X}_{t}|^{2}dt+4C^{2} M_{0}h+4C^{2}M_{1}h+4C^{2}h.$ $$

where $M_{0}$ and $M_{1}$ are upper bounds that relevant to $\mathbb{E}(|\bar{X}_{{T}}|^{2})$ because of $\mathbb{E}(|X_{{0}}|^{2})<\infty$. $[T]:=\max\left\{Mh\leq T\right\}$. Finally, based on the Gronwall’s inequality.

$$ $\mathcal{W}^{2}_{2}(\mu_{T},\bar{\mu}_{T})\leq\mathbb{E}|X_{T}-\bar{X}_{T}|^{2 }\leq Ch\exp{(8L^{2}T)}.$ $$

∎

Proof for Theorem 2.6.

Before proving this theorem we introduce the following definitions and lemmas. Some of these definitions and lemmas are taken from []

###### Definition 8.2 .

The function $f:\Theta\to\mathbb{R}$ : is lower semicontinuous at $x_{0}$ if

$$ ${\lim{\inf}}_{x\to x_{0}}f(x)\geq f(x_{0}).$ (15) $$

###### Definition 8.3 .

A sequence of functions $f_{n}:\Theta\to\mathbb{R}$ is said to epi-converge to $f:\Theta\to\mathbb{R}$ if for all $\theta\in\Theta$

$$ $\begin{cases}\lim\inf_{n\to\infty}f_{n}(\theta_{n})\geq f(\theta)&\text{for every sequence }\theta_{n}\to\theta,\\ \lim\sup_{n\to\infty}f_{n}(\theta_{n})\leq f(\theta)&\text{for some sequence } \theta_{n}\to\theta.\end{cases}$ (16) $$

###### Lemma 8.4 .

The sequence $f_{n}:\Theta\to\mathbb{R}$ epi-converges to $f:\Theta\to\mathbb{R}$ if and only if

$$ $\begin{cases}\liminf_{n\to\infty}\inf_{\theta\in\mathcal{K}}f_{n}(\theta)\geq \inf_{\theta\in\mathcal{K}}f(\theta)&\textit{for every compact set }\mathcal{K }\subset\Theta,\\ \limsup_{n\to\infty}\inf_{\theta\in\mathcal{O}}f_{n}(\theta)\leq\inf_{\theta \in\mathcal{O}}f(\theta)&\textit{for every open set }\mathcal{O}\subset\Theta. \end{cases}$ (17) $$

###### Lemma 8.5 .

Varadarajan’s theorem: If $X_{1},\ldots,X_{n}$ are i.i.d. $X\sim P$ on a separable metric space, then $\hat{\mu}^{n}\stackrel{{\scriptstyle w}}{{\to}}\mu^{*}$ $\mathbb{P}$-almost surely where $\hat{\mu}^{n}$ is empirical measure.

###### Lemma 8.6 .

Attainment of a minimum: Suppose $f:\Theta\to\mathbb{R}$ is lower semicontinues, level-bounded and proper. Then the value $\inf f$ is finite and the
set $\mathop{\rm argmin}f$ is nonempty and compact.

###### Lemma 8.7 .

The properties of epi-convergence:
If $f_{1}^{n}\leq f^{n}\leq f_{2}^{n}$ with $f^{n}_{1}\overset{\mathrm{epi}}{\to}f$ and $f^{n}_{2}\overset{\mathrm{epi}}{\to}f$, then $f^{n}\overset{\mathrm{epi}}{\to}f$.

###### Lemma 8.8 .

The limits of inf: Suppose $f_{n}\overset{\mathrm{epi}}{\to}f$ with $-\infty<\inf f<\infty$. Then $\inf f^{n}\to\inf f$ if and only if there exists for every $\varepsilon>0$ a compact set $B\subset R^{n}$ along with an index set $\mathcal{N}$ such that

$$ $\inf_{B}f^{n}~{}\leq~{}\inf f^{n}+\varepsilon\quad\textit{ for all }n\in \mathcal{N}.$ $$

###### Proof.

Based on Assumptions 2.4 and []. we can conclude the map is lower semicontinuous. i,e, $\mathcal{W}_{2}(\bar{\mu}_{T}^{\mathcal{P}_{\theta}},\nu)\leq{\lim\inf}_{N\to
\infty}\mathcal{W}_{2}(\bar{\mu}_{T}^{\mathcal{P}^{N}_{\theta}},\nu)$.

Then, we aim to prove $\mathcal{P_{\theta}}\mapsto\mathcal{W}_{2}(\hat{\mu}^{N},\bar{\mu}_{T}^{
\mathcal{P}_{\theta}})$ epi converges to $\mathcal{P_{\theta}}\mapsto\mathcal{W}_{2}(\mu^{*},\bar{\mu}_{T}^{\mathcal{P}_
{\theta}})$ $\mathbb{P}$-almost surely.

Firstly, we can observe the following inequality:

$$ $\mathcal{W}_{2}(\mu^{*},\bar{\mu}_{T}^{\mathcal{P}_{\theta}})-\mathcal{W}_{2}( \mu^{*},\hat{\mu}_{T}^{N})\leq\mathcal{W}_{2}(\hat{{\mu}}_{T}^{N},\bar{\mu}_{T }^{\mathcal{P}_{\theta}}),$ (18) $$

and

$$ $\mathcal{W}_{2}(\hat{{\mu}}_{T}^{N},\bar{\mu}_{T}^{\mathcal{P}_{\theta}})\leq \mathcal{W}_{2}(\hat{{\mu}}_{T}^{N},\hat{\bar{\mu}}_{T}^{N,\mathcal{P}_{\theta }})+\mathcal{W}_{2}(\hat{\bar{\mu}}_{T}^{N,\mathcal{P}_{\theta}},\bar{\mu}_{T} ^{\mathcal{P}_{\theta}}).$ (19) $$

In inequality (18), the function $\mathcal{P}_{\theta}\mapsto\mathcal{W}_{2}(\mu^{*},\bar{\mu}_{T}^{\mathcal{P}_
{\theta}})$ epi-converges to $\mathcal{P}_{\theta}\mapsto\mathcal{W}_{2}(\mu^{*},\bar{\mu}_{T}^{\mathcal{P}_
{\theta}})$ because it is independent of $N$, and the function $\mathcal{P}_{\theta}\mapsto\mathcal{W}_{2}(\mu^{*},\hat{\mu}_{T}^{N})$ epi-converges to 0 0 as $N\to\infty$, due to Assumption 2.3.

In inequality (19), the function $\mathcal{P}_{\theta}\mapsto\mathcal{W}_{2}(\hat{{\mu}}_{T}^{N},\hat{\bar{\mu}}
_{T}^{N,\mathcal{P}_{\theta}})$ epi converge to 0, because of assumptions about the loss function. Finally we aim to proof that the function $\mathcal{P}_{\theta}\mapsto\mathcal{W}_{2}(\hat{\bar{\mu}}_{T}^{N,\mathcal{P}_
{\theta}},\bar{\mu}_{T}^{\mathcal{P}_{\theta}})$ epi converge to $\mathcal{P}_{\theta}\mapsto\mathcal{W}_{2}(\mu^{\mathcal{P^{*}_{\theta}}},\bar
{\mu}_{T}^{\mathcal{P}_{\theta}})$. Combining Assumption 2.2 $\mathcal{W}_{2}(\mu^{\mathcal{P^{*}_{\theta}}},\nu)=\mathcal{W}_{2}(\mu^{*},\nu)$ and Lemma 8.7 we can conclude that $\mathcal{P_{\theta}}\mapsto\mathcal{W}_{2}(\hat{\mu}^{N},\bar{\mu}_{T}^{
\mathcal{P}_{\theta}})$ epi converges to $\mathcal{P_{\theta}}\mapsto\mathcal{W}_{2}(\mu^{*},\bar{\mu}_{T}^{\mathcal{P}_
{\theta}})$ $\mathbb{P}$-almost surely.

In the following we demonstrate that the function $\mathcal{P}_{\theta}\mapsto\mathcal{W}_{2}(\hat{\bar{\mu}}_{T}^{N,\mathcal{P}_
{\theta}},\bar{\mu}_{T}^{\mathcal{P}_{\theta}})$ epi converge to $\mathcal{P}_{\theta}\mapsto\mathcal{W}_{2}(\mu^{\mathcal{P^{*}_{\theta}}},\bar
{\mu}_{T}^{\mathcal{P}_{\theta}})$.

According to Lemma 8.4, we go on to verify two inequalities. For a compact set $\mathcal{K}$, by the lower semicontinuous of the map $\mathcal{P}_{\theta}\mapsto\mathcal{W}_{2}(\nu,\bar{\mu}_{T}^{\mathcal{P}_{
\theta}})$, we have

$$ $\inf_{\mathcal{P}_{\theta}\in\mathcal{K}}\mathcal{W}_{2}(\hat{\bar{\mu}}_{T}^{ N,\mathcal{P}^{N}_{\theta}},\bar{\mu}_{T}^{\mathcal{P}_{\theta}})=\mathcal{W}_ {2}(\hat{\bar{\mu}}_{T}^{N,\mathcal{P}^{N}_{\theta}},\bar{\mu}_{T}^{\mathcal{P }^{N}_{\theta}}),\textit{ for some }\mathcal{P}^{N}_{\theta}\in\mathcal{K}.$ (20) $$

$$ $\displaystyle\liminf_{N\to\infty}\inf_{\mathcal{P}_{\theta}\in\mathcal{K}} \mathcal{W}_{2}(\hat{\bar{\mu}}_{T}^{N,\mathcal{P}_{\theta}},\bar{\mu}_{T}^{ \mathcal{P}^{N}_{\theta}})$ $\displaystyle\overset{\text{1}}{=}\liminf_{N\to\infty}\mathcal{W}_{2}(\hat{ \bar{\mu}}_{T}^{N,\mathcal{P}_{\theta}},\bar{\mu}_{T}^{\mathcal{P}^{N}_{\theta }})$ (21) $\displaystyle\overset{\text{2}}{=}\lim_{k\to\infty}\mathcal{W}_{2}(\hat{\bar{ \mu}}_{T}^{N_{k},\mathcal{P}_{\theta}},\bar{\mu}_{T}^{\mathcal{P}^{N_{k}}_{ \theta}})$ $\displaystyle\overset{\text{3}}{=}\lim_{m\to\infty}\mathcal{W}_{2}(\hat{\bar{ \mu}}_{T}^{N_{k_{m}},\mathcal{P}_{\theta}},\bar{\mu}_{T}^{\mathcal{P}^{N_{k_{m }}}_{\theta}}),$ $$

$$ $\displaystyle\lim_{m\to\infty}\mathcal{W}_{2}(\hat{\bar{\mu}}_{T}^{N_{k_{m}}, \mathcal{P}_{\theta}},\bar{\mu}_{T}^{\mathcal{P}^{N_{k_{m}}}_{\theta}})$ $\displaystyle=\liminf_{m\to\infty}\mathcal{W}_{2}(\hat{\bar{\mu}}_{T}^{N_{k_{m }},\mathcal{P}_{\theta}},\bar{\mu}_{T}^{\mathcal{P}^{N_{k_{m}}}_{\theta}})$ (22) $\displaystyle\overset{\text{4}}{\geq}\mathcal{W}_{2}(\mu^{\mathcal{P^{*}_{ \theta}}},\bar{\mu}_{T}^{\mathcal{P}_{\theta_{0}}})$ $\displaystyle\geq\inf_{\mathcal{P}_{\theta}\in\mathcal{K}}\mathcal{W}_{2}(\mu^ {\mathcal{P^{*}_{\theta}}},\bar{\mu}_{T}^{\mathcal{P}_{\theta}}),$ $$

where the "1" holds due to the substitution for the equation 20, "2" holds due to the definition of the infimum, "3" holds due to the compactness of the $\mathcal{K}$, "4" holds due to the Lemma 8.5 and lower semicontinuous.

For an open set $\mathcal{O}$,

$$ $\inf_{\mathcal{P}_{\theta}\in\mathcal{K}}\mathcal{W}_{2}(\hat{\bar{\mu}}_{T}^{ N,\mathcal{P}_{\theta}},\bar{\mu}_{T}^{\mathcal{P}_{\theta}})\leq\mathcal{W}_{ 2}(\hat{\bar{\mu}}_{T}^{N,\mathcal{P}_{\theta}},\bar{\mu}_{T}^{\mathcal{P}^{N} _{\theta}}),\textit{ exist }\mathcal{P}^{N}_{\theta}\in\mathcal{O}.$ (23) $$

$$ $\displaystyle\limsup_{N\to\infty}\inf_{\mathcal{P}_{\theta}\in\mathcal{O}} \mathcal{W}_{2}(\hat{\bar{\mu}}_{T}^{N,\mathcal{P}_{\theta}},\bar{\mu}_{T}^{ \mathcal{P}_{\theta}})$ $\displaystyle\overset{\text{1}}{\leq}\limsup_{N\to\infty}\mathcal{W}_{2}(\hat{ \bar{\mu}}_{T}^{N,\mathcal{P}_{\theta}},\bar{\mu}_{T}^{\mathcal{P}^{N}_{\theta }})$ (24) $\displaystyle\overset{\text{2}}{\leq}\limsup_{N\to\infty}\mathcal{W}_{2}(\hat{ \bar{\mu}}_{T}^{N,\mathcal{P}_{\theta}},\mu^{\mathcal{P^{*}_{\theta}}})+ \limsup_{N\to\infty}\mathcal{W}_{2}(\mu^{\mathcal{P^{*}_{\theta}}},\bar{\mu}_{ T}^{\mathcal{P}^{N}_{\theta}}),$ $$

where "1" holds due to a substitution for equation 23, "2" holds due to triangle inequality.

$$ $\displaystyle\limsup_{N\to\infty}\mathcal{W}_{2}(\hat{\bar{\mu}}_{T}^{N, \mathcal{P}_{\theta}},\mu^{\mathcal{P^{*}_{\theta}}})+\limsup_{N\to\infty} \mathcal{W}_{2}(\mu^{\mathcal{P^{*}_{\theta}}},\bar{\mu}_{T}^{\mathcal{P}^{N}_ {\theta}})$ $\displaystyle=\limsup_{N\to\infty}\mathcal{W}_{2}(\mu^{*},\bar{\mu}_{T}^{ \mathcal{P}^{N}_{\theta}})$ (25) $\displaystyle=\inf_{\mathcal{P}_{\theta}\in\mathcal{O}}\mathcal{W}_{2}(\mu^{*} ,\bar{\mu}_{T}^{\mathcal{P}_{\theta}}),$ $$

Hence, $\mathcal{P_{\theta}}\mapsto\mathcal{W}_{2}(\hat{\mu}^{N},\bar{\mu}_{T}^{
\mathcal{P}_{\theta}})$ epi converges to $\mathcal{P_{\theta}}\mapsto\mathcal{W}_{2}(\mu^{*},\bar{\mu}_{T}^{\mathcal{P}_
{\theta}})$ $\mathbb{P}$-almost surely holds.

By the definition of epi convergence, Theorem 8.6, and Assumption 2.3.
According to Assumption 2.2, we can find an $E$ satisfying that for a $\omega\in E$, where $\mathbb{P}(E)=1$
Based on the above process, we can conclude that

$$ $\limsup_{N\to\infty}\inf_{\mathcal{P}_{\theta}\in\mathcal{P}_{\Theta}}\mathcal {W}_{2}(\hat{\mu}^{N}(\omega),\bar{\mu}_{T}^{\mathcal{P}_{\theta}})\leq\inf_{ \mathcal{P}_{\theta}\in\mathcal{P}_{\Theta}}\mathcal{W}_{2}(\mu^{*},\bar{\mu}_ {T}^{\mathcal{P}_{\theta}})=\mathcal{P}_{\theta}^{*}.$ (26) $$

By conditioning on the loss function and the definition of the infimum, there exists $N_{1}(\omega)$ such that

for $\epsilon>0$,

$$ $\inf_{\mathcal{P}_{\theta}\in\mathcal{P}_{\Theta}}\mathcal{W}_{2}(\hat{\mu}^{N }(\omega),\bar{\mu}_{T}^{\mathcal{P}_{\theta}})\leq\inf_{\mathcal{P}_{\theta} \in\mathcal{P}_{\Theta}}\mathcal{W}_{2}(\hat{\mu}^{n}(\omega),\bar{\mu}_{T}^{ \mathcal{P}_{\theta}})+\mathcal{W}_{2}(\mu^{n}(\omega),\hat{\mu}^{N,{\mathcal{ P}_{\theta}}})\leq\mathcal{P}_{\theta}^{*}+\epsilon/2$ (27) $$

and the set

$$ $\{\mathcal{P}_{\theta}\in\mathcal{P}_{\Theta}:\mathcal{W}_{2}(\hat{\mu}^{N}( \omega),\bar{\mu}_{T}^{\mathcal{P}_{\theta}})\leq\mathcal{P}_{\theta}^{*}+ \epsilon/2\},$ (28) $$

is non-empty for all $n\geq N_{1}(\omega)$. By the triangle inequality,

$$ $\mathcal{W}_{2}(\mu^{*},\bar{\mu}_{T}^{\mathcal{P}_{\theta}})\leq\mathcal{W}_{ 2}(\mu^{*},\hat{\mu}^{N}(\omega))+\mathcal{W}_{2}(\hat{\mu}^{N}(\omega),\bar{ \mu}_{T}^{\mathcal{P}_{\theta}}),$ (29) $$

and 2.3, there exists $N_{2}(\omega)$ such that

$$ $\mathcal{W}_{2}(\mu^{*},\hat{\mu}^{N}(\omega))\leq\epsilon/2,$ (30) $$

Then

$$ $\mathcal{W}_{2}(\mu^{*},\bar{\mu}_{T}^{\mathcal{P}_{\theta}})\leq\epsilon+ \mathcal{P}_{\theta}^{*}.$ (31) $$

This means that:

$$ $\{\mathcal{P}_{\theta}\in\mathcal{P}_{\Theta}:\mathcal{W}_{2}(\hat{\mu}^{N}( \omega),\bar{\mu}_{T}^{\mathcal{P}_{\theta}})\leq\mathcal{P}_{\theta}^{*}+ \epsilon/2\}\subset B(\epsilon).$ (32) $$

Let $N_{0}=\max{(N_{1}(\omega),N_{2}(\omega))}$ when $N\geq N_{0}$ we have

$$ $\inf_{\mathcal{P}_{\theta}\in\mathcal{P}_{\Theta}}\mathcal{W}_{2}(\hat{\mu}^{N }(\omega),\bar{\mu}_{T}^{\mathcal{P}_{\theta}})=\inf_{\mathcal{P}_{\theta}\in B (\epsilon)}\mathcal{W}_{2}(\hat{\mu}^{N}(\omega),\bar{\mu}_{T}^{\mathcal{P}_{ \theta}}).$ (33) $$

According to 8.8, we can get the result:

$$ $\inf\limits_{\mathcal{P}_{\theta}\in\mathcal{P}_{\Theta}}\mathcal{W}_{2}(\hat{ \mu}^{N}(\omega),\bar{\mu}_{T}^{\mathcal{P}_{\theta}})\to\inf\limits_{\mathcal {P}_{\theta}\in\mathcal{P}_{\Theta}}\mathcal{W}_{2}(\mu^{*},\bar{\mu}_{T}^{ \mathcal{P}_{\theta}}),$ $$

and $N_{0}(\omega)=N$ the sets ${\mathop{\rm argmin}}_{\mathcal{P}_{\theta}\in\mathcal{P}_{\Theta}}\mathcal{W}
_{2}(\hat{\mu}^{N}(\omega),\bar{\mu}_{T}^{\mathcal{P}_{\theta}})$ are non-empty form a bounded sequence with

$$ $\operatorname*{lim}_{N\to\infty}\operatorname{sup}{\mathop{\rm argmin}}_{ \mathcal{P}_{\theta}\in\mathcal{P}_{\Theta}}{\mathcal{W}}_{2}(\hat{\mu}^{N}( \omega),\bar{\mu}_{T}^{\mathcal{P}_{\theta}})\subset\mathop{\rm argmin}{ \mathcal{W}}_{2}(\mu^{*},\bar{\mu}_{T}^{\mathcal{P}_{\theta}}),$ $$

further we have:

$$ $\mathcal{P}_{\theta}^{N}\to\mathcal{P}_{\theta}^{*}\text{ as }N\to\infty, \epsilon_{l}\to 0,\mathbb{P}\text{-almost surely}.$ (34) $$

∎

Next we give the proof of Theorem 2.8.

###### Theorem 8.9 .

If the SDE has strong solution 8.1,
then the solution to the corresponding backward partial differential equation (PDE) can represent the expectation of the terminal distribution of the SDE.

$$ $\frac{\partial u(x,t)}{\partial t}+\sum_{i=1}^{d}\frac{\partial u(x,t)}{ \partial x_{i}}b_{i}(x,t)+\frac{1}{2}\sum_{i=1}^{d}\sum_{j=1}^{d}\frac{ \partial^{2}u(x,t)}{\partial x_{i}\partial x_{j}}(\sigma(x,t)\sigma(x,t)^{T})_ {i,j}=0\text{,}$ $$

$$ $u(x,T)=f(x).$ $$

The soluton of PDE at $(x_{0},0)$ is $\mathbb{E}[f(X_{T})|X_{0}=x_{0}]$,i.e. $u(x_{0},0)=\mathbb{E}[f(X_{T})|X_{0}=x_{0}]$.

###### Proof.

According to Itô’s formula

$$ $du(X_{t},t)=\left[\frac{\partial u(X_{t},t)}{\partial t}+\sum_{i=1}^{d}\frac{ \partial u(X_{t},t)}{\partial x_{i}}b_{i}(X_{t},t)+\frac{1}{2}\sum_{i=1}^{d} \sum_{j=1}^{d}\frac{\partial^{2}u(X_{t},t)}{\partial x_{i}\partial x_{j}}( \sigma(X_{t},t)\sigma(X_{t},t)^{T})_{i,j}\right]dt$ $$

$$ $+\left[\sum_{r=1}^{d}\sum_{i=1}^{d}\frac{\partial u(X_{t},t)}{\partial x_{i}} \sigma_{i,r}(X_{t},t)\right]dW_{t}^{r}\text{,}$ $$

where $W_{t}^{r}$ is the rth component of $W_{t}$. The first part is based on the equality inside the PDE being set to zero. Integrate from $[0,T]$ on both sides.

$$ $u(X_{T},T)-u(X_{0},0)=f(X_{T})-u(X_{0},0)=\int_{0}^{T}\left[\sum_{r=1}^{d}\sum _{i=1}^{d}\frac{\partial u(X_{t},t)}{\partial x_{i}}\sigma_{i,r}(X_{t},t) \right]dW_{t}^{r}\text{,}$ $$

Taking the conditional expectation on both sides while fixing $X_{0}$ and utilizing the properties of Itô integration as a martingale, we get

$$ $u(x_{0},0)=\mathbb{E}[f(X_{T})|X_{0}=x_{0}].$ $$

∎

###### Theorem 8.1 .

###### Proof.

###### Definition 8.2 .

###### Definition 8.3 .

###### Lemma 8.4 .

###### Lemma 8.5 .

###### Lemma 8.6 .

###### Lemma 8.7 .

###### Lemma 8.8 .

###### Proof.

###### Theorem 8.9 .

###### Proof.

## 9 Backgrounds and Tables of the experiment

### 9.1 Experiments on the Ising model

Traditional MCMC algorithms face limitations when dealing with discrete random variables and complex functions $f$, which results in high variance. Consequently, accurate estimates often require a large number of points in the Markov chain, especially in larger models. This issue is particularly prevalent in random graph models [, , ].

Ising model: Assume a sample space $\Omega$, Hamiltonian function $H:\Omega\rightarrow\{0\}\cup[1,\infty)$ , and inverse temperature parameter $\beta\in\mathbb{R}$ , referred to as inverse temperature. The Gibbs distribution on $\Omega,H(\cdot)$ , and $\beta$ is then characterized by probability law
$\forall x\in\Omega:\pi_{\beta}(x)\doteq\frac{1}{Z(\beta)}\exp(-\beta H(x))$
Here $Z(\beta)$ is the normalizing constant or Gibbs partition function (GPF) of the distribution, with
$Z(\beta)\doteq\sum_{x\in\Omega}\exp(-\beta H(x))$. Specifically, we considered Ising model on 2D lattices: It has $n\times n$ dimensions and a total of $n^{2}$ random variables, each of which takes the values +1,-1 with Hamiltonian function $H(x)=-\sum_{(i,j)\in E}\mathbb{I}(x(i)=x(j))$. For $\beta_{0}$ the results are easy to compute and for $\left[\beta_{1},\beta_{2}\right]$ between we can use the PPE-method, we do not use the Tpa-Method [], which is an algorithm on splitting the region $\left[\beta_{1},\beta_{2}\right]$. Specifically we can compute the following
$\mathbb{E}F=\mathbb{E}{\exp(-\frac{\beta_{2}-\beta_{1}}{2}H(X_{\beta_{1}}))}$ and $\mathbb{E}G=\mathbb{E}{\exp(\frac{\beta_{2}-\beta_{1}}{2}H(X_{\beta_{2}}))}$.$X_{\beta_{i}}$ is a Gibbs distribution obeying parameter $\beta_{i}$.
$Q=\frac{\mathbb{E}G}{\mathbb{E}F}=\frac{Z(\beta_{1})}{Z(\beta_{2})}$.We set $\beta_{1}=-0.02$ and $\beta_{2}=0$. Then we can find $Z(\beta_{1})$ based on the fact that $Z(\beta_{2})=Z(0)$. So we need to estimate two mathematical expectations and we propose two ways to approximate this expectation. In this Experiment, for a definite temperature $\beta$, the distribution on the random graph is often easy to approximate, but the complexity of the exponent in the target expectation and also the function $H(x)$ can lead to the need for a large sample size to reduce the variance when MCMC deals with this problem. Our approach demonstrates superior efficiency in dealing with the distribution on a random graph, particularly when considering the complexities introduced by the target expectation exponent and the function $H(x)$. The implementation of our method, FKEE, stands out in handling larger-sized graphs ($n\geq 6$) where traditional MCMC and its variants, as found in [], face challenges due to sample complexity. We have the following two methods:

First method: direct approximation of the overall part of the expectation. That is, we consider the approximate stochastic process $H_{\beta}(X)$, which is a one-dimensional problem. We generate the chain using the same method as in [] and compute the value $H_{\beta}(X_{t})$ under each moment. The diffusion bridge model and the Feynman-Kac model are then used to estimate the expectation separately. In the diffusion bridge model, we generated the same number of Brownian motions at the same number of moments and then calculated the loss at each moment to train. The Feynman-Kac model uses the already established diffusion bridge model to get an estimate of the expectation by solving the PDE.

The second approach better exemplifies the substantial improvement in harnessing Markov chains facilitated by our method. It highlights the remarkable flexibility embedded in our approach. Specifically, we directly approximate the distribution on a random graph, conceptualizing this graph as an $n^{2}$ random variable $(X_{1},X_{2},...,X_{n^{2}})$, with each variable assuming two discrete values. A Markov chain is executed to obtain a sizable sample of random variables, and we subsequently approximate this $n^{2}$ dimensional distribution using a diffusion bridge model. However, since we are using a continuous model via SDE to obtain $Y_{T}$, which cannot accurately approximate a discrete random variable with values of $\left\{0,1\right\}$, we employ the sigmoid function in the output $Y_{T}$. The loss for $X_{T}$ is then computed. Finally, when using the diffusion model, we apply post-processing to obtain the output value, i.e., $torch.round()$. In the case of the Feynman-Kac model, we set the boundary conditions to $u(x,T)=p(H(round(sigmoid(x))))$, where $p$ is $exp(-\beta/2*(x))$. In other words, we set the composite function $p(H(round(sigmoid))$ to $f$ in the boundary.

Table 2 and Table 3 are one table. We have separated them for ease of presentation, and they have the same rows. In Table 2 and Table 3, $wi,vi,z$ represent the values of the corresponding $\mathbb{E}F,\mathbb{E}G,Q$ estimated using the corresponding estimators, respectively. $true_wi,true_vi,true_z$ indicate the corresponding true values. The $error_wi,error_vi,error_z$ represent the squared error using the corresponding estimators. The terms $w_{i}$ sample points and $v_{i}$ sample points refer to the number of sampled points utilized by the estimator. The terms $w_{i}$ time and $v_{i}$ time refer to the time taken by the estimator, measured in seconds. MCMC method we employed to generate samples follows the same approach as used in [https://github.com/zysophia/Doubly_Adaptive_MCMC](https://github.com/zysophia/Doubly_Adaptive_MCMC).

At the same time we compare with the method RelMeanEst in []. MCMC-C is the method RelMeanEst, MCMC-R is the empirical mean taken using the samples obtained from resampling, and MCMC-T is the estimate of the expectation obtained using the established diffusion bridge. And the number of data points used indicates the number of points in the Markov chain used. To be fair, we lower the threshold in MCMC-C to reduce its algorithmic complexity. Because only a small number of sample points are used in MCMC-R and MCMC-T. sample points means the number of points sampled from the Markov chain. Note that when $n\geq 6$ is in the MCMC-C method due to the larger complexity we do not discuss it. We only compare MCMC-R and MCMC-T. Note: GPU types: the first of these uses Tesla P100 while the second uses Tesla V100 when $n\geq 6$. The two methods are shown in Table 2 and Table 3. Above the horizontal is the first method below the second method. We can find the performance of PINN. In the high-dimensional case ($d=n^{2}=225$).

**Table 2: Comparison of different MCMC Expectation Estimator**
| Method | $n$ | $wi$ | $vi$ | $z$ | $true_wi$ | $true_vi$ | $true_q$ |
| --- | --- | --- | --- | --- | --- | --- | --- |
| MCMC-C | 2 | 0.9706396 | 1.0306606 | 1.0618365 | 0.9654024 | 1.0357122 | 1.072778 |
| MCMC-R | 2 | 0.9550395 | 1.0308957 | 1.0794273 | 0.9654024 | 1.0357122 | 1.072778 |
| MCMC-T | 2 | 0.9626546 | 1.0333116 | 1.073398 | 0.9654024 | 1.0357122 | 1.072778 |
| MCMC-C | 3 | 0.9340726 | 1.0744393 | 1.1502738 | 0.9226402 | 1.0834867 | 1.174333 |
| MCMC-R | 3 | 0.9269992 | 1.0774463 | 1.1622948 | 0.9226402 | 1.0834867 | 1.174333 |
| MCMC-T | 3 | 0.9283546 | 1.0795156 | 1.1628268 | 0.9226402 | 1.0834867 | 1.174333 |
| MCMC-C | 4 | 0.8844253 | 1.1470378 | 1.2969301 | 0.8641533 | 1.1563625 | 1.338233 |
| MCMC-R | 4 | 0.8686283 | 1.159192 | 1.3345087 | 0.8641533 | 1.1563625 | 1.338233 |
| MCMC-T | 4 | 0.8692993 | 1.1552249 | 1.328915 | 0.8641533 | 1.1563625 | 1.338233 |
| MCMC-R | 2 | 0.9950697 | 1.00498 | 1.0099593 | 0.9654024 | 1.0357122 | 1.072778 |
| MCMC-T | 2 | 0.9735975 | 1.0444663 | 1.0727907 | 0.9654024 | 1.0357122 | 1.072778 |
| MCMC-R | 3 | 0.9949918 | 1.0049603 | 1.0100187 | 0.9226402 | 1.0834867 | 1.174333 |
| MCMC-T | 3 | 0.9188372 | 1.0843412 | 1.1801233 | 0.9226402 | 1.0834867 | 1.174333 |
| MCMC-R | 4 | 0.9951742 | 1.0050679 | 1.0099418 | 0.8599499 | 1.1563472 | 1.344668 |
| MCMC-T | 4 | 0.8664092 | 1.1570783 | 1.3354871 | 0.8599499 | 1.1563472 | 1.344668 |
| MCMC-R | 6 | 0.9949221 | 1.0049597 | 1.0100888 | 0.7163408 | 1.3985122 | 1.9523 |
| MCMC-T | 6 | 0.6953082 | 1.3970394 | 2.0092378 | 0.7163408 | 1.3985122 | 1.9523 |
| MCMC-R | 8 | 0.9950855 | 1.0050602 | 1.010024 | 0.5468445 | 1.8348543 | 3.3553494 |
| MCMC-T | 8 | 0.5683886 | 1.9384431 | 3.4104189 | 0.5468445 | 1.8348543 | 3.3553494 |
| MCMC-R | 10 | 0.9949943 | 1.0050541 | 1.0101104 | 0.3853279 | 2.60382 | 6.7574135 |
| MCMC-T | 10 | 0.3684352 | 2.833073 | 7.6894751 | 0.3853279 | 2.60382 | 6.7574135 |
| MCMC-R | 15 | 0.9949888 | 1.0050285 | 1.0100903 | 0.1135434 | 8.894777 | 78.3381355 |
| MCMC-T | 15 | 0.1181741 | 10.4130456 | 88.1161667 | 0.1135434 | 8.894777 | 78.3381355 |

**Table 3: Comparison of different MCMC Expectation Estimator**
| $error_wi$ | $error_vi$ | $error_q$ | $wi$ sample points | $vi$ sample points | $wi$ time (s) | $vi$ time (s) |
| --- | --- | --- | --- | --- | --- | --- |
| 2.74E-05 | 2.55E-05 | 0.000119716 | 3157 | 3157 | 0.29448 | 0.28554 |
| 0.00010739 | 2.32E-05 | 4.42E-05 | 100 | 100 | 14.421 | 9.702 |
| 7.55E-06 | 5.76E-06 | 3.84E-07 | 100 | 100 | 11.671 | 11.789 |
| 0.0001307 | 8.19E-05 | 0.000578845 | 30700 | 30700 | 3.40208 | 3.43238 |
| 1.90E-05 | 3.65E-05 | 0.000144918 | 2000 | 2000 | 20.163 | 19.544 |
| 3.27E-05 | 1.58E-05 | 0.000132393 | 2000 | 2000 | 236.109 | 233.225 |
| 0.000410954 | 8.70E-05 | 0.00170593 | 11383 | 11383 | 0.922 | 0.93 |
| 2.00E-05 | 8.01E-06 | 1.39E-05 | 2000 | 2000 | 26.917 | 26.696 |
| 2.65E-05 | 1.29E-06 | 8.68E-05 | 2000 | 2000 | 284.784 | 285.609 |
| 0.000880149 | 0.000944468 | 0.003946189 | 500 | 500 | 7.562 | 5.728 |
| 6.72E-05 | 7.66E-05 | 1.61E-10 | 2000 | 2000 | 15.809 | 15.456 |
| 0.005234754 | 0.006166395 | 0.026999189 | 500 | 500 | 7.748 | 7.547 |
| 1.45E-05 | 7.30E-07 | 3.35E-05 | 2000 | 2000 | 32.482 | 32.562 |
| 0.018285611 | 0.022885427 | 0.112041629 | 500 | 500 | 10.762 | 10.573 |
| 4.17E-05 | 5.35E-07 | 8.43E-05 | 2000 | 2000 | 59.597 | 58.064 |
| 0.077607541 | 0.15488357 | 0.887761945 | 500 | 500 | 25.194 | 23.99 |
| 0.00044237 | 2.17E-06 | 0.003241913 | 2000 | 2000 | 302.916 | 303.36 |
| 0.200919994 | 0.688558248 | 5.500551232 | 500 | 500 | 64.449 | 62.8 |
| 0.000464148 | 0.010730639 | 0.00303265 | 2000 | 2000 | 516.126 | 509.603 |
| 0.371693119 | 2.556052403 | 33.03149292 | 500 | 500 | 91.826 | 90.796 |
| 0.000285363 | 0.052556938 | 0.868738826 | 2000 | 2000 | 889.684 | 906.444 |
| 0.776945993 | 62.24813139 | 5979.626574 | 500 | 500 | 165.086 | 166.38 |
| 2.14E-05 | 2.305139542 | 95.60989415 | 2000 | 2000 | 1968.976 | 1997.663 |

### 9.2 Effects of diffusion bridge model

This method is also applicable for estimating integrals in high dimensions, particularly in tandem with the Monte Carlo method, when the target distribution is easily samplable. We conducted a comparison in estimating the distribution of the bridge using a bridge constructed from partially sampled high-quality samples. This approach enables continuous sampling of the target distribution by utilizing a well-established bridge. To illustrate, we simulated a diffusion bridge model (DBM) to approximate the distribution of a target variable $Y=(X_{1},X_{2},X_{3})$, where $X_{1}\sim N(1,2)+Beta(4,2)$, $X_{2}\sim N(-1,2)+Gamma(1,2)$, and $X_{3}\sim N(3,2)+geometric(0.5)$. We sampled 500 points from the target distribution and employed DBM matching to obtain an SDE. Subsequently, we compared the distribution of the generated tracks to the target distribution. Continuing the target distribution sampling using the constructed bridge, we sampled an additional 500 points and compared the differences between the resampled samples and the original target distribution. The specific parameters include $T=0.2$, time step size $h=0.025$. The DBM is trained for 300 epochs with a learning rate of 0.001, using the Adam optimizer and Wasserstein distance as the loss. This method facilitates the construction of a pair of target distributions amenable to sampling. The expectation $\mathbb{E}(f(X))$ of the target distribution can be obtained by utilizing FCM.

Figure: Figure 1: Comparison of the probability density functions of the generated and resampled paths and target distributions for each dimension. Two of the blue lines are the mean of the experience of the target sample and the mean of the experience of the re-generation sample, respectively.
Refer to caption: extracted/5704698/fig5.png

Figure: Figure 2: Generated tracks
Refer to caption: extracted/5704698/fig4.png

### 9.3 Other baseline Experiments

Since there are too many variants for the MCMC sampler, and our aim in this paper is to estimate the expectation rather than focusing on the selection aspect of the sampler, we consider one of the simplest LMCMC (Langevin diffusion model). It is worth noting, however, that we are using the unadjusted LMCMC here.

$$ $dX_{t}=b(x)dt+dW_{t},$ $$

where $b(x)=\frac{1}{2}\nabla_{x}\log p(x)$ and SDE solver is

$$ $X_{t+h}=X_{t}+b(X_{t})h+W_{t+h}-W_{t}\quad X_{0}=x_{0}$ $$

In all the experiments below, the interpretation of parameters is as follows: Total points in time: $D$,
Initial value: $X_{0}=x_{0}$, Brownian motion: $W_{t}$, Time Series: $t_{0},t_{1},\dots,t_{D}=T$, Euler-Maruyama method step size: $h$, Number of paths simulated: $N$

We consider the MCMC algorithm and our method to sample the density function of the target and obtain the corresponding expectation. In the MCMC algorithm configuration, we use the Langevin MCMC to get independent samples. Here we use only the value of $X_{T}$ at the terminal moment to estimate the expectation. We use the same paths in LDM+FCM, but with a different way of computing expectations.

As an illustration, we consider a one-dimensional SDE, where we define the target distribution as $p(x)=C\exp(\frac{-(x-1)^{2}}{2})$, corresponding to the drift coefficient of the LDM being $\mu(x)=\frac{1-x}{2}$. We evaluate the expectation of $\mathbb{E}(X_{10})$. To decrease the error of the Euler-Maruyama method, we use a small step size of $h=0.01$ and iterate 1000 steps to obtain the final path. We repeat the experiment $M=30$ times. During the training process, we extract points from each path every 100 points and add them to the training process, instead of using all the points on the path.

We examine an extreme case, employing a very limited number of paths ($N=5$) to estimate the true expectation $\mathbb{E}(X_{10})$. In Figure 1, we present the empirical distributions obtained through two different methods. The results obtained by LDM+FCM outperform Langevin MCMC, validating that paths can offer more informative outcomes. By incorporating gradient information from path points and integrating it into PINN for training, our method demonstrates lower variance under the same experimental configuration, significantly enhancing the efficiency of the MCMC algorithm with appropriate optimization. Although we utilize Unadjusted Langevin MCMC, our method provides unbiased estimates. This is attributed to the fact that the bias in Unadjusted Langevin MCMC stems from the numerical SDE solver, while our method does not necessitate high accuracy in $X_{t}$; we are more concerned with the precision of the corresponding $(b,\sigma)$ on $X_{t}$. Unlike direct sampling using the SDE method, which requires a highly precise SDE solver [], such precision is unnecessary in our method. We only require accurate estimations at each point on the path for the coefficients of the drift and diffusion terms.

Figure: Figure 3: The empirical distribution of $\mathbb{E}_{estimated}(X_{10})$
Refer to caption: extracted/5704698/fig2.png

For other cases, which can be handled by our method, we add an example of a broader computation of the expectation of a stable distribution in the absence of the corresponding convergence result for MCMC. for example:

$$ $dX_{t}=\frac{1}{2}h^{2}\frac{1-2X_{t}}{X_{t}^{\frac{1}{2}}(1-X_{t})^{\frac{1}{ 2}}}dt+2h{X_{t}^{\frac{1}{4}}(1-X_{t})^{\frac{1}{4}}}dW_{t}$ $$

where $X_{0}=0.5,\quad\mathbb{E}X_{1}=0.5$

Figure: Figure 4: The empirical distribution of $\mathbb{E}_{estimated}(X_{1})$
Refer to caption: extracted/5704698/fig3.png

In this example, we use this method to estimate mathematical expectations in high dimensions. We consider a normal distribution with independent and identical marginal distributions as follows: $p(x)=C\exp(-0.5(x-0.2)^{2})$ for each dimension. When $g(x_{1},\dots,x_{d})=x_{1}+\dots+x_{d}$, We use a smaller number of paths ($N=50,100$). We use the Euler-Maruyama method with a step size of 0. 1 for 100 iterations and calculate the internal loss function of the PDE every 10 points. In the case of $d=5,10$ ,we employ a 2-layer neural network with 108 units per layer and a tanh activation function. In the case of $d=20$, we use a 2-layer neural network with 526 units per layer, set $N=100$, and compute the internal loss function of the PDE every 20 points. We also repeate the experiment $M=30$ times by using different random number seeds and measure the average time required to estimate the mathematical expectation each time with GPU (Tesla P100). In the training process we train 400 epochs by using the Adam optimizer with a learning rate of 0.001.

We compute $\mathbb{E}(g(X_{1},X_{2},\dots,X_{d}))$ where $X_{i}\sim N(0.2,1)$ and estimate its error.
The errors we use is

$$ $\textbf{Absolute value error}=\frac{1}{M}\sum_{i=1}^{M}|\mathbb{E}^{i}_{ estimated}(g(X_{1},X_{2},\dots,X_{d}))-\mathbb{E}(g(X_{1},X_{2},\dots,X_{d})))|$ $$

and

$$ $\textbf{Square Error}=\frac{1}{M}\sum_{i=1}^{M}|\mathbb{E}^{i}_{estimated}(g(X _{1},X_{2},\dots,X_{d}))-\mathbb{E}_{mean}(g(X_{1},X_{2},\dots,X_{d})))|^{2}$ $$

where

$$ $\mathbb{E}_{mean}(g(X_{1},X_{2},\dots,X_{d})))=\frac{1}{M}\sum_{i=1}^{M} \mathbb{E}^{i}_{estimated}(g(X_{1},X_{2},\dots,X_{d}))$ $$

The method is LDM+FCM and we compare the results of this method with those obtained by the Langevin MCMC (LMCMC in short).

**Table 4: Comparison of different methods**
| Method | Dimension($d$) | paths($N$) | Absolute value error | Square Error | GPU time |
| --- | --- | --- | --- | --- | --- |
| LMCMC | 5 | 50 | 2. 927620e-01 | 1.253495e-01 | $\times$ |
| LDM+FCM | 5 | 50 | 1. 031084e-01 | 1.600998e-02 | 29. 62s |
| LMCMC | 10 | 50 | 4. 696985e-01 | 3.077161e-01 | $\times$ |
| LDM+FCM | 10 | 50 | 3. 330310e-01 | 1.382318e-01 | 46. 79s |
| LMCMC | 20 | 100 | 3. 630368e-01 | 1.863313e-01 | $\times$ |
| LDM+FCM | 20 | 100 | 2. 959023e-01 | 1.042063e-01 | 49. 74s |

### 9.4 Potential applications and future work

The independence of samples: Obtaining accurate expectations with entirely unknown sample independence remains a significant challenge in the real world, particularly in stochastic optimization algorithms or loss functions, where independent sampling is frequently required for estimation. The independence of samples plays a crucial role in machine learning, and its violation can significantly impact the performance and validity of machine learning models. Many machine learning models rely on the assumption of independent and identically distributed (i.i.d) samples. Non-independent samples can introduce dependencies that the model may mistakenly learn as patterns. Nonlinear mathematical expectations play a critical role in such non-iid scenarios []. However, methods like using Max-Mean Monte Carlo for calculating nonlinear mathematical expectations are often challenging. This is because we need to partition the dataset into parts where the samples are independent and then calculate the linear mathematical expectation for each part. Finally, we take the largest to get the nonlinear mathematical expectation. Our approach provides a completely new way to consider the use of Stochastic Differential Equations (SDEs) with G-Brownian motions. The diffusion bridge model is constructed using the same method and then solved directly using the Feynman-Kac model in the case of nonlinear mathematical expectations. This avoids problems such as data grouping.

Representation learning and Distributional regression learning: In the theory of statistical learning, we assume $X\sim P_{X}$ and $Y\sim P_{Y}$. A basic loss function is $l=\mathbb{E}[h_{\theta}(X)-Y)^{2})]$ and $l=\mathbb{E}[\textbf{Corss Entropy}[h_{\theta}(X),Y)]$ where $h_{\theta}$ is model, and we often need to sample a portion of the sample $\{x_{i},y_{i}\}^{N}_{i=1}$, and then optimise the empirical loss function $l=\frac{1}{N}\sum_{i=1}^{N}(h_{\theta}(x_{i})-y_{i})^{2}$ and $l=\frac{1}{N}\sum_{i=1}^{N}\textbf{Corss Entropy}(h_{\theta}(x_{i}),y_{i})$. But in the case where the sample size does not fully cover the distribution of the corresponding totality, because the loss function is obtained by sampling a portion of the dataset, the loss function that we obtain tends to be biased, or has a large variance. When we have a high quality diffusion bridge that can accurately approximate the distribution of the target $(P_{X},P_{Y})$, which most of the current diffusion bridge models can do. We can achieve this by configuring the boundary conditions in the Feynman-Kac model to be $f(x,y)=(h_{\theta}(x)-y)^{2}$ and $\textbf{Cross Entropy}(h_{\theta}(x),y)$. We then replace the empirical loss function with the PDE loss and the PDE loss at the boundary. This approach may enable us to enhance the learning of the Representation of a Distribution. This is because the diffusion bridge model captures information about the entire distribution rather than just the local distribution of specific points. When estimating expectations, we incorporate the PDE loss function, which contains gradient information regarding the diffusion bridge coefficients. The coefficients of the diffusion bridge tend to exhibit correlations with the target distribution. In this case, the number of points required for the diffusion bridge coefficients is often significantly smaller than the number of points $N$ directly sampled from the data. Finally, we can use the trained diffusion bridge model to perform some basic statistical learning tasks.

Variational Inference: Due to the extensive application of mathematical expectations in machine learning and probabilistic statistics, we are unable to comprehensively demonstrate all relevant methods in this paper. We will consider applying these methods to important domains, such as estimating the evidence lower bound (ELBO) in Black-Box Variational [] Inference. We often need to use reparameterization techniques to estimate the $\textbf{ELBO}=\mathbb{E}_{q(z|\phi)}\log p(x,z)-\log q(z|\phi)$ with small bias, but we can consider using a diffusion bridge to approximate the target distribution $q(z|\phi)$, and select $f(z)=\log p(x,z)-\log q(z|\phi)$, where $\phi$ can be designed as a trainable parameter. In this way, we can modify our optimization objective from ELBO to $-u(x_{0},t_{0})+\texttt{PDE loss}+\texttt{boundary loss}$, which can achieve lower variance and GPU acceleration.

## 10 Algorithms

For $X$ taking values in $\mathbb{R}^{d}$, if the marginal distribution is not independent, we employ one-dimensional Wasserstein distances []. In this case the Sinkhorn algorithm [] can be used to address the optimal transport problem in $d$-dimensions.

Backpropagation: By encapsulating the computation of the 2-Wasserstein algorithm into an nn.Module and implementing it using PyTorch, we can retain the computation graph during the calculation. This allows us to obtain precise gradients using automatic differentiation methods.

Figure: Algorithm 3 Diffusion bridge model (DBM)

The algorithm of the Feynman-Kac model is similar to that of PINN. It is mainly a matter of using the diffusion coefficients obtained earlier and solving the corresponding PDE for the data points. Not all points on the paths need to be included in the training in this algorithm. This is the same as the training of PINN, where we only need to sample a fraction of the points to get the solution.

Figure: Algorithm 4 Feynman-Kac model (FCM)

## Peer Reviews

confidence: 4
contribution: 1
initial_rating: 1
presentation: 1
questions: In Theorem 2.1, what does "Linear growth" mean?

In Assumption 2.4, what does "D is the metric of the parameter" mean?

In Theorem 2.6, should "we exist" be "there exists a set"?

In Algorithm 1, can the authors clarify what is the main difference between $X_t$ and $X_i$? The distinction is not clear to me. 

How are the integrals in Algorithm 2 computed? Are the authors able to evaluate the integrals explicitly? If so, they should indicate how and why they are able to do so.

How does the computational cost of this approach compare to other MCMC approaches?

One of the main criticisms posed about MCMC methods is that they are not optimal since they spend quite some time in the burn-in phase (lines 56, 82). However, for the proposed method to work the authors assume that they have access to samples from a distribution that are obtained via an MCMC, that has already gone through a burn-in period (lines 188, 192). How can the authors support their claim that this method is better (line 107) than MCMC if it is still spending a similar amount of samples in burn-in?
soundness: 1
strengths: The goals stated in the introduction are bold, and quite interesting. Obtaining results in this line would prove quite useful in general for ML and statistics.

I appreciate that the authors include small introductions for the Euler-Maruyama method and physics informed neural networks in the appendix. However their existence should be indicated in the main text.
summary: Doing MCMC is hard (time consuming, and somewhat wasteful because of the burn-in period). The authors propose a post-processing method using the samples from some MCMC procedures, which admit an Itô decomposition, to approximate moment estimates of the desired sampling density. They use a denoising technique based on physics informed neural networks as part of the post-processing mechanism.
weaknesses: Presentation is bad throughout. There are plenty of typos. The authors do not use parenthetical citations and instead insert them in the text which makes for a less pleasant reading experience. 

The notation introduced in line 232 definitely needs improvement, I do not understand which side is supposed to be the one that will be used later. Even then, it is unclear what is being defined, as there are two definitions for $\hat{\mu}\_{t\_i}$ .

In Assumption 2.2, it is not clear what $\mu^{\mathcal{P}_\theta}$ means.

The notation in Algorithms 1 and 3 should be introduced before the algorithms. For example, it is not clear where $Y_T$ comes from, and why it is required.

Table 1 is impossible for me to parse. I invite the authors to mimic the conciseness of their own Table 5 for summarizing the numerical results.

There is a link to GitHub page, which has been confirmed to not belong to the authors, but I do not see a good reason for why the authors would want to include a link to a GitHub that is not theirs. Usually a citation to the original article is enough.

The proofs in the main text for Theorems 2.1, 2.6 and 2.7 should either refer directly to the Appendix where they are proved, or be proved right there.

Regarding Theorem 2.8, the comment in the 'proof' space makes me think the Authors were not the first to prove it, in which case they should indicate it explicitly; otherwise it would be plagiarizing.

Currently Section 4 is quite lacking, including the aforementioned Table 1 which I cannot comprehend (by the way, it is missing a reasonable caption). 

A proposed method like this should be thoroughly tested, which in the current state of the paper it has not been. The methods the authors refer for comparison should include appropriate references.

---

confidence: 3
contribution: 2
initial_rating: 3
presentation: 2
questions: The authors mention in the Discussion that their method requires the boundary conditions of the PDE to satisfy a smoothness condition, specifically that f is in C^2, and that this requirement broadens the scope of their approach. However, it seems that C^2 smoothness could be more restrictive than a Lipschitz assumption. Could the authors clarify how they view this requirement as less restrictive? Additionally, could they discuss any potential limitations this might introduce for functions that are not in C^2?
soundness: 2
strengths: A significant strength of this work is the innovative linking of the diffusion model to high-dimensional partial differential equations (PDEs), with Physics-Informed Neural Networks (PINNs) effectively employed to overcome the curse of dimensionality in solving these PDEs.
summary: The authors propose the Feynman-Kac Operator Expectation Estimator (FKEE) to approximate the target distribution E[f(X)]. This estimator contains two parts: (1) A diffusion bridge model with parameters optimized to minimize the Wasserstein distance to the target distribution, and (2) a method based on the Feynman–Kac equation, formulated as a partial differential equation (PDE) and solved approximately using Physics-Informed Neural Networks (PINNs), which employ a least-squares approach. The experiments focus on approximating the partition function in a random graph model.
weaknesses: 1. The Feynman–Kac model (Algorithm 2) with the PINN solver lacks a convergence or error estimate, which would be valuable for assessing its accuracy and reliability.

2. In the experiments, the authors claim that "using fewer points on the Markov chain achieves higher accuracy in approximating expectations." However, it is unclear if this result generalizes beyond the specific example provided, as it appears quite context-dependent.

---

confidence: 4
contribution: 2
initial_rating: 5
presentation: 2
questions: 1. I don't know why and when MCMC is required to impose complex constraints on the distribution and performance function. Some references are suggested.

2. I don't see when MCMC is not the optimal decoding method. It appears to me that burn-in is not a significant limitation and only affects the performance on a negligible scale and can be easily fixed via a large stepsize in the beginning for warm-up.

3. Discussions on the limitations would be preferred.
soundness: 2
strengths: the idea of using Feyman-Kac to approximate the expectation is interesting.
summary: The authors proposed to leverage the Feynmann-Kac equation via Physically Informed Neural Networks (PINN) to approximate the target Mathematical Expectation efficiently and heuristically without causing a large variance.
weaknesses: The scalability of the algorithm w.r.t. dimension is not verified sufficiently. d=20 is too small. There are no real-world simulations.

The authors criticize the large variance issue by the MCMC method but fail to justify theoretically why the proposed method yields a lower variance. The empirical support is limited.

NIT: Theorem 2.1: the discretization error by Growall inequality is weak and exponentially dependent on time. Girsanov can be used to fix it.

---

confidence: 5
contribution: 2
initial_rating: 1
presentation: 2
questions: + What are the limitations  of your methods regarding dimension? 

+ How does your algorithm compare with other neural ODE/SDE, diffusion models, bridges models in the literature?
soundness: 2
strengths: Developing new sampling methods using ideas from generative model is an active area with lots of promising recent advances.
summary: The paper presents two generative models to use in the context of sampling: a diffusion bridge with $W^2$-loss and algorithm based on solving the Feynmann-Kac PDE using PINNS.
weaknesses: + There is no really new theoretical contributions in this paper. All the theorems mentioned in the paper are standard results. The diffusion bridge model has been used in various iterations in countless papers in the literature, for example by Doucet and collaborators.  

+ There is no indication that the algorithms presented here will scale up with dimension, especially if using the W^2 loss, so the claim that this improves on MCMC seem somewhat overblown.

+ Using a PINNs to solve the Feynmann-Kac equation is very unlikely to work in high dimension as PINNs are usually are no easy to to train.   

+ Lack of experiments on high-dimensional data sets.

+ No head-to-head comparisons with state of the art algorithms.

---

confidence: 3
contribution: 2
initial_rating: 3
presentation: 2
questions: Can the author disentangle their works from past literature on neural SDE and Feynman-Kac's use of averaging observables?

Can the authors provide evidence that their universal sampler outperforms the existing techniques, including those in Blechshmidt
et. al.(2021)? Can they provide some canonical well-excepted benchmark at which they excel over others?
soundness: 3
strengths: The work addresses an important fundamental topic. It provides an interesting combination of DNN-based techniques.
summary: The manuscript combines several DNN-based numerical techniques, associated with solving and learning SDEs and their PDE counterparts, with the aim of sampling from distributions and averaging over observables. Specifically, a type of neural-SDE is put forward and trained to converge to the desired distribution. Next averages of functions are estimated using the Feynman-Kac formula, wherein a PDE is derived from the learned SDE which “calculates” the averages for us given that one sets the boundary condition to be the function one wishes to average. This complex problem is solved using a PINN approach. Some numerical results on random graphs are reported in a table.

While the paper suggests an interesting technical path to explore in this important subdomain MC samplers— it is not clear from reading the manuscript how much of this is incremental work stringing together a series of known results or rather an original and meaningful step forward. While I am only a causal user of Monte-Carlo, my strong feeling is that the benchmarks shown are insufficient to prove that
this specific combination of techniques outperforms existing ones (including DNNs/PINN-based ones and more advanced MC techniques). Turning to the techniques themselves, it seems that neural-SDEs preceded the current work [Tzen, Ragidsky 2019], and that using Feynman-Kac formula with a combination of DNNs to estimate observables on an SDE has also been done before [Blechshmidt, Ernst 2021 and refs therein]. Reading into the latter work, it seems that DNNs have been used slightly differently than in the current work to solve the PDE associated with the SDE, but is this a conceptual change given the common use PINNs to solve PDE? Does it hold the key to any SOTA results? The related work section (which has been delegated to the appendix) and the general causal referencing to related works (such as Tzen et. al.) leave the non-expert reader with little understanding of the true novelty of the current results.

The presentation of the work also leaves a gap between conceptual claims and practical contributions. It also feels fragmented and the
common use of signposts and bold notation only worsens this in my mind. The conceptual claims, which are sometimes grand, are hard to substantiate. For instance, “Establishing a Link Between Sampling Methods and High-Dimensional Partial Differential Equations”, taken at face value, can hardly be attributed to the current work with all the knowledge on SDE, Fokker-Plank Equations, and Feynman-Kac formula. Also “Expanding the Scope of Expectation Estimators“, feels vague. Is there a current scope of expectation estimators? What results in the current work expand this scope in a way that others can't? Finally, in their introduction, the authors allude to the fact that the authors have an affirmative answer to the question “Is it possible to unify most existing MCMC algorithms into a cohesive framework to create a universal sampler for expectation estimation?“--- This is such a rich and complex problem that providing an affirmative answer would clearly violate various prevalent complexity theory assumptions. For instance, can the authors show that their sampler solves the Ising Spin-Glass problem? Can the authors even solve the much simpler case of the 2d Ising model at the phase transition and compute long-range observables? Does their technique outperform various ones used in physics to overcome sampling problems such as replica-exchange Monte Carlo?

The above issues, concerning the entanglement with previous works, evidence of going beyond SOTA, and its portrayed grand scope, prevent me from recommending it for publication in ICLR.
weaknesses: Relation to past literature is left vague. It does not show SOTA on a sufficient set of canonical problems with a suitable comparison to
other techniques. Its claims feel too broad and general. The writing style can be improved.

Meta-review: **Summary of Discussion:**  
The reviewers appreciated the novel ideas proposed in the paper, such as using the Feynman-Kac equation (FKE) and Physics-Informed Neural Networks (PINNs) to estimate mathematical expectations. However, the submission has significant issues that hinder its acceptance:  

1. **Lack of Novelty and Theoretical Depth:** The reviewers noted that the paper primarily combines existing techniques (diffusion bridges, PINNs, and MCMC sampling). While innovative connections were drawn, the theoretical contributions were incremental and lacked depth. Key claims, such as reduced variance and scalability to high dimensions, were not sufficiently justified either empirically or theoretically.

2. **Insufficient Empirical Evidence:**  
   - The experiments were limited to relatively low-dimensional settings (e.g., \(d=20\) or Ising model with \(15^2\) sites).  
   - There was a lack of real-world applications or benchmarks against state-of-the-art methods on well-established datasets.  
   - Standard deviations and error analysis were missing, leaving readers unable to assess the robustness of results.  

3. **Overstated Claims:**  
   - Grand claims, such as “universal samplers” and “expanding the scope of expectation estimators,” were not well-substantiated.  
   - Connections to broader MCMC frameworks were vague and not explicitly supported by evidence.  

4. **Presentation and Clarity Issues:**  
   - Poor organization made the paper difficult to follow. Important notations and algorithms were introduced without sufficient explanation.  
   - The related work section, relegated to the appendix, left readers unclear about the novelty of the contributions relative to prior work.  

5. **Rebuttal Limitations:** While the authors clarified some issues during the discussion, fundamental concerns about novelty, experimental rigor, and clarity remain unresolved. Addressing these would require a substantial overhaul of the paper.  

**Conclusion:**  
The paper addresses an important topic but lacks the necessary clarity, rigor, and empirical support to meet the standards of ICLR. A thorough revision with a stronger theoretical foundation, robust experiments, and clearer presentation is recommended.

Decision: Reject