Diffusion-PINN Sampler

## Abstract

Abstract Recent success of diffusion models has inspired a surge of interest in developing sampling techniques using reverse diffusion processes.
However, accurately estimating the drift term in the reverse stochastic differential equation (SDE) solely from the unnormalized target density poses significant challenges, hindering existing methods from achieving state-of-the-art performance.
In this paper, we introduce the Diffusion-PINN Sampler (DPS), a novel diffusion-based sampling algorithm that estimates the drift term by solving the governing partial differential equation of the log-density of the underlying SDE marginals via physics-informed neural networks (PINN).
We prove that the error of log-density approximation can be controlled by the PINN residual loss, enabling us to establish convergence guarantees of DPS.
Experiments on a variety of sampling tasks demonstrate the effectiveness of our approach, particularly in accurately identifying mixing proportions when the target contains isolated components.

## 1 Introduction

Sampling from unnormalized distributions is a fundamental yet challenging task encountered across various scientific disciplines such as Bayesian statistics, computational physics, chemistry, and biology (Liu & Liu, ; Stoltz et al., ).
Markov chain Monte Carlo (MCMC) and variational inference (VI) have historically been the go-to methods for this problem. However, these approaches exhibit limitations when dealing with complex target distributions (e.g., distributions with multimodality or heavy tails).
Recently, the success of diffusion models for generative modeling (Song et al., ; Ho et al., ; Nichol & Dhariwal, ; Kingma et al., ) have sparked considerable interest in tackling the sampling problem using the reverse diffusion processes that transport a given prior density to the target, governed by stochastic differential equations (SDE).

In diffusion-based generative models, the score function in the drift term of the reverse SDE is learned based on score matching techniques (Hyvärinen & Dayan, ; Vincent, ) that require samples from the target data distribution.
However, for sampling tasks, we only have access to an unnormalized density function $\pi$, making it challenging to estimate the score function for the reverse SDE.
From a stochastic optimal control perspective (Tzen & Raginsky, ; De Bortoli et al., ), several VI methods that parameterize the drift term with neural network approximation have been proposed (Zhang & Chen, ; Berner et al., ; Vargas et al., , ).
Nevertheless, these approaches face challenges such as instability during training, the computational complexity associated with differentiating through SDE solvers, and mode collapse issues arising from training objectives based on reverse Kullback-Leibler (KL) divergences (Zhang & Chen, ; Vargas et al., ).
On the other hand, Huang et al. () proposed a scheme based on the connection between score matching and non-parametric posterior mean estimation.
More specifically, they use MCMC estimation of the scores to potentially alleviate the numerical bias intrinsic in parametric estimation methods such as neural networks.
However, this method also introduces noise in the estimates and requires repetitive posterior sampling in each time step of the reverse SDE.
Overall, despite their potential, diffusion-based sampling methods have not yet achieved state-of-the-art performance.

In addition to its connection with posterior mean estimation, the score function has also been shown to evolve according to a partial differential equation known as the *score Fokker Planck equation* (score FPE) (Lai et al., ).
This discovery has led to a novel regularization technique for enhancing score function estimation in diffusion models (Lai et al., ; Deveney et al., ).
In this paper, we adopt this strategy for diffusion-based sampling methods.
While the score function can be recovered by solving the score FPE using the score of target distribution $\pi$ as the initial condition, we demonstrate that it may fail to identify correct mixing proportions when $\pi$ has isolated components, a common limitation of score-based methods (Wenliang, ; Zhang et al., ).
To remedy this issue, we propose to solve the log-density FPE, a similar partial differential equation for the log-density function, using physics-informed neural networks (PINN) (Raissi et al., ; Wang et al., ).
The estimated log-density function is then integrated into the reverse SDE, leading to a novel sampling algorithm termed *Diffusion-PINN Sampler* (DPS).
We prove that the error of log-density estimation can be controlled by the PINN residual loss, which allows us to ensure convergence guarantee of DPS based on established results for score-based generative models (Chen et al., , ; Benton et al., ).
Experiments on a variety of sampling tasks provide compelling numerical evidence for the superiority of our method compared to other baseline methods.

## 2 Related Works

Recently, several works have explored the combination of Physics-Informed Neural Networks (PINN) and sampling techniques.
For instance, Máté & Fleuret (); Fan et al. (); Tian et al. () address the continuity equation using PINN based on ODEs and achieve flow-based sampling through a linear interpolation (i.e., annealing) path between the target distribution and a simple prior, such as a Gaussian distribution.
Besides, Berner et al. () (in the appendix of their paper) and Sun et al. () propose solving the log-density Hamilton–Jacobi–Bellman (HJB) equation via PINN to develop a SDE-based sampling algorithm.
However, both approaches lack comprehensive numerical investigation and thorough theoretical analysis.
In contrast, our work investigates a limitation of score-based Fokker-Planck equations (FPE) in identifying the mixing proportions of multi-modal distributions, introduces novel computational techniques for solving PDEs via PINN in the context of diffusion-based sampling, and provides the first complete theoretical analysis of the algorithm.

## 3 Background

##### Notations.

Throughout the paper, $\Omega\subset\mathbb{R}^{d}$ denotes a bounded and closed domain.
For simplicity, we do not distinguish a probabilistic measure from its density function.
We use $\bm{x}=(x_{1},\cdots,x_{d})^{\prime}$ to denote a vector in $\mathbb{R}^{d}$ and $\|\bm{x}\|=\sqrt{x_{1}^{2}+\cdots x_{d}^{2}}$ stands for the $L^{2}$-norm.
Let $\nu$ denote a probability measure on $\mathbb{R}^{d}$, for any $\bm{f}:\mathbb{R}^{d}\to\mathbb{R}^{m}$, we denote $\|\bm{f}(\cdot)\|^{2}_{L^{2}(\Omega;\nu)}:=\int_{\Omega}\|\bm{f}(\bm{x})\|^{2}\mathrm{\leavevmode\nobreak\ d}\nu(\bm{x})$.
For any $\bm{f}:\mathbb{R}^{d}\times[0,T]\to\mathbb{R}^{m}$ , we define $\|\bm{f}_{t}(\cdot)\|_{L^{2}(\Omega;\nu)}^{2}:=\int_{\Omega}\|\bm{f}_{t}(\bm{x})\|^{2}\mathrm{\leavevmode\nobreak\ d}\nu(\bm{x})$ as a function of $t\in[0,T]$.
For any $\bm{F}=(F_{1},\cdots,F_{d})^{\prime}:\mathbb{R}^{d}\to\mathbb{R}^{d}$, we denote the divergence of $\bm{F}$ by $\nabla\cdot\bm{F}:=\sum_{i=1}^{d}\partial_{x_{i}}F_{i}$.
For any $F:\mathbb{R}^{d}\to\mathbb{R}$, we denote the Laplacian of $F$ by $\Delta F:=\sum_{i=1}^{d}\partial^{2}_{x_{i}}F$.

##### Diffusion models.

In diffusion models, noise is progressively added to the training samples via a forward stochastic process described by the following stochastic differential equation (SDE)

$$ $\mathrm{\leavevmode\nobreak\ d}\bm{x}_{t}=\bm{f}(\bm{x}_{t},t)\mathrm{\leavevmode\nobreak\ d}t+g(t)\mathrm{\leavevmode\nobreak\ d}\bm{B}_{t},\quad\bm{x}_{0}\sim p_{0}(\cdot),\quad 0\leqslant t\leqslant T,$ (1) $$

where $p_{0}(\cdot)$ is the data distribution, $\bm{B}_{t}$ is a standard Brownian motion, and $\bm{f}(\bm{x}_{t},t)$ and $g(t)$ are the drift and diffusion coefficients respectively.
The derivatives of the log-density of the forward marginals, i.e., *scores*, are learned via score matching techniques (Vincent, ; Song et al., ) and new samples from the data distribution can be obtained by simulating the following reverse process

$$ $\mathrm{\leavevmode\nobreak\ d}\bm{x}_{t}=\left[\bm{f}(\bm{x}_{t},t)-g^{2}(t)\nabla_{\bm{x}_{t}}\log p_{t}(\bm{x}_{t})\right]\mathrm{\leavevmode\nobreak\ d}t+g(t)\mathrm{\leavevmode\nobreak\ d}\bar{\bm{B}}_{t},\quad\bm{x}_{T}\sim p_{T}(\cdot),$ (2) $$

where $p_{t}(\cdot)$ is the probability density of $\bm{x}_{t}$ and $\bar{\bm{B}}_{t}$ is a standard Brownian motion from $T$ to 0 0.

##### Physics-informed neural networks (PINN).

PINN is a deep learning method for solving partial differential equations (PDEs) (Raissi et al., ).
Consider the following general form of PDE

$$ $\displaystyle\mathcal{L}u(\bm{x})=\varphi(\bm{x}),\quad$ $\displaystyle\bm{x}\in\Omega\subseteq\mathbb{R}^{d},$ (3a) $\displaystyle\mathcal{B}u(\bm{x})=\psi(\bm{x}),\quad$ $\displaystyle\bm{x}\in\partial\Omega,$ (3b) $$

where $\mathcal{L}$ and $\mathcal{B}$ are the differential operators on domain $\Omega$ and boundary $\partial\Omega$, respectively.
PINN seeks an approximate solution using deep model $u_{\theta}(\bm{x})$ by minimizing the $L^{2}$ PINN residual losses

$$ $\displaystyle\ell_{\Omega}(u_{\theta})$ $\displaystyle:=\left\|\mathcal{L}u_{\theta}(\bm{x})-\varphi(\bm{x})\right\|^{2}_{L^{2}(\Omega;\nu)},$ (4a) $\displaystyle\ell_{\partial\Omega}(u_{\theta})$ $\displaystyle:=\left\|\mathcal{B}u_{\theta}(\bm{x})-\psi(\bm{x})\right\|^{2}_{L^{2}(\Omega;\nu)},$ (4b) $$

where $\nu$ is a probability measure for collocation points generation, often taken to be the uniform distribution on $\Omega$.
The two terms $\ell_{\Omega}(u)$ and $\ell_{\partial\Omega}(u)$ in Eq. ([4](#S3.E4)) reflect the approximation error on $\Omega$ and $\partial\Omega$ respectively.
In practice, the losses in Eq. ([4](#S3.E4)) can be optimized by gradient-based methods with Monte Carlo gradient estimation.

##### Fokker Planck equation.

The evolution of the density $p_{t}(\bm{x})$ associated with the forward SDE ([1](#S3.E1)) is governed by the Fokker Planck equation (FPE) (Øksendal, )

$$ $\partial_{t}p_{t}(\bm{x})=\underbrace{\frac{1}{2}g^{2}(t)\Delta p_{t}(\bm{x})-\nabla\cdot\left[\bm{f}(\bm{x},t)p_{t}(\bm{x})\right]}_{:=\mathcal{L}_{\textrm{FPE}}p_{t}(\bm{x})}.$ (5) $$

Recently, Lai et al. () derive an equivalent system of PDEs for the log density $\log p_{t}(\bm{x})$ and score $\nabla_{\bm{x}}\log p_{t}(\bm{x})$, termed as the log-density Fokker Planck equation (log-density FPE) and the score Fokker Planck equation (score FPE) respectively, as summarized in Theorem [1](#Thmtheorem1) (the proof can be found in Appendix [A.1](#A1.SS1)).

###### Theorem 1 (Log-density FPE and score FPE; Proposition 3.1 in Lai et al. ( 2023 ) ) .

Assume the density $p_{t}(\bm{x})$ is sufficiently smooth on $\mathbb{R}^{d}\times[0,T]$. Then for all $(\bm{x},t)\in\mathbb{R}^{d}\times[0,T]$, the log-density $u_{t}(\bm{x}):=\log p_{t}(\bm{x})$ satisfies the PDE

$$ $\partial_{t}u_{t}(\bm{x})=\underbrace{\frac{1}{2}g^{2}(t)\Delta u_{t}(\bm{x})+\frac{1}{2}g^{2}(t)\left\|\nabla_{\bm{x}}u_{t}(\bm{x})\right\|^{2}-\bm{f}(\bm{x},t)\cdot\nabla_{\bm{x}}u_{t}(\bm{x})-\nabla\cdot\bm{f}(\bm{x},t)}_{:=\mathcal{L}_{\emph{{L-FPE}}}u_{t}(\bm{x})},$ (6) $$

and the score $\bm{s}_{t}(\bm{x}):=\nabla_{\bm{x}}\log p_{t}(\bm{x})$ satisfies the PDE

$$ $\partial_{t}\bm{s}_{t}(\bm{x})=\underbrace{\nabla_{\bm{x}}\left[\frac{1}{2}g^{2}(t)\nabla\cdot\bm{s}_{t}(\bm{x})+\frac{1}{2}g^{2}(t)\left\|\bm{s}_{t}(\bm{x})\right\|^{2}-\bm{f}(\bm{x},t)\cdot\bm{s}_{t}(\bm{x})-\nabla\cdot\bm{f}(\bm{x},t)\right]}_{:=\mathcal{L}_{\emph{{S-FPE}}}\bm{s}_{t}(\bm{x})}.$ (7) $$

###### Theorem 1 (Log-density FPE and score FPE; Proposition 3.1 in Lai et al. ( 2023 ) ) .

## 4 Diffusion-PINN Sampler

We consider sampling from a probability density $\pi(\bm{x})=\mu(\bm{x})/Z$ with $\bm{x}\in\mathbb{R}^{d}$, where $\mu(\bm{x})$ has an analytical form and $Z=\int_{\mathbb{R}^{d}}\mu(\bm{x})\mathrm{d}\bm{x}$ is the intractable normalizing constant.
Throughout, we only consider the forward process ([1](#S3.E1)) with an explicit conditional density of $\bm{x}_{t}|\bm{x}_{0}\sim\pi_{t|0}(\cdot|\bm{x}_{0})$.
We denote by $\pi_{t}$ the marginal density of $\bm{x}_{t}$ associated with ([1](#S3.E1)) from $\bm{x}_{0}\sim\pi_{0}=\pi$.

Inspired by diffusion models, sampling can be performed by simulating a reverse process ([8](#S4.E8)) targeting at $\pi(\bm{x})$, given an accurate estimate of the perturbed scores $\bm{s}_{\theta}(\bm{x},t)\approx\nabla_{\bm{x}}\log\pi_{t}(\bm{x})$,

$$ $\mathrm{\leavevmode\nobreak\ d}\bm{x}_{t}=\left[\bm{f}(\bm{x}_{t},t)-g^{2}(t)\bm{s}_{\theta}(\bm{x}_{t},t)\right]\mathrm{\leavevmode\nobreak\ d}t+g(t)\mathrm{\leavevmode\nobreak\ d}\bar{\bm{B}}_{t},\quad\bm{x}_{T}\sim\pi_{\textrm{prior}},$ (8) $$

where $\pi_{\textrm{prior}}$ denotes the stationary distribution of the forward process ([1](#S3.E1)) and $T$ is large enough such that $\pi_{T}\approx\pi_{\textrm{prior}}$.
However, unlike generative models, sampling tasks lack training data from $\pi$, which hinders the application of denoising score matching for perturbed score estimation.
In this section, we propose to solve the log-density FPE ([6](#S3.E6)) with PINN to estimate the perturbed scores.
While the score FPE can also be used for this purpose, we find that it may fail to learn the mixing proportions properly when the target contains isolated modes.

### 4.1 Failure of Score FPE

Consider the case where the target is a mixture of Gaussians (MoG) with two distant modes.
The following example shows that, for two MoGs with the same modes but different weights, the Fisher divergence(^*^**The Fisher divergence is defined as $F(\nu_{1},\nu_{2}):=\mathbb{E}_{\bm{x}\sim\nu_{1}}\left[\left\|\nabla_{\bm{x}}\log\nu_{1}(\bm{x})-\nabla_{\bm{x}}\log\nu_{2}(\bm{x})\right\|^{2}\right]$ for two probability measures $\nu_{1}$ and $\nu_{2}$.) between them can be arbitrarily small but the KL divergence between them remains large when the two modes are sufficiently separated.
See Figure [1](#S4.F1) (left) for an illustration of this phenomenon.
More general theoretical results can be found in Appendix [A.2](#A1.SS2).

###### Example 1 .

For any $\tau>0$, there exists $M_{\tau}(d)>0$ such that the following holds. For every $\bm{a}\in\mathbb{R}^{d}$ satisfied $\|\bm{a}\|\geqslant M_{\tau}(d)$, $w_{1},w_{2},\tilde{w}_{1},\tilde{w}_{2}\geqslant 0.1$, $w_{1}+w_{2}=1$, and $\tilde{w}_{1}+\tilde{w}_{2}=1$, MoG $\pi^{M}=w_{1}\mathcal{N}(\bm{a},I_{d})+w_{2}\mathcal{N}(-\bm{a},I_{d})$ and $\tilde{\pi}^{M}=\tilde{w}_{1}\mathcal{N}(\bm{a},I_{d})+\tilde{w}_{2}\mathcal{N}(-\bm{a},I_{d})$ satisfy

$$ $\operatorname{KL}(\pi^{M}\|\tilde{\pi}^{M})\geqslant w_{1}\log\frac{w_{1}}{\tilde{w}_{1}}+w_{2}\log\frac{w_{2}}{\tilde{w}_{2}}-\tau,\text{ but }F(\pi^{M},\tilde{\pi}^{M})<\tau,$ (9) $$

where $F(\pi^{M},\tilde{\pi}^{M})$ denotes the Fisher divergence between $\pi^{M}$ and $\tilde{\pi}^{M}$.

###### Example 1 .

#### 4.1.1 Solving Score FPE Struggles to Learn the Weights

Let $\pi^{M},\tilde{\pi}^{M}$ be the MoGs in Example [1](#Thmexample1). For any $t\in[0,T]$, $\pi^{M}_{t}$ denotes the marginal distribution of $\bm{x}_{t}$ associated with the forward process ([1](#S3.E1)) from $\bm{x}_{0}\sim\pi^{M}$. We denote $\bm{s}_{t}^{M}(\bm{x}):=\nabla_{\bm{x}}\log\pi^{M}_{t}(\bm{x})$ which is the solution to ([7](#S3.E7)) with $\bm{s}_{0}^{M}(\bm{x})=\nabla_{\bm{x}}\log\pi^{M}(\bm{x})$.
Similarly, we define $\tilde{\pi}^{M}_{t}$ and $\tilde{\bm{s}}_{t}^{M}(\bm{x})$.

Figure: Figure 1: Left: KL divergence, Fisher divergence, and log-density error between $\pi^{M}$ and $\tilde{\pi}^{M}$ as functions of $w_{1}$, where $\tilde{w}_{1}=0.2$ and $\bm{a}=(-5,-5)^{\prime}$. Middle/Right: The evolution of log-density error/Fisher divergence along the forward process respectively. The forward process achieves standard Gaussian at $t=1$.
Refer to caption: /html/2410.15336/assets/x1.png

Consider solving score-FPE ([7](#S3.E7)) using the following PINN residual loss

$$ $\ell_{\textrm{S-res}}\left(\bm{s};\bm{x},t\right):=\left\|\partial_{t}\bm{s}_{t}(\bm{x})-\mathcal{L}_{\textrm{S-FPE}}\bm{s}_{t}(\bm{x})\right\|^{2},$ (10) $$

Though $\pi^{M}$ and $\tilde{\pi}^{M}$ are equipped with different weights, their scores both satisfy the PDE ([7](#S3.E7)) such that $\ell_{\textrm{S-res}}(\bm{s}^{M};\bm{x},t)=\ell_{\textrm{S-res}}(\tilde{\bm{s}}^{M};\bm{x},t)=0$ for any $(\bm{x},t)\in\mathbb{R}^{d}\times[0,T]$.
The PINN approach, therefore, can only distinguish $\pi^{M}$ and $\tilde{\pi}^{M}$ through the initial condition.
However, Example [1](#Thmexample1) shows that the difference between $\bm{s}_{0}^{M}(\bm{x})$ and $\tilde{\bm{s}}_{0}^{M}(\bm{x})$ can be arbitrarily small, indicating the difficulty of correctly identifying the weights by solving the score FPE.
Figure [1](#S4.F1) (right) shows that the perturbed score can not tell the difference of weights until the every end of the forward process.
On the other hand, it is noticeable that the perturbed log-density distinguishes the weights well throughout the forward process (Figure [1](#S4.F1), middle).
This suggests us to solve log-density FPE and compute the scores by taking the gradient of the approximated log-density.

### 4.2 Solving Log-density FPE

To estimate the perturbed scores, we consider solving log-density FPE with initial condition:

$$ $\displaystyle\partial_{t}u_{t}(\bm{x})$ $\displaystyle=\mathcal{L}_{\textrm{L-FPE}}u_{t}(\bm{x}),$ (11a) $\displaystyle u_{0}(\bm{x})$ $\displaystyle=\log\mu(\bm{x}),$ (11b) $$

where the exact solution is $u^{*}_{t}(\bm{x})=\log\mu_{t}(\bm{x}):=\log\pi_{t}(\bm{x})+\log Z$ (which induces the same score as $\nabla_{\bm{x}}u^{*}_{t}(\bm{x})=\nabla_{\bm{x}}\log\pi_{t}(\bm{x})$).
In what follows, we describe how to find an approximation $u_{\theta}(\bm{x},t)$ to $u^{*}_{t}(\bm{x})$ within the PINN framework.

##### Target-informed parameterization.

To incorporate the initial condition ([11b](#S4.E11.2)), we use the following parameterization for the log-density function

$$ $u_{\theta}(\bm{x},t)=\frac{T-t}{T}\log\mu(\bm{x})+\frac{t}{T}\times\mathrm{NN}_{\theta}(\bm{x},t),\quad\forall(\bm{x},t)\in\mathbb{R}^{d}\times[0,T],$ (12) $$

where $\mathrm{NN}_{\theta}(\bm{x},t):\mathbb{R}^{d}\times[0,T]\to\mathbb{R}$ is a deep neural network.
This parameterization satisfies the initial condition ([11b](#S4.E11.2)), thus we only need to consider the PINN residual loss induced by ([11a](#S4.E11.1)).
Similar strategy is also used in consistency models (Song et al., ).

##### Underlying distribution for collocation points.

When training PINN, it is very important to collect proper collocation points $(\bm{x}_{t},t)\in\mathbb{R}^{d}\times[0,T]$ where $\bm{x}_{t}\sim\nu_{t}$. We expect samples from $\nu_{t}$ to cover the high-density domain of $\pi_{t}$ where PINN can provide a good approximation.
To achieve this, we first generate samples $\bm{x}_{0}\sim\nu_{0}$ by running a short chain of Langevin Monte Carlo (LMC) for $\pi$ so that $\nu_{0}$ covers the high density domain of $\pi$.
Given $\bm{x}_{0}\sim\nu_{0}$, we obtain $\bm{x}_{t}\sim\nu_{t}$ by sampling from the conditional distribution of the forward process given $\bm{x}_{0}$, namely, $\bm{x}_{t}|\bm{x}_{0}\sim\pi_{t|0}(\cdot|\bm{x}_{0})$.

##### Training objective.

One useful property of the forward process ([1](#S3.E1)) is that $\bm{x}_{T}\sim\pi_{T}\approx\pi_{\textrm{prior}}$ when $T$ is large.
In practice, we may use this property to further regularize the PINN residual loss, leading to the following training objective:

$$ $\displaystyle L_{\textrm{train}}(u_{\theta}):=$ $\displaystyle\ \mathbb{E}_{t\sim\mathcal{U}[0,T]}\mathbb{E}_{\bm{x}_{t}\sim\nu_{t}}\left[\beta^{2}(t)\cdot\left\|\partial_{t}u_{\theta}(\bm{x}_{t},t)-\mathcal{L}_{\textrm{L-FPE}}u_{\theta}(\bm{x}_{t},t)\right\|^{2}\right]$ (13) $\displaystyle\ \quad\quad+\lambda\cdot\mathbb{E}_{\bm{z}\sim\pi_{\textrm{prior}}}\left[\ell_{\textrm{reg}}(u_{\theta};T,\bm{z})\right],$ $$

where $\ell_{\textrm{reg}}(u_{\theta};T,\bm{z}):=\left\|\nabla_{\bm{z}}u_{\theta}(\bm{z},T)-\nabla_{\bm{z}}\log\pi_{\textrm{prior}}(\bm{z})\right\|^{2}$ denotes the regularization term, $\beta(t)$ is a weight function and $\lambda$ is a regularization coefficient.
We seek a good approximation $u_{\theta}(\bm{x},t)$ by minimizing ([13](#S4.E13)) via stochastic optimization methods where the stochastic gradient is computed by Monte Carlo estimation.
Our algorithm is summarized in Algorithm [1](#alg1).

Once $u_{\theta}(\bm{x},t)$ is learned, the induced score approximation is then substituted into the reverse process ([8](#S4.E8)), resulting in a new variant of diffusion-based sampling method that we call *Diffusion-PINN Sampler* (DPS).

Figure: Algorithm 1 : Solving Log-density FPE via PINN

##### Hutchinson’s trick for the gradient of the PINN residual.

Hutchinson’s trace estimator provides a stochastic method for estimating the trace of any square matrix and is commonly used in Laplacian estimation.
However, directly using Hutchinson’s trick here can result in biased gradient estimation.
To address this issue, we propose a novel variant of Hutchinson’s trick that allows unbiased gradient estimation.
Recall that the PINN residual can be decomposed as

$$ $\partial_{t}u_{\theta}-\mathcal{L}_{\textrm{L-FPE}}u_{\theta}:=\underbrace{\partial_{t}u_{\theta}+\bm{f}\cdot\nabla_{\bm{x}}u_{\theta}+\nabla\cdot\bm{f}-\frac{g^{2}(t)}{2}\|\nabla_{\bm{x}}u_{\theta}\|^{2}}_{:=\mathcal{L}_{\textrm{I}}u_{\theta}}-\frac{g^{2}(t)}{2}\Delta u_{\theta}.$ $$

Using this decomposition, the PINN residual loss $\|\partial_{t}u_{\theta}-\mathcal{L}_{\textrm{L-FPE}}u_{\theta}\|^{2}$ has the following gradient,

$$ $\displaystyle\nabla_{\theta}\left\|\partial_{t}u_{\theta}-\mathcal{L}_{\textrm{L-FPE}}u_{\theta}\right\|^{2}=2\left(\mathcal{L}_{\textrm{I}}u_{\theta}-\frac{g^{2}(t)}{2}\Delta u_{\theta}\right)\nabla_{\theta}\left(\mathcal{L}_{\textrm{I}}u_{\theta}-\frac{g^{2}(t)}{2}\Delta u_{\theta}\right)$ $\displaystyle=$ $\displaystyle 2\left(\mathcal{L}_{\textrm{I}}u_{\theta}-\frac{g^{2}(t)}{2}\cdot\mathbb{E}_{v_{1}}\left[v_{1}^{\top}\nabla_{\bm{x}}\left(v_{1}^{\top}\nabla_{\bm{x}}u_{\theta}\right)\right]\right)\nabla_{\theta}\left(\mathcal{L}_{\textrm{I}}u_{\theta}-\frac{g^{2}(t)}{2}\cdot\mathbb{E}_{v_{2}}\left[v_{2}^{\top}\nabla_{\bm{x}}\left(v_{2}^{\top}\nabla_{\bm{x}}u_{\theta}\right)\right]\right)$ $\displaystyle=$ $\displaystyle\mathbb{E}_{v_{1},v_{2}}\left[2\left(\mathcal{L}_{\textrm{I}}u_{\theta}-\frac{g^{2}(t)}{2}\cdot v_{1}^{\top}\nabla_{\bm{x}}\left(v_{1}^{\top}\nabla_{\bm{x}}u_{\theta}\right)\right)\nabla_{\theta}\left(\mathcal{L}_{\textrm{I}}u_{\theta}-\frac{g^{2}(t)}{2}\cdot v_{2}^{\top}\nabla_{\bm{x}}\left(v_{2}^{\top}\nabla_{\bm{x}}u_{\theta}\right)\right)\right],$ $$

where $v_{1}$ and $v_{2}$ are independent and satisfy $\mathbb{E}_{v_{1}}[v_{1}v_{1}^{\top}]=\mathbb{E}_{v_{2}}[v_{2}v_{2}^{\top}]=\bm{I}_{d}$.
Therefore, the following objective yields an unbiased gradient estimate of the PINN residual loss,

$$ $\mathbb{E}_{v_{1},v_{2}}\left[\textrm{Detach}\left(2\left(\mathcal{L}_{\textrm{I}}u_{\theta}-\frac{g^{2}(t)}{2}\cdot v_{1}^{\top}\nabla_{\bm{x}}\left(v_{1}^{\top}\nabla_{\bm{x}}u_{\theta}\right)\right)\right)\left(\mathcal{L}_{\textrm{I}}u_{\theta}-\frac{g^{2}(t)}{2}\cdot v_{2}^{\top}\nabla_{\bm{x}}\left(v_{2}^{\top}\nabla_{\bm{x}}u_{\theta}\right)\right)\right].$ $$

## 5 Theoretical Guarantees

##### Notations.

Let us denote $e_{t}(\bm{x}):=u_{\theta}(\bm{x},t)-u_{t}^{*}(\bm{x})$ and $r_{t}(\bm{x}):=\partial_{t}u_{\theta}(\bm{x},t)-\mathcal{L}_{\textrm{L-FPE}}u_{\theta}(\bm{x},t)$.
For any $C\in\mathbb{R}$ and $t\in[0,T]$, we define the weighted PINN objective on $\Omega$ as

$$ $L_{\textrm{PINN}}(t;C):=\int_{0}^{t}e^{C(t-s)}\|r_{s}(\cdot)\|^{2}_{L^{2}(\Omega;\nu_{s})}\mathrm{\leavevmode\nobreak\ d}s,$ (15) $$

where $\{\nu_{t}\}_{t=0}^{T}$ denotes the underlying distribution for collocation points introduced in Section [4.2](#S4.SS2) which satisfies the FPE $\partial_{t}\nu_{t}(\bm{x})=\mathcal{L}_{\textrm{FPE}}\nu_{t}(\bm{x})$.

### 5.1 Approximation Error of PINN for Log-density FPE

In this section, we provide an upper bound on the approximation error of PINN for the log-density FPE ([6](#S3.E6)) on a constrained domain $\Omega$.
Namely, we control $\|e_{t}(\cdot)\|^{2}_{L^{2}(\Omega;\nu_{t})}$ and $\|\nabla_{\bm{x}}e_{t}(\cdot)\|^{2}_{L^{2}(\Omega;\nu_{t})}$ by the residual loss $\|r_{t}(\cdot)\|^{2}_{L^{2}(\Omega;\nu_{t})}$ and the weighted PINN objective ([15](#S5.E15)).
We make the following assumptions.

###### Assumption 1 .

$u^{*}$ and $u_{\theta}$ are the same on the boundary, i.e., $u^{*}_{t}(\bm{x})=u_{\theta}(\bm{x},t)$ on $\partial\Omega\times[0,T]$.

###### Assumption 2 .

For any $t\in\left[0,T\right]$, $g^{2}(t)$ is bounded: $m_{1}\leqslant g^{2}(t)\leqslant M_{1}$ for some $m_{1},M_{1}>0$.

###### Assumption 3 .

$\log\nu_{t}(\bm{x}),u^{*}_{t}(\bm{x}),u_{\theta}(\bm{x},t)\in\mathcal{C}^{2}(\Omega\times[0,T])$.

Assumption [1](#Thmassumption1) is necessary for us to ensure the uniqueness of the solution to ([6](#S3.E6)) on $\Omega$, which is also considered in Deveney et al. (); Wang et al. ().
Assumption [2](#Thmassumption2), [3](#Thmassumption3) are also considered in Deveney et al. ().
Based on Assumption [3](#Thmassumption3), there exists $B^{\nu}_{0},B^{*}_{0},\widehat{B}_{0},B^{\nu}_{1},B^{*}_{1},\widehat{B}_{1}\in\mathbb{R}_{+}$ and $B^{\nu}_{2},B^{*}_{2},\widehat{B}_{2}\in\mathbb{R}$ depended on $\Omega$ such that for any $(\bm{x},t)\in\Omega\times[0,T]$, we have

$$ $\displaystyle|\partial_{t}\log\nu_{t}(\bm{x})|\leqslant B^{\nu}_{0},\quad|\partial_{t}u^{*}_{t}(\bm{x})|\leqslant B^{*}_{0},\quad|\partial_{t}u_{\theta}(\bm{x},t)|\leqslant\widehat{B}_{0},$ $\displaystyle\|\nabla_{\bm{x}}\log\nu_{t}(\bm{x})\|^{2}\leqslant B^{\nu}_{1},\quad\|\nabla_{\bm{x}}u^{*}_{t}(\bm{x})\|^{2}\leqslant B^{*}_{1},\quad\|\nabla_{\bm{x}}u_{\theta}(\bm{x},t)\|^{2}\leqslant\widehat{B}_{1},$ $\displaystyle\Delta\log\nu_{t}(\bm{x})\leqslant B^{\nu}_{2},\quad\Delta u^{*}_{t}(\bm{x})\geqslant B^{*}_{2},\quad\Delta u_{\theta}(\bm{x},t)\geqslant\widehat{B}_{2}.$ $$

In practice, using weights clipping strategy as in Arjovsky et al. (), we can control the regularity of neural network approximation $u_{\theta}(\bm{x},t)$, thus bound the constants $\widehat{B}_{0},\widehat{B}_{1},\widehat{B}_{2}$.

We summarize our main results in the following theorem.
The proof is deferred to Appendix [A.3](#A1.SS3), which generalizes the framework in Deveney et al. ().

###### Theorem 2 .

Suppose that Assumption [1](#Thmassumption1), [2](#Thmassumption2), and [3](#Thmassumption3) hold.
We further assume that $u_{\theta}(\bm{x},0)=u^{*}_{0}(\bm{x})$(^†^††This is a reasonable assumption due to the target-informed parameterization introduced in Section [4.2](#S4.SS2).) for any $\bm{x}\in\Omega$.
Then for any positive constant $\varepsilon>0$, the following holds for any $0\leqslant t\leqslant T$,

$$ $\|e_{t}(\cdot)\|^{2}_{L^{2}(\Omega;\nu_{t})}\leqslant\varepsilon L_{\emph{{PINN}}}(t;C_{1}(\varepsilon)).$ (16) $$

Moreover, for any $0\leqslant t\leqslant T$,

$$ $m_{1}\|\nabla_{\bm{x}}e_{t}(\cdot)\|^{2}_{L^{2}(\Omega;\nu_{t})}\leqslant\varepsilon\|r_{t}(\cdot)\|^{2}_{L^{2}(\Omega;\nu_{t})}+C_{3}(\varepsilon)L_{\emph{{PINN}}}(t;C_{1}(\varepsilon))+C_{2}\sqrt{\varepsilon L_{\emph{{PINN}}}(t;C_{1}(\varepsilon))}.$ (17) $$

where $C_{2}:=2\sqrt{2}(\widehat{B}_{0}^{2}+B_{0}^{*2})^{1/2}$, $C_{3}(\varepsilon):=\varepsilon(C_{1}(\varepsilon)+B^{\nu}_{0})$, and $C_{1}(\varepsilon)$ is a constant depended on $B^{\nu}_{1}$, $B^{*}_{1}$, $\widehat{B}_{1}$, $B^{\nu}_{2}$, $B^{*}_{2}$, $\widehat{B}_{2}$ and $m_{1}$, $M_{1}$.

###### Remark 1 .

The results of Wang et al. () show that the
$L^{2}$-error cannot be bounded by the PINN residual with universal constants independent of the approximate solution.
Therefore, some natural continuity assumptions (Assumption [3](#Thmassumption3)) about the approximate solution are necessary to control the $L^{2}$-error by the PINN residual.
It is noted that this continuity assumption can be satisfied by regularizing the neural network via weight clipping (Arjovsky et al., ), and would not sacrifice much approximation accuracy as the true solution is initialized as the log-density of the target and follows the diffusion process (e.g., the OU process) that would only become smoother as time evolves.
Moreover, our upper bound of $L^{2}$-error depends on continuity constants rather than an universal bound.
In this regard, our analysis aligns with the results of Wang et al. (), but with a more flexible bound based on some natural continuity assumption in the context of diffusion-based sampling.

###### Assumption 1 .

###### Assumption 2 .

###### Assumption 3 .

###### Theorem 2 .

###### Remark 1 .

### 5.2 Convergence of Diffusion-PINN Sampler

In this section, we present our convergence analysis of DPS based on Theorem [2](#Thmtheorem2) and the analysis of score-based generative modeling in Chen et al. ().
Following Chen et al. (, ), we focus on the forward process with $\bm{f}(\bm{x},t)=-\frac{1}{2}\bm{x}$ and $g(t)\equiv 1$, which is driven by

$$ $\mathrm{\leavevmode\nobreak\ d}\bm{x}_{t}=-\frac{1}{2}\bm{x}_{t}\mathrm{\leavevmode\nobreak\ d}t+\mathrm{\leavevmode\nobreak\ d}\bm{B}_{t},\quad\bm{x}_{0}\sim\pi,\quad 0\leqslant t\leqslant T.$ (18) $$

In practice, we use a discrete-time approximation for the reverse process.
Let $0=t_{0}<\cdots<t_{N}=T$ be the discretization points and $h_{k}:=t_{k}-t_{k-1}$ be the step size for $1\leqslant k\leqslant N$.
Let $t_{k}^{\prime}:=T-t_{N-k}$ for $0\leqslant k\leqslant N$ be the corresponding discretization points in the reverse SDE.
In our analysis, we consider the exponential integrator scheme which leads to the following sampling dynamics for $0\leqslant k\leqslant N-1$,

$$ $\mathrm{\leavevmode\nobreak\ d}\widehat{\bm{y}}_{t}=\left(\frac{1}{2}\widehat{\bm{y}}_{t}+\bm{s}_{T-t_{k}^{\prime}}(\widehat{\bm{y}}_{t_{k}^{\prime}})\right)\mathrm{\leavevmode\nobreak\ d}t+\mathrm{\leavevmode\nobreak\ d}\bm{B}_{t},\quad\widehat{\bm{y}}_{0}\sim\mathcal{N}(\bm{0},\bm{I}_{d}),\quad t\in[t_{k}^{\prime},t_{k+1}^{\prime}],$ (19) $$

where $\bm{s}_{t}(\bm{x})\approx\nabla_{\bm{x}}\log\pi_{t}(\bm{x})$ denotes the score approximation.
Let $\widehat{\pi}_{T}$ denotes the distribution of $\widehat{\bm{y}}_{T}$ from ([19](#S5.E19)).
We summarize all the assumptions we need as follows.

###### Assumption 4 .

The target distribution admits a density $\pi\in\mathcal{C}^{2}(\mathbb{R}^{d})$ where $\nabla_{\bm{x}}\log\pi(\bm{x})$ is $K$-Lipschitz and has the finite second moment, i.e., $M_{2}:=\mathbb{E}_{\pi}\left[\|\bm{x}\|^{2}\right]<\infty$.

###### Assumption 5 .

For any $\delta>0$, there exists bounded $\Omega$ such that $\int_{\Omega^{c}}\pi_{t}(\bm{x})\|\nabla_{\bm{x}}\log\pi_{t}(\bm{x})\|^{2}\mathrm{\leavevmode\nobreak\ d}\bm{x}\leqslant\delta$ for any $t\in[0,T]$.

###### Assumption 6 .

For any $(\bm{x},t)\in\Omega\times[0,T]$, there exists $R_{t}\geqslant 0$ depended on $t$, such that $\frac{\pi_{t}(\bm{x})}{\nu_{t}(\bm{x})}\leqslant R_{t}$.

Theorem [3](#Thmtheorem3) summarizes our main theoretical results of DPS.
The proof can be found in Appendix [A.4](#A1.SS4), which is based on the convergence results of score-based generative modeling in Chen et al. () and our Theorem [2](#Thmtheorem2).

###### Theorem 3 .

Suppose that $T\geqslant 1$, $K\geqslant 2$, and Assumptions [1](#Thmassumption1)-[6](#Thmassumption6) hold.
For any $\delta>0$, let $\Omega$ be chosen as in Assumption [5](#Thmassumption5).
For any positive constant $\varepsilon>0$, we further assume that $u_{\theta}(\bm{x},t)$ satisfies

$$ $\varepsilon\sum_{k=1}^{N}h_{k}R_{t_{k}}\|r_{t_{k}}(\cdot)\|^{2}_{L^{2}(\Omega;\nu_{t_{k}})}\leqslant\delta_{1}\quad{\textrm{and}}\quad\varepsilon\sum_{k=1}^{N}h_{k}R_{t_{k}}L_{\emph{{PINN}}}(t_{k};C_{1}(\varepsilon))\leqslant\delta_{2}.$ (20) $$

Then there is a universal constant $\alpha\geqslant 2$ such that, using step size
$h_{k}=h\min\{\max\{t_{k},1/(4K)\},1\}$, $0<h\leqslant 1/(\alpha d)$, and $\bm{s}_{t}(\bm{x})=\nabla_{\bm{x}}u_{\theta}(\bm{x},t)\cdot\mathbbm{1}\{\bm{x}\in\Omega\}$ in ([19](#S5.E19)), we have the following upper bound on the KL divergence between the target and the approximate distribution

$$ $\operatorname{KL}(\pi\|\widehat{\pi}_{T})\lesssim(d+M_{2})\cdot e^{-T}+d^{2}h(\log K+T)+T\delta+\delta_{1}+C_{5}(\varepsilon)\delta_{2}+C_{2}\sqrt{\sum_{k=1}^{N}h_{k}R_{t_{k}}\delta_{2}},$ (21) $$

where $C_{5}(\varepsilon):=C_{1}(\varepsilon)+B^{\nu}_{0}$, $C_{2}$ and $C_{1}(\varepsilon)$ are defined in Theorem [2](#Thmtheorem2).

###### Assumption 4 .

###### Assumption 5 .

###### Assumption 6 .

###### Theorem 3 .

### 5.3 Theoretical Comparison between Different Sampling Methods for Collocation Points

In practice, we typically lack prior knowledge of the high-probability regions of the diffusion path starting from the target distribution. As a result, specifying a sufficiently large support for uniform sampling of collocation points, becomes challenging and inefficient, especially in high-dimensional settings.
In contrast, we employ a more sophisticated strategy for generating collocation points that integrates Langevin Monte Carlo (LMC) with the forward pass (see Section [4.2](#S4.SS2) for details).
In this section, we theoretically demonstrate our strategy for collocation points generation has a faster convergence for DPS.

Similar to Theorem [2](#Thmtheorem2) and [3](#Thmtheorem3), we can establish theoretical guarantees for uniformly sampled collocation points.
However, these guarantees differ from those in Theorem [2](#Thmtheorem2) and [3](#Thmtheorem3), since the distribution of uniformly sampled collocation points does not depend on $t$ and does not satisfy the corresponding FPE.
This distinction necessitates a slightly different analysis, requiring an additional assumption compared to the previous analyses (see more details in Section [5.3.1](#S5.SS3.SSS1)).
Our findings are summarized in Theorem [4](#Thmtheorem4) and [5](#Thmtheorem5), with the corresponding proofs provided in Appendix [A.5](#A1.SS5) and Appendix [A.6](#A1.SS6).

#### 5.3.1 Approximation Error of PINN for Log-density FPE with Uniformly Sampled Collocation Points

In this section, we control the approximation error of PINN for the log-density FPE ([6](#S3.E6)) on a constrained domain $\Omega$ with uniformly sampled collocation points $\nu_{t}\sim{\rm Unif}(\Omega)$.
First, we redefine the weighted PINN objective on $\Omega$ as follows and state an additional assumption for our analysis in Assumption [7](#Thmassumption7),

$$ $L_{\textrm{PINN}}^{\textrm{Unif}}(t;C):=\int_{0}^{t}e^{C(t-s)}\left\|r_{s}(\cdot)\right\|^{2}_{L^{2}(\Omega)}{\textnormal{d}}s.$ $$

###### Assumption 7 .

For any $(\bm{x},t)\in\Omega\times[0,T]$, $\nabla\cdot\bm{f}(\bm{x},t)\leqslant\mathrm{m}_{2}$ for some $\mathrm{m}_{2}\in\mathbb{R}$.

###### Theorem 4 .

Suppose that Assumption [1](#Thmassumption1), [2](#Thmassumption2), [3](#Thmassumption3) and [7](#Thmassumption7) hold.
We further assume that $u_{\theta}(\bm{x},0)=u^{*}_{0}(\bm{x})$ for any $\bm{x}\in\Omega$.
Then for any positive constant ${\varepsilon}>0$, the following holds for any $t\in[0,T]$,

$$ $\|e_{t}(\cdot)\|^{2}_{L^{2}(\Omega)}\leqslant{\varepsilon}L_{\emph{{PINN}}}^{\emph{{Unif}}}(t;C_{1}^{\emph{{U}}}({\varepsilon})).$ (22) $$

Moreover, for any $t\in[0,T]$,

$$ $m_{1}\|\nabla_{\bm{x}}e_{t}(\cdot)\|^{2}_{L^{2}(\Omega)}\leqslant{\varepsilon}\|r_{t}(\cdot)\|^{2}_{L^{2}(\Omega)}+{\varepsilon}\cdot C_{1}^{\emph{{U}}}({\varepsilon})L_{\emph{{PINN}}}^{\emph{{Unif}}}(t;C_{1}^{\emph{{U}}}({\varepsilon}))+C_{2}\sqrt{{\varepsilon}L_{\emph{{PINN}}}^{\emph{{Unif}}}(t;C_{1}^{\emph{{U}}}({\varepsilon}))},$ (23) $$

where $C_{1}^{\emph{{U}}}({\varepsilon}):=\frac{1}{{\varepsilon}}+\mathrm{m}_{2}-\frac{m_{1}}{2}\left(B_{2}^{*}+\widehat{B}_{2}\right)$ and $C_{2}$ is defined in Theorem [2](#Thmtheorem2).

###### Remark 2 .

Indeed, we only need the assumption of the lower bound on $g^{2}(t)$ from Assumption [2](#Thmassumption2) and the continuity assumptions for $u_{t}^{*}(\bm{x})$ and $u_{\theta}(\bm{x},t)$ from Assumption [3](#Thmassumption3) in Theorem [4](#Thmtheorem4).

###### Assumption 7 .

###### Theorem 4 .

###### Remark 2 .

#### 5.3.2 Convergence of Diffusion-PINN Sampler with Uniformly Sampled Collocation Points

In this section, within the similar setting in section [5.2](#S5.SS2), we present the convergence analysis of DPS when the collocation points are uniformly sampled from $\Omega$.

###### Theorem 5 .

Suppose that $T\geqslant 1,K\geqslant 2$, and Assumption [1](#Thmassumption1), [2](#Thmassumption2), [3](#Thmassumption3), [4](#Thmassumption4), [5](#Thmassumption5), and [7](#Thmassumption7) hold.
For any $\delta>0$, let $\Omega$ be chosen as in Assumption [5](#Thmassumption5).
For any positive constant ${\varepsilon}>0$, we further assume that $u_{\theta}(\bm{x},t)$ satisfies the following(^‡^‡‡Here, we contain the term $\textrm{Vol}(\Omega)$ since the PINN residual objective used for uniform collocation points is given by $\|r_{t}(\cdot)\|^{2}_{L^{2}(\Omega)}/\textrm{Vol}(\Omega)$.),

$$ $\displaystyle\ {\varepsilon}\sum_{k=1}^{N}h_{k}\max_{\bm{x}\in\Omega}\left\{\pi_{t_{k}}(\bm{x})\right\}\cdot\left\|r_{t_{k}}(\cdot)\right\|^{2}_{L^{2}(\Omega)}\leqslant\delta_{1}\cdot\emph{{Vol}}(\Omega),$ (24) $\displaystyle\ {\varepsilon}\sum_{k=1}^{N}h_{k}\max_{\bm{x}\in\Omega}\left\{\pi_{t_{k}}(\bm{x})\right\}\cdot L_{\emph{{PINN}}}^{\emph{{Unif}}}(t_{k};C_{1}^{\emph{{U}}}({\varepsilon}))\leqslant\delta_{2}\cdot\emph{{Vol}}(\Omega).$ $$

Then there is a universal constant $\alpha\geqslant 2$ such that, using step size $h_{k}:=h\min\{\max\{t_{k},\frac{1}{4K}\}\}$ for $0<h\leqslant\frac{1}{\alpha d}$, and $\bm{s}_{t}(\bm{x})=\nabla_{\bm{x}}u_{\theta}(\bm{x},t)\cdot\bm{1}\{\bm{x}\in\Omega\}$,
we have the following upper bound on the KL divergence between the target and the approximate distribution,

$$ $\displaystyle\operatorname{KL}\left(\pi\|\widehat{\pi}_{T}\right)\lesssim$ $\displaystyle\ \left(d+M_{2}\right)\cdot e^{-T}+d^{2}h\left(\log K+T\right)+T\delta+\left(\delta_{1}+C_{1}^{\emph{{U}}}({\varepsilon})\delta_{2}\right)\cdot\emph{{Vol}}(\Omega)$ $\displaystyle\ \quad+C_{2}\sqrt{\sum_{k=1}^{N}h_{k}\max_{\bm{x}\in\Omega}\left\{\pi_{t_{k}}(\bm{x})\right\}\cdot\delta_{2}\cdot\emph{{Vol}}(\Omega)}.$ $$

where $C_{1}^{\emph{{U}}}({\varepsilon})$ and $C_{2}$ are defined in Theorem [4](#Thmtheorem4).

Specifically, our results indicate that employing LMC and the forward pass for sampling collocation points is advantageous over uniform sampling. This is because, in the uniform case, the KL bound includes a factor proportional to the volume of the support, ${\rm Vol}(\Omega)$, which can be prohibitively large in high dimensions. In contrast, our bound depends on the density ratio $\pi_{t}/\nu_{t}$, which is more manageable due to LMC and converges to $1$ as $t$ increases, thanks to the forward process.

###### Theorem 5 .

## 6 Numerical Experiments

In this section, we conduct experiments on various sampling tasks to demonstrate the effectiveness and efficiency of the Diffusion-PINN Sampler (DPS) compared to previous methods.
Our sampling tasks includes 9-Gaussians ($d=2$), Rings ($d=2$), Funnel (Neal, ) ($d=10$), and Double-well ($d=30$), which are commonly used to evaluate diffusion-based sampling algorithms (Zhang & Chen, ; Berner et al., ; Grenioux et al., ).
For multimodal distributions, the modes are designed to be well-separated, with challenging mixing proportions between different modes (see more details in Appendix [C.2](#A3.SS2)).
For DPS, we employ a time-rescaled forward process and use a weight function $\beta(t)=2(1-t)$ for the PINN residual loss to improve numerical stability.
To generate collocation points for each task, we run a short chain of LMC with a relatively large step size for better coverage of the high-density domain.
For 9-Gaussians, Rings, and Double-well, the PINN residual loss alone suffices for good performance, so we set the regularization coefficient $\lambda=0$.
For Funnel, however, regularization proves helpful, and we set $\lambda=1$ (details in Section [6.3](#S6.SS3)).
More details on experiment settings can be found in Appendix [C.3](#A3.SS3).

Figure: Figure 2: Comparison between solving log-density FPE by PINN and denoising score matching on score estimation.
Refer to caption: /html/2410.15336/assets/x2.png

##### Baselines.

We benchmark DPS performance against a wide range of strong baseline methods.
For MCMC methods, we consider the Langevin Monte Carlo (LMC).
As for sampling methods using reverse diffusion, we include RDMC (Huang et al., ) and SLIPS (Grenioux et al., ).
We also compare with the VI-based PIS (Zhang & Chen, ) and DIS (Berner et al., ).
See Appendix [C.1](#A3.SS1) for more details.

### 6.1 Score Estimation

We first evaluate the accuracy of score function estimates obtained by solving the log-density FPE (Algorithm [1](#alg1)).
To do that, we conduct an experiment on the 9-Gaussians target $\pi$ where we know the ground truth scores throughout the entire forward process.
Figure [2](#S6.F2) shows the $L^{2}(\pi)$ error of the score estimation for our method compared to denoising score matching (Vincent, ; Song et al., ).
We see clearly that our method provides more accurate score estimation than denoising score matching.

Figure: Figure 3: Sampling performance of different methods for 9-Gaussians ($d=2$), Rings ($d=2$), Funnel ($d=10$), and Double-well ($d=30$).
Refer to caption: /html/2410.15336/assets/x3.png

**Table 1: KL divergence ($\downarrow$) to the ground truth obtained by different methods. Bold font indicates the best results. We use the KL divergence of the first two dimensions for Funnel ($d=10$) and the KL divergence of the first five dimensions for Double-well ($d=30$). All the KL divergence is computed by the ITE package (Szabó, ).**
| Target | LMC | RDMC | SLIPS | PIS | DIS | DPS (ours) |
| --- | --- | --- | --- | --- | --- | --- |
| 9-Gaussians | $1.6568_{\pm 0.0189}$ | $1.0844_{\pm 0.0132}$ | $0.0901_{\pm 0.0071}$ | $2.0042_{\pm 0.0203}$ | $2.2758_{\pm 0.0240}$ | $\mathbf{0.0131_{\pm 0.0093}}$ |
| Rings | $2.4754_{\pm 0.0302}$ | $0.7487_{\pm 0.0073}$ | $0.4127_{\pm 0.0144}$ | $2.6985_{\pm 0.0290}$ | $2.3433_{\pm 0.0275}$ | $\mathbf{0.0176_{\pm 0.0059}}$ |
| Funnel | $0.1908_{\pm 0.0156}$ | $2.0250_{\pm 0.0364}$ | $0.1971_{\pm 0.0133}$ | $0.4377_{\pm 0.0199}$ | $0.2383_{\pm 0.0169}$ | $\mathbf{0.0846_{\pm 0.0122}}$ |
| Double-well | $0.1915_{\pm 0.0122}$ | $1.5735_{\pm 0.0162}$ | $0.4840_{\pm 0.0145}$ | $0.0969_{\pm 0.0114}$ | $0.6796_{\pm 0.0139}$ | $\mathbf{0.0273_{\pm 0.0113}}$ |

**Table 2: $L^{2}$ error ($\downarrow$) of the mixing proportions estimation when sampling multimodal target distributions using different methods. Bold font indicates the best results. All the estimation is computed with 1,000 samples.**
| Target | LMC | RDMC | SLIPS | PIS | DIS | DPS (ours) |
| --- | --- | --- | --- | --- | --- | --- |
| 9-Gaussians | $0.5199_{\pm 0.0159}$ | $0.1313_{\pm 0.0099}$ | $0.0018_{\pm 0.0005}$ | $0.4893_{\pm 0.0110}$ | $0.7268_{\pm 0.0146}$ | $\mathbf{0.0006_{\pm 0.0003}}$ |
| Rings | $0.6005_{\pm 0.0251}$ | $0.0537_{\pm 0.0035}$ | $0.2471_{\pm 0.0144}$ | $0.8016_{\pm 0.0194}$ | $0.5233_{\pm 0.0194}$ | $\mathbf{0.0006_{\pm 0.0006}}$ |
| Double-well | $0.0673_{\pm 0.0082}$ | $0.2154_{\pm 0.0075}$ | $0.1645_{\pm 0.0113}$ | $0.0044_{\pm 0.0011}$ | $0.0684_{\pm 0.0035}$ | $\mathbf{0.0004_{\pm 0.0002}}$ |

### 6.2 Sample Quality

In this section, we compare DPS with the aforementioned baseline methods on various target distributions.
We use KL divergence to evaluate the quality of samples provided by different methods in low dimensional problems (9-Gaussians, Rings), and use the projected KL divergence instead for Funnel and Double-well that are problems with relatively higher dimensions.
The results are reported in Table [1](#S6.T1).
Figure [3](#S6.F3) visualizes the samples from different methods.
We clearly see that DPS provides the best approximation accuracy and sample quality among all methods.
Although we use LMC to generate collocation points, DPS greatly outperforms LMC, indicating the power of diffusion-based sampling methods with learned score functions.

For multimodal distributions, we estimate the mixing proportions for different modes using samples generated by different methods, and evaluate the estimation accuracy in terms of $L^{2}$ error to the true weights.
The results are shown in Table [2](#S6.T2).
It is clear that DPS provides accurate weights estimation while other baselines tend to struggle to learn the weights.

### 6.3 Ablation Study

In this section, we compare the performance of score estimation between solving the score FPE and the log-density FPE, and investigate the effect of regularization in DPS.

We first solve the corresponding score FPE and log-density FPE for a MoG with two distant modes: $\pi^{M}=0.2\mathcal{N}((-5,-5)^{\prime},\bm{I}_{2})+0.8\mathcal{N}((5,5)^{\prime},\bm{I}_{2})$.
The left plot in Figure [4](#S6.F4) show the PINN residual loss and the score estimation error as functions of the number of iterations.
We see that for the score FPE, the score approximation error decreases rapidly at first but quickly levels off, while the PINN residual loss continues to decrease with more iterations.
In contrast, when solving the log-density FPE, the PINN residual loss and the score approximation error decrease consistently, resulting in more accurate score approximation overall.
The middle and right plots in Figure [4](#S6.F4) display the histogram based on samples generated from the reverse SDE using the score estimates from both methods, together with the true marginal density. We observe that the score FPE-based method fails to identify the correct mixing proportions, whereas the log-density FPE-based method successfully recovers the correct weights.

Next, we solve the log-density PFE with different regularization coefficients $\lambda$ on the Funnel target.
Figure [5](#S6.F5) (left) shows the KL divergence for various $\lambda$ as a function of the number of iterations.
We see that, compared to the non-regularized case ($\lambda=0$), both the convergence speed and overall approximation accuracy have been greatly improved when regularization is applied.
The middle and right plots in Figure [5](#S6.F5) show the samples generated from DPS with $\lambda=1$ and $\lambda=0$ respectively.
With regularization, DPS provides a better fit to the target distribution, more accurately capturing the thickness in the tails. This indicates that regularization could be beneficial for heavy-tail distributions.

Figure: Figure 4: Left: PINN residual loss and score approximation error during solving score/log-density FPE; Middle/Right: Marginals of the first dimension from DPS by solving score/log-density FPE for MoG with two modes.
Refer to caption: /html/2410.15336/assets/x4.png

Figure: Figure 5: Left: KL divergence to the ground truth during solving log-density FPE with different regularization for Funnel. Middle/Right: Sampling performance of DPS with/without regularization for Funnel.
Refer to caption: /html/2410.15336/assets/x5.png

## 7 Conclusion

In this work, we proposed *Diffusion-PINN Sampler* (DPS), a novel method that leverages Physics-Informed Neural Networks (PINN) and diffusion models for accurate sampling from complex target distributions.
By solving the log-density FPE that governs the evolution of the log-density of the underlying SDE marginals via PINN, DPS demonstrates accurate sampling capabilities even for distributions with multiple modes or heavy tails, and it excels in identifying mixing proportions when the target features isolated modes.
The control of log-density estimation error via PINN residual loss ensures convergence guarantees to the target distribution, building upon established results for score-based diffusion models.
We demonstrated the effectiveness of our approach on multiple numerical examples.
Limitations are discussed in Appendix [B](#A2).

## Acknowledgements

This work was supported by National Natural Science Foundation of China (grant no. 12201014 and grant no. 12292983).
The research of Cheng Zhang was support in part by National Engineering Laboratory for Big Data Analysis and Applications, the Key Laboratory of Mathematics and Its Applications (LMAM) and the Key Laboratory of Mathematical Economics and Quantitative Finance (LMEQF) of Peking University.
Zhekun Shi is partially supported by the elite undergraduate training program of School of Mathematical Sciences in Peking University.

## Appendix A Proofs

### A.1 Proof of Theorem 1

###### Proof of Theorem 1 .

Recall that $p_{t}(\bm{x})$ denotes the marginal density of $\bm{x}_{t}$ following the forward process ([1](#S3.E1)),
and satisfies

$$ $\partial_{t}p_{t}(\bm{x})=\frac{1}{2}g^{2}(t)\Delta p_{t}(\bm{x})-\nabla\cdot\left[\bm{f}(\bm{x},t)p_{t}(\bm{x})\right].$ (25) $$

Therefore, the log-density $u_{t}(\bm{x}):=\log p_{t}(\bm{x})$ satisfies

$$ $\partial_{t}u_{t}(\bm{x})=\frac{\partial_{t}p_{t}(\bm{x})}{p_{t}(\bm{x})}=\frac{1}{2}g^{2}(t)\frac{\Delta p_{t}(\bm{x})}{p_{t}(\bm{x})}-\frac{\nabla\cdot\left[\bm{f}(\bm{x},t)p_{t}(\bm{x})\right]}{p_{t}(\bm{x})}.$ (26) $$

Note that we have the identities

$$ $\displaystyle\Delta p_{t}(\bm{x})=\nabla\cdot\left[p_{t}(\bm{x})\nabla_{\bm{x}}u_{t}(\bm{x})\right]=\nabla_{\bm{x}}p_{t}(\bm{x})\cdot\nabla_{\bm{x}}u_{t}(\bm{x})+p_{t}(\bm{x})\Delta u_{t}(\bm{x}),$ (27) $\displaystyle\nabla\cdot\left[\bm{f}(\bm{x},t)p_{t}(\bm{x})\right]=\nabla_{\bm{x}}p_{t}(\bm{x})\cdot\bm{f}(\bm{x},t)+p_{t}(\bm{x})\left[\nabla\cdot\bm{f}(\bm{x},t)\right].$ $$

Plug ([27](#A1.E27)) into ([26](#A1.E26)), we have

$$ $\partial_{t}u_{t}(\bm{x})=\frac{1}{2}g^{2}(t)\Delta u_{t}(\bm{x})+\frac{1}{2}g^{2}(t)\left\|\nabla_{\bm{x}}u_{t}(\bm{x})\right\|^{2}-\bm{f}(\bm{x},t)\cdot\nabla_{\bm{x}}u_{t}(\bm{x})-\nabla\cdot\bm{f}(\bm{x},t).$ $$

Since $\log p_{t}(\bm{x})$ is sufficiently smooth, we can swap the order of differentiation and get

$$ $\partial_{t}\bm{s}_{t}(\bm{x})=\partial_{t}\nabla_{\bm{x}}u_{t}(\bm{x})=\nabla_{\bm{x}}\partial_{t}u_{t}(\bm{x}).$ $$

Hence, the theorem is proved.
∎

###### Proof of Theorem 1 .

### A.2 Omitted Proof in Example 1

##### Notations.

For two probability measures $\nu_{1}$ and $\nu_{2}$ in $\mathbb{R}^{d}$, we define the $L^{2}(p)$ error of their scores as $\mathrm{SE}_{p}(\nu_{1}\|\nu_{2}):=\mathbb{E}_{\bm{x}\sim p}[\|\nabla_{\bm{x}}\log\nu_{1}(\bm{x})-\nabla_{\bm{x}}\log\nu_{2}(\bm{x})\|^{2}]$ where $p$ also denotes a probability measure.
Note that if we choose $p=\nu_{1}$, we have $\mathrm{SE}_{\nu_{1}}(\nu_{1}\|\nu_{2})=F(\nu_{1},\nu_{2})$ where $F(\nu_{1},\nu_{2})$ denotes the Fisher divergence between $\nu_{1}$ and $\nu_{2}$.
For any $\bm{a}\in\mathbb{R}^{d}$, we denote $\gamma_{\bm{a}}(\bm{x}):=\exp(-\|\bm{x}-\bm{a}\|^{2}/2)$.
For simplify, we denote $\mathbb{E}_{\bm{x}\sim\mathcal{N}(\bm{a},I_{d})}[\cdot]$ by $\mathbb{E}_{\gamma_{\bm{a}}}[\cdot]$.
Thus the probability density of $\mathcal{N}(\bm{a},I_{d})$ is $p(\bm{x})=\gamma_{\bm{a}}(\bm{x})/(\sqrt{2\pi})^{d}$.
For the MoG
$\pi^{M}=w_{1}\mathcal{N}(\bm{a}_{1},I_{d})+w_{2}\mathcal{N}(\bm{a}_{2},I_{d})$, the score is given by

$$ $\nabla_{\bm{x}}\log\pi^{M}(\bm{x})=\frac{w_{1}\bm{a}_{1}\gamma_{\bm{a}_{1}}(\bm{x})+w_{2}\bm{a}_{2}\gamma_{\bm{a}_{2}}(\bm{x})}{w_{1}\gamma_{\bm{a}_{1}}(\bm{x})+w_{2}\gamma_{\bm{a}_{2}}(\bm{x})}-\bm{x}.$ (28) $$

Then we show our general results in Theorem [6](#Thmtheorem6) where we state a lower bound of $\mathrm{KL}(\pi^{M}\|\tilde{\pi}^{M})$ and an upper bound of $\mathrm{SE}_{p}(\pi^{M}\|\tilde{\pi}^{M})$.

###### Theorem 6 .

Consider two MoGs in $\mathbb{R}^{d}$: $\pi^{M}=w_{1}\mathcal{N}(\bm{a}_{1},I_{d})+w_{2}\mathcal{N}(\bm{a}_{2},I_{d})$, $\tilde{\pi}^{M}=\tilde{w}_{1}\mathcal{N}(\bm{a}_{1},I_{d})+\tilde{w}_{2}\mathcal{N}(\bm{a}_{2},I_{d})$ where $\bm{a}_{1},\bm{a}_{2}\in\mathbb{R}^{d}$, $w_{1},w_{2},\tilde{w}_{1},\tilde{w}_{2}>0$, $w_{1}+w_{2}=1$ and $\tilde{w}_{1}+\tilde{w}_{2}=1$. Then $\operatorname{KL}\left(\pi^{M}\|\tilde{\pi}^{M}\right)$ is lower bounded by

$$ $\displaystyle\operatorname{KL}\left(\pi^{M}\|\tilde{\pi}^{M}\right)\geqslant$ $\displaystyle\ w_{1}\left(\log w_{1}-\log\left(\tilde{w}_{1}+\exp\left(-\frac{\left\|\bm{a}_{1}-\bm{a}_{2}\right\|^{2}}{4}\right)\right)\right)$ (29) $\displaystyle\ +w_{2}\left(\log w_{2}-\log\left(\tilde{w}_{2}+\exp\left(-\frac{\left\|\bm{a}_{1}-\bm{a}_{2}\right\|^{2}}{4}\right)\right)\right)$ $\displaystyle\ -\left(\log 4+d\right)\exp\left(\frac{d}{2}\log 2-\frac{\left\|\bm{a}_{1}-\bm{a}_{2}\right\|^{2}}{64}\right).$ $$

Let $p(\bm{x})$ denotes any distribution that is absolutely continuous w.r.t. $\mu$, then $\mathrm{SE}_{p}(\pi^{M}\|\tilde{\pi}^{M})$ is upper bounded by

$$ $\displaystyle\mathrm{SE}_{p}(\pi^{M}\|\tilde{\pi}^{M})\leqslant$ $\displaystyle\ 2\exp\left(-\frac{\left\|\bm{a}_{1}-\bm{a}_{2}\right\|^{2}}{2}\right)\left[\frac{w_{2}^{2}}{w_{1}^{2}}+\frac{\tilde{w}_{2}^{2}}{\tilde{w}_{1}^{2}}+\frac{w_{1}^{2}}{w_{2}^{2}}+\frac{\tilde{w}_{1}^{2}}{\tilde{w}_{2}^{2}}\right]\left\|\bm{a}_{1}-\bm{a}_{2}\right\|^{2}$ (30) $\displaystyle\ +8\left[\left\|\bm{a}_{1}\right\|^{2}+\left\|\bm{a}_{2}\right\|^{2}\right]\int_{\Omega_{3}}p(\bm{x})\mathrm{\leavevmode\nobreak\ d}\bm{x},$ $$

where $\Omega_{1}=\left\{\bm{x}\in\mathbb{R}^{d}:\left\|\bm{x}-\bm{a}_{1}\right\|\leqslant\frac{\left\|\bm{a}_{1}-\bm{a}_{2}\right\|}{4}\right\}$, $\Omega_{2}=\left\{\bm{x}\in\mathbb{R}^{d}:\left\|\bm{x}-\bm{a}_{2}\right\|\leqslant\frac{\left\|\bm{a}_{1}-\bm{a}_{2}\right\|}{4}\right\}$, and $\Omega_{3}=\Omega_{1}^{c}\bigcap\Omega_{2}^{c}$.

###### Remark 3 .

If we choose $p(\bm{x})=\pi^{M}(\bm{x})$ in Theorem [6](#Thmtheorem6), the Fisher divergence $F(\pi^{M},\tilde{\pi}^{M})$ is upper bounded by

$$ $\displaystyle F(\pi^{M},\tilde{\pi}^{M})\leqslant$ $\displaystyle\ 2\exp\left(-\frac{\left\|\bm{a}_{1}-\bm{a}_{2}\right\|^{2}}{2}\right)\left[\frac{w_{2}^{2}}{w_{1}^{2}}+\frac{\tilde{w}_{2}^{2}}{\tilde{w}_{1}^{2}}+\frac{w_{1}^{2}}{w_{2}^{2}}+\frac{\tilde{w}_{1}^{2}}{\tilde{w}_{2}^{2}}\right]\left\|\bm{a}_{1}-\bm{a}_{2}\right\|^{2}$ (31) $\displaystyle\ +8\left[\left\|\bm{a}_{1}\right\|^{2}+\left\|\bm{a}_{2}\right\|^{2}\right]\exp\left(\frac{d}{2}\log 2-\frac{\left\|\bm{a}_{1}-\bm{a}_{2}\right\|^{2}}{64}\right),$ $$

where we use the following inequality

$$ $\displaystyle\int_{\Omega_{3}}\pi^{M}(\bm{x})\mathrm{\leavevmode\nobreak\ d}\bm{x}=$ $\displaystyle w_{1}\int_{\Omega_{3}}\frac{1}{(\sqrt{2\pi})^{d}}\gamma_{\bm{a}_{1}}(\bm{x})\mathrm{\leavevmode\nobreak\ d}\bm{x}+w_{2}\int_{\Omega_{3}}\frac{1}{(\sqrt{2\pi})^{d}}\gamma_{\bm{a}_{2}}(\bm{x})\mathrm{\leavevmode\nobreak\ d}\bm{x}$ $\displaystyle\leqslant$ $\displaystyle\ w_{1}\exp\left(\frac{d}{2}\log 2-\frac{\left\|\bm{a}_{1}-\bm{a}_{2}\right\|^{2}}{64}\right)+w_{2}\exp\left(\frac{d}{2}\log 2-\frac{\left\|\bm{a}_{1}-\bm{a}_{2}\right\|^{2}}{64}\right)$ $\displaystyle=$ $\displaystyle\exp\left(\frac{d}{2}\log 2-\frac{\left\|\bm{a}_{1}-\bm{a}_{2}\right\|^{2}}{64}\right).$ $$

Thus Example [1](#Thmexample1) holds naturally.

###### Proof of Theorem 6 .

We first prove ([29](#A1.E29)).
We can decompose $\operatorname{KL}\left(\pi^{M}\|\tilde{\pi}^{M}\right)$ as

$$ $\displaystyle\operatorname{KL}\left(\pi^{M}\|\tilde{\pi}^{M}\right)=\mathbb{E}_{\pi^{M}}\left[\log\left(\frac{\pi^{M}(\bm{x})}{\tilde{\pi}^{M}(\bm{x})}\right)\right]$ (32) $\displaystyle=$ $\displaystyle\ w_{1}\mathbb{E}_{\gamma_{\bm{a}_{1}}}\left[\log\left(\frac{w_{1}\gamma_{\bm{a}_{1}}(\bm{x})+w_{2}\gamma_{\bm{a}_{2}}(\bm{x})}{\tilde{w}_{1}\gamma_{\bm{a}_{1}}(\bm{x})+\tilde{w}_{2}\gamma_{\bm{a}_{2}}(\bm{x})}\right)\right]+w_{2}\mathbb{E}_{\gamma_{\bm{a}_{2}}}\left[\log\left(\frac{w_{1}\gamma_{\bm{a}_{1}}(\bm{x})+w_{2}\gamma_{\bm{a}_{2}}(\bm{x})}{\tilde{w}_{1}\gamma_{\bm{a}_{1}}(\bm{x})+\tilde{w}_{2}\gamma_{\bm{a}_{2}}(\bm{x})}\right)\right].$ $$

Note that

$$ $\displaystyle\mathbb{E}_{\gamma_{\bm{a}_{1}}}\left[\log\left(\frac{w_{1}\gamma_{\bm{a}_{1}}(\bm{x})+w_{2}\gamma_{\bm{a}_{2}}(\bm{x})}{\tilde{w}_{1}\gamma_{\bm{a}_{1}}(\bm{x})+\tilde{w}_{2}\gamma_{\bm{a}_{2}}(\bm{x})}\right)\right]=$ $\displaystyle\mathbb{E}_{\gamma_{\bm{a}_{1}}}\left[\log\left(\frac{w_{1}+w_{2}\gamma_{\bm{a}_{2}}(\bm{x})/\gamma_{\bm{a}_{1}}(\bm{x})}{\tilde{w}_{1}+\tilde{w}_{2}\gamma_{\bm{a}_{2}}(\bm{x})/\gamma_{\bm{a}_{1}}(\bm{x})}\right)\right]$ (33) $\displaystyle\geqslant$ $\displaystyle\log w_{1}-\mathbb{E}_{\gamma_{\bm{a}_{1}}}\left[\log\left(\tilde{w}_{1}+\tilde{w}_{2}\frac{\gamma_{\bm{a}_{2}}(\bm{x})}{\gamma_{\bm{a}_{1}}(\bm{x})}\right)\right].$ $$

Let $\widetilde{\Omega}_{1}=\{\bm{x}\in\mathbb{R}^{d}:\|\bm{x}-\bm{a}_{1}\|\leqslant\frac{\left\|\bm{a}_{1}-\bm{a}_{2}\right\|}{4}\}$, $\widetilde{\Omega}_{2}=\widetilde{\Omega}_{1}^{c}\bigcap\{\bm{x}\in\mathbb{R}^{d}:\tilde{w}_{2}\gamma_{\bm{a}_{2}}(\bm{x})/\gamma_{\bm{a}_{1}}(\bm{x})\leqslant\tilde{w}_{1}\}$, and $\widetilde{\Omega}_{3}=(\widetilde{\Omega}_{1}\bigcup\widetilde{\Omega}_{2})^{c}=\widetilde{\Omega}_{1}^{c}\bigcap\widetilde{\Omega}_{2}^{c}$.
Then for any $\bm{x}\in\widetilde{\Omega}_{1}$, we have $\left\|\bm{x}-\bm{a}_{2}\right\|\geqslant\left\|\bm{a}_{1}-\bm{a}_{2}\right\|-\left\|\bm{x}-\bm{a}_{1}\right\|\geqslant\frac{3}{4}\left\|\bm{a}_{1}-\bm{a}_{2}\right\|$, thus $\gamma_{\bm{a}_{2}}(\bm{x})/\gamma_{\bm{a}_{1}}(\bm{x})=\exp\left(\frac{\left\|\bm{x}-\bm{a}_{1}\right\|^{2}-\left\|\bm{x}-\bm{a}_{2}\right\|^{2}}{2}\right)\leqslant\exp\left(-\frac{\left\|\bm{a}_{1}-\bm{a}_{2}\right\|^{2}}{4}\right)$.
Then we have

$$ $\displaystyle\int_{\widetilde{\Omega}_{1}}\frac{1}{\left(\sqrt{2\pi}\right)^{d}}\gamma_{\bm{a}_{1}}\left(\bm{x}\right)\log\left(\tilde{w}_{1}+\tilde{w}_{2}\frac{\gamma_{\bm{a}_{2}}(\bm{x})}{\gamma_{\bm{a}_{1}}(\bm{x})}\right)\mathrm{\leavevmode\nobreak\ d}\bm{x}$ (34) $\displaystyle\leqslant$ $\displaystyle\ \log\left(\tilde{w}_{1}+\tilde{w}_{2}\exp\left(-\frac{\left\|\bm{a}_{1}-\bm{a}_{2}\right\|^{2}}{4}\right)\right)\leqslant\log\left(\tilde{w}_{1}+\exp\left(-\frac{\left\|\bm{a}_{1}-\bm{a}_{2}\right\|^{2}}{4}\right)\right).$ $$

Note that

$$ $\displaystyle\int_{\widetilde{\Omega}_{1}^{c}}\frac{1}{\left(\sqrt{2\pi}\right)^{d}}\gamma_{\bm{a}_{1}}\left(\bm{x}\right)\mathrm{\leavevmode\nobreak\ d}\bm{x}\leqslant$ $\displaystyle\int_{\widetilde{\Omega}_{1}^{c}}\exp\left(-\frac{\left\|\bm{a}_{1}-\bm{a}_{2}\right\|^{2}}{64}\right)\cdot\frac{1}{\left(\sqrt{2\pi}\right)^{d}}\exp\left(-\frac{\left\|\bm{x}-\bm{a}_{1}\right\|^{2}}{4}\right)\mathrm{\leavevmode\nobreak\ d}\bm{x}$ (35) $\displaystyle\leqslant$ $\displaystyle\exp\left(\frac{d}{2}\log 2-\frac{\left\|\bm{a}_{1}-\bm{a}_{2}\right\|^{2}}{64}\right),$ $$

and

$$ $\displaystyle\int_{\widetilde{\Omega}_{1}^{c}}\frac{1}{\left(\sqrt{2\pi}\right)^{d}}\gamma_{\bm{a}_{1}}\left(\bm{x}\right)\frac{\left\|\bm{x}-\bm{a}_{1}\right\|^{2}}{2}\mathrm{\leavevmode\nobreak\ d}\bm{x}$ (36) $\displaystyle\leqslant$ $\displaystyle\int_{\widetilde{\Omega}_{1}^{c}}\exp\left(-\frac{\left\|\bm{a}_{1}-\bm{a}_{2}\right\|^{2}}{64}\right)\cdot\frac{1}{\left(\sqrt{2\pi}\right)^{d}}\exp\left(-\frac{\left\|\bm{x}-\bm{a}_{1}\right\|^{2}}{4}\right)\frac{\left\|\bm{x}-\bm{a}_{1}\right\|^{2}}{2}\mathrm{\leavevmode\nobreak\ d}\bm{x}$ $\displaystyle\leqslant$ $\displaystyle\exp\left(\frac{d}{2}\log 2-\frac{\left\|\bm{a}_{1}-\bm{a}_{2}\right\|^{2}}{64}\right)\cdot\mathbb{E}_{\bm{x}\sim\mathcal{N}(\bm{a}_{1},2I_{d})}\left[\frac{\left\|\bm{x}-\bm{a}_{1}\right\|^{2}}{2}\right]$ $\displaystyle=$ $\displaystyle\exp\left(\log d+\frac{d}{2}\log 2-\frac{\left\|\bm{a}_{1}-\bm{a}_{2}\right\|^{2}}{64}\right).$ $$

For every $\bm{x}\in\widetilde{\Omega}_{2}$, we have $\log\left(\tilde{w}_{1}+\tilde{w}_{2}\gamma_{\bm{a}_{2}}(\bm{x})/\gamma_{\bm{a}_{1}}(\bm{x})\right)\leqslant\log 2$. Using ([35](#A1.E35)) and ([36](#A1.E36)), we obtain

$$ $\displaystyle\int_{\widetilde{\Omega}_{2}}\frac{1}{\left(\sqrt{2\pi}\right)^{d}}\gamma_{\bm{a}_{1}}\left(\bm{x}\right)\log\left(\tilde{w}_{1}+\tilde{w}_{2}\frac{\gamma_{\bm{a}_{2}}(\bm{x})}{\gamma_{\bm{a}_{1}}(\bm{x})}\right)\mathrm{\leavevmode\nobreak\ d}\bm{x}\leqslant\log 2\cdot\exp\left(\frac{d}{2}\log 2-\frac{\left\|\bm{a}_{1}-\bm{a}_{2}\right\|^{2}}{64}\right).$ (37) $$

Similarly, for any $\bm{x}\in\widetilde{\Omega}_{3}$, we have $\log\left(\tilde{w}_{1}+\tilde{w}_{2}\gamma_{\bm{a}_{2}}(\bm{x})/\gamma_{\bm{a}_{1}}(\bm{x})\right)\leqslant\log 2+\frac{\left\|\bm{x}-\bm{a}_{1}\right\|^{2}}{2}$. Thus,

$$ $\displaystyle\int_{\widetilde{\Omega}_{3}}\frac{1}{\left(\sqrt{2\pi}\right)^{d}}\gamma_{\bm{a}_{1}}\left(\bm{x}\right)\log\left(\tilde{w}_{1}+\tilde{w}_{2}\frac{\gamma_{\bm{a}_{2}}(\bm{x})}{\gamma_{\bm{a}_{1}}(\bm{x})}\right)\mathrm{\leavevmode\nobreak\ d}\bm{x}\leqslant\left(\log 2+d\right)\exp\left(\frac{d}{2}\log 2-\frac{\left\|\bm{a}_{1}-\bm{a}_{2}\right\|^{2}}{64}\right).$ (38) $$

Putting ([34](#A1.E34)), ([37](#A1.E37)), and ([38](#A1.E38)) together, $\mathbb{E}_{\gamma_{\bm{a}_{1}}}\left[\log\left(\tilde{w}_{1}+\tilde{w}_{2}\gamma_{\bm{a}_{2}}(\bm{x})/\gamma_{\bm{a}_{1}}(\bm{x})\right)\right]$ is upper bounded by

$$ $\displaystyle\mathbb{E}_{\gamma_{\bm{a}_{1}}}\left[\log\left(\tilde{w}_{1}+\tilde{w}_{2}\frac{\gamma_{\bm{a}_{2}}(\bm{x})}{\gamma_{\bm{a}_{1}}(\bm{x})}\right)\right]$ (39) $\displaystyle=$ $\displaystyle\left(\int_{\widetilde{\Omega}_{1}}+\int_{\widetilde{\Omega}_{2}}+\int_{\widetilde{\Omega}_{3}}\right)\frac{1}{\left(\sqrt{2\pi}\right)^{d}}\gamma_{\bm{a}_{1}}\left(\bm{x}\right)\log\left(\tilde{w}_{1}+\tilde{w}_{2}\frac{\gamma_{\bm{a}_{2}}(\bm{x})}{\gamma_{\bm{a}_{1}}(\bm{x})}\right)\mathrm{\leavevmode\nobreak\ d}\bm{x}$ $\displaystyle\leqslant$ $\displaystyle\log\left(\tilde{w}_{1}+\exp\left(-\frac{\left\|\bm{a}_{1}-\bm{a}_{2}\right\|^{2}}{4}\right)\right)+\left(\log 4+d\right)\exp\left(\frac{d}{2}\log 2-\frac{\left\|\bm{a}_{1}-\bm{a}_{2}\right\|^{2}}{64}\right).$ $$

Plugging ([39](#A1.E39)) into ([33](#A1.E33)), we have

$$ $\displaystyle\ w_{1}\mathbb{E}_{\gamma_{\bm{a}_{1}}}\left[\log\left(\frac{w_{1}\gamma_{\bm{a}_{1}}(\bm{x})+w_{2}\gamma_{\bm{a}_{2}}(\bm{x})}{\tilde{w}_{1}\gamma_{\bm{a}_{1}}(\bm{x})+\tilde{w}_{2}\gamma_{\bm{a}_{2}}(\bm{x})}\right)\right]$ (40) $\displaystyle\geqslant$ $\displaystyle\ w_{1}\left[\log w_{1}-\log\left(\tilde{w}_{1}+\exp\left(-\frac{\left\|\bm{a}_{1}-\bm{a}_{2}\right\|^{2}}{4}\right)\right)\right]$ $\displaystyle\quad-w_{1}\left(\log 4+d\right)\exp\left(\frac{d}{2}\log 2-\frac{\left\|\bm{a}_{1}-\bm{a}_{2}\right\|^{2}}{64}\right).$ $$

Similarly, we have

$$ $\displaystyle\ w_{2}\mathbb{E}_{\gamma_{\bm{a}_{2}}}\left[\log\left(\frac{w_{1}\gamma_{\bm{a}_{1}}(\bm{x})+w_{2}\gamma_{\bm{a}_{2}}(\bm{x})}{\tilde{w}_{1}\gamma_{\bm{a}_{1}}(\bm{x})+\tilde{w}_{2}\gamma_{\bm{a}_{2}}(\bm{x})}\right)\right]$ (41) $\displaystyle\geqslant$ $\displaystyle\ w_{2}\left(\log w_{2}-\log\left(\tilde{w}_{2}+\exp\left(-\frac{\left\|\bm{a}_{1}-\bm{a}_{2}\right\|^{2}}{4}\right)\right)\right)$ $\displaystyle\quad-w_{2}\left(\log 4+d\right)\exp\left(\frac{d}{2}\log 2-\frac{\left\|\bm{a}_{1}-\bm{a}_{2}\right\|^{2}}{64}\right).$ $$

Plugging ([40](#A1.E40)) and ([41](#A1.E41)) into ([32](#A1.E32)), we obtain the lower bound ([29](#A1.E29)) in Theorem [6](#Thmtheorem6).
Then we prove ([30](#A1.E30)).
Using ([28](#A1.E28)), we obtain

$$ $\displaystyle\nabla_{\bm{x}}\log\pi^{M}(\bm{x})-\nabla_{\bm{x}}\log\tilde{\pi}^{M}(\bm{x})$ (42) $\displaystyle=$ $\displaystyle\frac{w_{1}\bm{a}_{1}\gamma_{\bm{a}_{1}}(\bm{x})+w_{2}\bm{a}_{2}\gamma_{\bm{a}_{2}}(\bm{x})}{w_{1}\gamma_{\bm{a}_{1}}(\bm{x})+w_{2}\gamma_{\bm{a}_{2}}(\bm{x})}-\frac{\tilde{w}_{1}\bm{a}_{1}\gamma_{\bm{a}_{1}}(\bm{x})+\tilde{w}_{2}\bm{a}_{2}\gamma_{\bm{a}_{2}}(\bm{x})}{\tilde{w}_{1}\gamma_{\bm{a}_{1}}(\bm{x})+\tilde{w}_{2}\gamma_{\bm{a}_{2}}(\bm{x})}$ $\displaystyle=$ $\displaystyle\frac{w_{1}\bm{a}_{1}+w_{2}\bm{a}_{2}\gamma_{\bm{a}_{2}}(\bm{x})/\gamma_{\bm{a}_{1}}(\bm{x})}{w_{1}+w_{2}\gamma_{\bm{a}_{2}}(\bm{x})/\gamma_{\bm{a}_{1}}(\bm{x})}-\frac{\tilde{w}_{1}\bm{a}_{1}+\tilde{w}_{2}\bm{a}_{2}\gamma_{\bm{a}_{2}}(\bm{x})/\gamma_{\bm{a}_{1}}(\bm{x})}{\tilde{w}_{1}+\tilde{w}_{2}\gamma_{\bm{a}_{2}}(\bm{x})/\gamma_{\bm{a}_{1}}(\bm{x})}$ $\displaystyle=$ $\displaystyle\frac{w_{1}\bm{a}_{1}\gamma_{\bm{a}_{1}}(\bm{x})/\gamma_{\bm{a}_{2}}(\bm{x})+w_{2}\bm{a}_{2}}{w_{1}\gamma_{\bm{a}_{1}}(\bm{x})/\gamma_{\bm{a}_{2}}(\bm{x})+w_{2}}-\frac{\tilde{w}_{1}\bm{a}_{1}\gamma_{\bm{a}_{1}}(\bm{x})/\gamma_{\bm{a}_{2}}(\bm{x})+\tilde{w}_{2}\bm{a}_{2}}{\tilde{w}_{1}\gamma_{\bm{a}_{1}}(\bm{x})/\gamma_{\bm{a}_{2}}(\bm{x})+\tilde{w}_{2}}.$ $$

Recall that $\Omega_{1}=\{\bm{x}\in\mathbb{R}^{d}:\left\|\bm{x}-\bm{a}_{1}\right\|\leqslant\frac{\left\|\bm{a}_{1}-\bm{a}_{2}\right\|}{4}\}$, $\Omega_{2}=\{\bm{x}\in\mathbb{R}^{d}:\left\|\bm{x}-\bm{a}_{2}\right\|\leqslant\frac{\left\|\bm{a}_{1}-\bm{a}_{2}\right\|}{4}\}$ and $\Omega_{3}=\Omega_{1}^{c}\bigcap\Omega_{2}^{c}$.
For any $\bm{x}\in\Omega_{1}$, we can rewrite ([42](#A1.E42)) as

$$ $\displaystyle\nabla_{\bm{x}}\log\pi^{M}(\bm{x})-\nabla_{\bm{x}}\log\tilde{\pi}^{M}(\bm{x})$ (43) $\displaystyle=$ $\displaystyle\bm{a}_{1}+\frac{w_{2}(\bm{a}_{2}-\bm{a}_{1})\gamma_{\bm{a}_{2}}(\bm{x})/\gamma_{\bm{a}_{1}}(\bm{x})}{w_{1}+w_{2}\gamma_{\bm{a}_{2}}(\bm{x})/\gamma_{\bm{a}_{1}}(\bm{x})}-\left\{\bm{a}_{1}+\frac{\tilde{w}_{2}(\bm{a}_{2}-\bm{a}_{1})\gamma_{\bm{a}_{2}}(\bm{x})/\gamma_{\bm{a}_{1}}(\bm{x})}{\tilde{w}_{1}+\tilde{w}_{2}\gamma_{\bm{a}_{2}}(\bm{x})/\gamma_{\bm{a}_{1}}(\bm{x})}\right\}$ $\displaystyle=$ $\displaystyle\frac{w_{2}(\bm{a}_{2}-\bm{a}_{1})\gamma_{\bm{a}_{2}}(\bm{x})/\gamma_{\bm{a}_{1}}(\bm{x})}{w_{1}+w_{2}\gamma_{\bm{a}_{2}}(\bm{x})/\gamma_{\bm{a}_{1}}(\bm{x})}-\frac{\tilde{w}_{2}(\bm{a}_{2}-\bm{a}_{1})\gamma_{\bm{a}_{2}}(\bm{x})/\gamma_{\bm{a}_{1}}(\bm{x})}{\tilde{w}_{1}+\tilde{w}_{2}\gamma_{\bm{a}_{2}}(\bm{x})/\gamma_{\bm{a}_{1}}(\bm{x})}.$ $$

Note that $|\gamma_{\bm{a}_{2}}(\bm{x})/\gamma_{\bm{a}_{1}}(\bm{x})|^{2}=\exp(\|\bm{x}-\bm{a}_{1}\|^{2}-\|\bm{x}-\bm{a}_{2}\|^{2})\leqslant\exp(-\|\bm{a}_{1}-\bm{a}_{2}\|^{2}/2)$ for every $\bm{x}\in\Omega_{1}$. Then use ([43](#A1.E43)), we obtain

$$ $\displaystyle\int_{\Omega_{1}}\left\|\nabla_{\bm{x}}\log\pi^{M}(\bm{x})-\nabla_{\bm{x}}\log\tilde{\pi}^{M}(\bm{x})\right\|^{2}p(\bm{x})\mathrm{\leavevmode\nobreak\ d}\bm{x}$ (44) $\displaystyle=$ $\displaystyle\int_{\Omega_{1}}\left\|\frac{w_{2}(\bm{a}_{2}-\bm{a}_{1})\gamma_{\bm{a}_{2}}(\bm{x})/\gamma_{\bm{a}_{1}}(\bm{x})}{w_{1}+w_{2}\gamma_{\bm{a}_{2}}(\bm{x})/\gamma_{\bm{a}_{1}}(\bm{x})}-\frac{\tilde{w}_{2}(\bm{a}_{2}-\bm{a}_{1})\gamma_{\bm{a}_{2}}(\bm{x})/\gamma_{\bm{a}_{1}}(\bm{x})}{\tilde{w}_{1}+\tilde{w}_{2}\gamma_{\bm{a}_{2}}(\bm{x})/\gamma_{\bm{a}_{1}}(\bm{x})}\right\|^{2}p(\bm{x})\mathrm{\leavevmode\nobreak\ d}\bm{x}$ $\displaystyle\leqslant$ $\displaystyle\ 2\int_{\Omega_{1}}\left\|\frac{w_{2}(\bm{a}_{2}-\bm{a}_{1})}{w_{1}}\right\|^{2}\Big{|}\frac{\gamma_{\bm{a}_{2}}(\bm{x})}{\gamma_{\bm{a}_{1}}(\bm{x})}\Big{|}^{2}p(\bm{x})\mathrm{\leavevmode\nobreak\ d}\bm{x}+2\int_{\Omega_{1}}\left\|\frac{\tilde{w}_{2}(\bm{a}_{2}-\bm{a}_{1})}{\tilde{w}_{1}}\right\|^{2}\Big{|}\frac{\gamma_{\bm{a}_{2}}(\bm{x})}{\gamma_{\bm{a}_{1}}(\bm{x})}\Big{|}^{2}p(\bm{x})\mathrm{\leavevmode\nobreak\ d}\bm{x}$ $\displaystyle\leqslant$ $\displaystyle\ 2\exp\left(-\frac{\left\|\bm{a}_{1}-\bm{a}_{2}\right\|^{2}}{2}\right)\left[\frac{w_{2}^{2}}{w_{1}^{2}}+\frac{\tilde{w}_{2}^{2}}{\tilde{w}_{1}^{2}}\right]\left\|\bm{a}_{2}-\bm{a}_{1}\right\|^{2}.$ $$

Similarly, we obtain

$$ $\displaystyle\int_{\Omega_{2}}\left\|\nabla_{\bm{x}}\log\pi^{M}(\bm{x})-\nabla_{\bm{x}}\log\tilde{\pi}^{M}(\bm{x})\right\|^{2}p(\bm{x})\mathrm{\leavevmode\nobreak\ d}\bm{x}$ (45) $\displaystyle\leqslant$ $\displaystyle\ 2\exp\left(-\frac{\left\|\bm{a}_{1}-\bm{a}_{2}\right\|^{2}}{2}\right)\left[\frac{w_{1}^{2}}{w_{2}^{2}}+\frac{\tilde{w}_{1}^{2}}{\tilde{w}_{2}^{2}}\right]\left\|\bm{a}_{1}-\bm{a}_{2}\right\|^{2}.$ $$

Using ([42](#A1.E42)), we obtain that

$$ $\displaystyle\int_{\Omega_{3}}\left\|\nabla_{\bm{x}}\log\pi^{M}(\bm{x})-\nabla_{\bm{x}}\log\tilde{\pi}^{M}(\bm{x})\right\|^{2}p(\bm{x})\mathrm{\leavevmode\nobreak\ d}\bm{x}$ (46) $\displaystyle=$ $\displaystyle\int_{\Omega_{3}}\left\|\frac{w_{1}\bm{a}_{1}+w_{2}\bm{a}_{2}\gamma_{\bm{a}_{2}}(\bm{x})/\gamma_{\bm{a}_{1}}(\bm{x})}{w_{1}+w_{2}\gamma_{\bm{a}_{2}}(\bm{x})/\gamma_{\bm{a}_{1}}(\bm{x})}-\frac{\tilde{w}_{1}\bm{a}_{1}+\tilde{w}_{2}\bm{a}_{2}\gamma_{\bm{a}_{2}}(\bm{x})/\gamma_{\bm{a}_{1}}(\bm{x})}{\tilde{w}_{1}+\tilde{w}_{2}\gamma_{\bm{a}_{2}}(\bm{x})/\gamma_{\bm{a}_{1}}(\bm{x})}\right\|^{2}p(\bm{x})\mathrm{\leavevmode\nobreak\ d}\bm{x}$ $\displaystyle\leqslant$ $\displaystyle\ 4\int_{\Omega_{3}}\left\|\frac{w_{1}\bm{a}_{1}}{w_{1}+w_{2}\gamma_{\bm{a}_{2}}(\bm{x})/\gamma_{\bm{a}_{1}}(\bm{x})}\right\|^{2}p(\bm{x})\mathrm{\leavevmode\nobreak\ d}\bm{x}+4\int_{\Omega_{3}}\left\|\frac{\tilde{w}_{1}\bm{a}_{1}}{\tilde{w}_{1}+\tilde{w}_{2}\gamma_{\bm{a}_{2}}(\bm{x})/\gamma_{\bm{a}_{1}}(\bm{x})}\right\|^{2}p(\bm{x})\mathrm{\leavevmode\nobreak\ d}\bm{x}$ $\displaystyle\quad+\ 4\int_{\Omega_{3}}\left\|\frac{w_{2}\bm{a}_{2}\gamma_{\bm{a}_{2}}(\bm{x})/\gamma_{\bm{a}_{1}}(\bm{x})}{w_{1}+w_{2}\gamma_{\bm{a}_{2}}(\bm{x})/\gamma_{\bm{a}_{1}}(\bm{x})}\right\|^{2}p(\bm{x})\mathrm{\leavevmode\nobreak\ d}\bm{x}+4\int_{\Omega_{3}}\left\|\frac{\tilde{w}_{2}\bm{a}_{2}\gamma_{\bm{a}_{2}}(\bm{x})/\gamma_{\bm{a}_{1}}(\bm{x})}{\tilde{w}_{1}+\tilde{w}_{2}\gamma_{\bm{a}_{2}}(\bm{x})/\gamma_{\bm{a}_{1}}(\bm{x})}\right\|^{2}p(\bm{x})\mathrm{\leavevmode\nobreak\ d}\bm{x}$ $\displaystyle\leqslant$ $\displaystyle\ 8\left[\left\|\bm{a}_{1}\right\|^{2}+\left\|\bm{a}_{2}\right\|^{2}\right]\int_{\Omega_{3}}p(\bm{x})\mathrm{\leavevmode\nobreak\ d}\bm{x}.$ $$

Note that we have the following decomposition

$$ $\displaystyle\mathrm{SE}_{p}(\pi^{M}\|\tilde{\pi}^{M})=$ $\displaystyle\int_{\mathbb{R}^{d}}\left\|\nabla_{\bm{x}}\log\pi^{M}(\bm{x})-\nabla_{\bm{x}}\log\tilde{\pi}^{M}(\bm{x})\right\|^{2}p(\bm{x})\mathrm{\leavevmode\nobreak\ d}\bm{x}$ (47) $\displaystyle=$ $\displaystyle\left(\int_{\Omega_{1}}+\int_{\Omega_{2}}+\int_{\Omega_{3}}\right)\left\|\nabla_{\bm{x}}\log\pi^{M}(\bm{x})-\nabla_{\bm{x}}\log\tilde{\pi}^{M}(\bm{x})\right\|^{2}p(\bm{x})\mathrm{\leavevmode\nobreak\ d}\bm{x}.$ $$

Plug ([44](#A1.E44)), ([45](#A1.E45)), and ([46](#A1.E46)) into ([47](#A1.E47)), we obtain the upper bound ([30](#A1.E30)) in Theorem [6](#Thmtheorem6).
∎

###### Theorem 6 .

###### Remark 3 .

###### Proof of Theorem 6 .

### A.3 Proof of Theorem 2

First, we present the divergence theorem and Green’s first identity, which is very useful in our proof.
Then we state the Grönwall’s inequality used in our proof.
Finally, we state and prove Theorem [7](#Thmtheorem7) which includes Theorem [2](#Thmtheorem2) and sharper bounds when ([52](#A1.E52)) holds.

###### Lemma 1 (divergence theorem) .

Let $\mathbf{F}(\cdot):\Omega\to\mathbb{R}^{d}$, then
$\int_{\Omega}\nabla\cdot\mathbf{F}(\bm{x})\mathrm{\leavevmode\nobreak\ d}\bm{x}=\int_{\partial\Omega}\mathbf{F}\cdot\bm{n}\mathrm{\leavevmode\nobreak\ d}\bm{S}$.

###### Lemma 2 (Green’s first identity) .

Let $v(\cdot),u(\cdot):\Omega\to\mathbb{R}$, then it holds that

$$ $\int_{\Omega}\nabla_{\bm{x}}v\cdot\nabla_{\bm{x}}u\mathrm{\leavevmode\nobreak\ d}\bm{x}+\int_{\Omega}v\Delta u\mathrm{\leavevmode\nobreak\ d}\bm{x}=\int_{\partial\Omega}v\frac{\partial u}{\partial\bm{n}}\mathrm{\leavevmode\nobreak\ d}\bm{S}.$ $$

###### Lemma 3 (Grönwall’s inequality) .

Let $f(\cdot),\alpha(\cdot),\beta(\cdot):[0,T]\to\mathbb{R}$, and suppose that $\forall\ 0\leqslant t\leqslant T$,

$$ $f^{\prime}(t)\leqslant\alpha(t)+\beta(t)f(t).$ $$

Then we have $\forall\ 0\leqslant t\leqslant T$,

$$ $f(t)\leqslant e^{\int_{0}^{t}\beta(s)\mathrm{\leavevmode\nobreak\ d}s}f(0)+\int_{0}^{t}e^{\int_{s}^{t}\beta(r)\mathrm{\leavevmode\nobreak\ d}r}\alpha(s)\mathrm{\leavevmode\nobreak\ d}s.$ $$

###### Proof of Lemma 3 .

Consider $g(t)=e^{-\int_{0}^{t}\beta(s)\mathrm{\leavevmode\nobreak\ d}s}f(t),\forall\ 0\leqslant t\leqslant T$. Then we have $\forall\ 0\leqslant t\leqslant T$,

$$ $\displaystyle g^{\prime}(t)=$ $\displaystyle\ e^{-\int_{0}^{t}\beta(s)\mathrm{\leavevmode\nobreak\ d}s}f^{\prime}(t)-\beta(t)e^{-\int_{0}^{t}\beta(s)\mathrm{\leavevmode\nobreak\ d}s}f(t)$ (48) $\displaystyle=$ $\displaystyle\ e^{-\int_{0}^{t}\beta(s)\mathrm{\leavevmode\nobreak\ d}s}\left(f^{\prime}(t)-\beta(t)f(t)\right)$ $\displaystyle\leqslant$ $\displaystyle\ e^{-\int_{0}^{t}\beta(s)\mathrm{\leavevmode\nobreak\ d}s}\alpha(t).$ $$

Integrating ([48](#A1.E48)), we obtain

$$ $e^{-\int_{0}^{t}\beta(s)\mathrm{\leavevmode\nobreak\ d}s}f(t)\leqslant f(0)+\int_{0}^{t}e^{-\int_{0}^{s}\beta(r)\mathrm{\leavevmode\nobreak\ d}r}\alpha(s)\mathrm{\leavevmode\nobreak\ d}s.$ (49) $$

Hence, we complete our proof.
∎

###### Theorem 7 .

Suppose that Assumption [1](#Thmassumption1), [2](#Thmassumption2), and [3](#Thmassumption3) hold.
We further assume that $u_{\theta}(\bm{x},0)=u^{*}_{0}(\bm{x})$ for any $\bm{x}\in\Omega$.
Then for any positive constant $\varepsilon>0$, the following holds for any $0\leqslant t\leqslant T$,

$$ $\|e_{t}(\cdot)\|^{2}_{L^{2}(\Omega;\nu_{t})}\leqslant\varepsilon L_{\emph{{PINN}}}(t;C_{1}(\varepsilon)).$ (50) $$

Moreover, for any $0\leqslant t\leqslant T$,

$$ $m_{1}\|\nabla_{\bm{x}}e_{t}(\cdot)\|^{2}_{L^{2}(\Omega;\nu_{t})}\leqslant\varepsilon\|r_{t}(\cdot)\|^{2}_{L^{2}(\Omega;\nu_{t})}+C_{3}(\varepsilon)L_{\emph{{PINN}}}(t;C_{1}(\varepsilon))+C_{2}\sqrt{\varepsilon L_{\emph{{PINN}}}(t;C_{1}(\varepsilon))}.$ (51) $$

In addition, if there exists constant $\mathcal{C}_{\nu}(\Omega)>0$ such that the following holds for any $0\leqslant t\leqslant T$,

$$ $\|\nabla_{\bm{x}}e_{t}(\cdot)\|^{2}_{L^{2}(\Omega;\nu_{t})}\geqslant\mathcal{C}^{2}_{\nu}(\Omega)\|e_{t}(\cdot)\|^{2}_{L^{2}(\Omega;\nu_{t})}.$ (52) $$

Then for any positive constant $\varepsilon>0$, the following holds for any $0\leqslant t\leqslant T$,

$$ $\|e_{t}(\cdot)\|^{2}_{L^{2}(\Omega;\nu_{t})}\leqslant\varepsilon L_{\emph{{PINN}}}(t;C_{4}(\varepsilon)).$ (53) $$

Moreover, for any $0\leqslant t\leqslant T$,

$$ $m_{1}\|\nabla_{\bm{x}}e_{t}(\cdot)\|^{2}_{L^{2}(\Omega;\nu_{t})}\leqslant\varepsilon\|r_{t}(\cdot)\|^{2}_{L^{2}(\Omega;\nu_{t})}+C_{3}(\varepsilon)L_{\emph{{PINN}}}(t;C_{4}(\varepsilon))+C_{2}\sqrt{\varepsilon L_{\emph{{PINN}}}(t;C_{4}(\varepsilon))}.$ (54) $$

where $C_{2}:=2\sqrt{2}(\widehat{B}_{0}^{2}+B_{0}^{*2})^{1/2}$, $C_{3}(\varepsilon):=\varepsilon(C_{1}(\varepsilon)+B^{\nu}_{0})$, $C_{4}(\varepsilon):=C_{1}(\varepsilon)-m_{1}\mathcal{C}^{2}_{\nu}(\Omega)$ and

$$ $C_{1}(\varepsilon):=\frac{1}{\varepsilon}+\frac{M_{1}}{4}(B^{\nu}_{1}+2B^{*}_{1}+2\widehat{B}_{1})+c_{1}(B^{\nu}_{1}+B^{\nu}_{2})-\frac{c_{2}}{2}(B^{*}_{2}+\widehat{B}_{2}),$ $$

where

$$ $c_{1}:=\begin{cases}M_{1},\quad\text{if }B^{\nu}_{1}+B^{\nu}_{2}\geqslant 0.\\ m_{1},\quad\text{if }B^{\nu}_{1}+B^{\nu}_{2}<0.\end{cases},\quad c_{2}:=\begin{cases}m_{1},\quad\text{if }B^{*}_{2}+\widehat{B}_{2}\geqslant 0.\\ M_{1},\quad\text{if }B^{*}_{2}+\widehat{B}_{2}<0.\end{cases}.$ $$

###### Proof of Theorem 7 .

We first prove ([50](#A1.E50)) and ([53](#A1.E53)).
Note that $u_{t}^{*}(\bm{x})$ satisfies

$$ $\partial_{t}u^{*}_{t}(\bm{x})+\nabla_{\bm{x}}u^{*}_{t}(\bm{x})\cdot\bm{f}(\bm{x},t)+\nabla\cdot\bm{f}(\bm{x},t)-\frac{1}{2}g^{2}(t)\Delta u_{t}^{*}(\bm{x})-\frac{1}{2}g^{2}(t)\left\|\nabla_{\bm{x}}u_{t}^{*}(\bm{x})\right\|^{2}=0.$ (55) $$

And $u_{\theta}(\bm{x},t)$ satisfies

$$ $\partial_{t}u_{\theta}(\bm{x},t)+\nabla_{\bm{x}}u_{\theta}(\bm{x},t)\cdot\bm{f}(\bm{x},t)+\nabla\cdot\bm{f}(\bm{x},t)-\frac{1}{2}g^{2}(t)\Delta u_{\theta}(\bm{x},t)-\frac{1}{2}g^{2}(t)\left\|\nabla_{\bm{x}}u_{\theta}(\bm{x},t)\right\|^{2}=r_{t}(\bm{x}).$ (56) $$

Subtracting ([55](#A1.E55)) for $u^{*}$ from ([56](#A1.E56)) for $u_{\theta}$, we have

$$ $\partial_{t}e_{t}(\bm{x})+\nabla_{\bm{x}}e_{t}(\bm{x})\cdot\bm{f}(\bm{x},t)-\frac{1}{2}g^{2}(t)\left(\left\|\nabla_{\bm{x}}u_{\theta}(\bm{x},t)\right\|^{2}-\left\|\nabla_{\bm{x}}u^{*}_{t}(\bm{x})\right\|^{2}\right)-\frac{1}{2}g^{2}(t)\Delta e_{t}(\bm{x})=r_{t}(\bm{x}).$ (57) $$

Note that $\frac{1}{2}\partial_{t}e_{t}^{2}(\bm{x})=e_{t}(\bm{x})\partial_{t}e_{t}(\bm{x})$ and $\frac{1}{2}\nabla_{\bm{x}}e_{t}^{2}(\bm{x})=e_{t}(\bm{x})\nabla_{\bm{x}}e_{t}(\bm{x})$, then we obtain

$$ $\displaystyle\frac{1}{2}\partial_{t}e_{t}^{2}(\bm{x})=$ $\displaystyle\ \frac{1}{2}g^{2}(t)e_{t}(\bm{x})\left(\left\|\nabla_{\bm{x}}u_{\theta}(\bm{x},t)\right\|^{2}-\left\|\nabla_{\bm{x}}u^{*}_{t}(\bm{x})\right\|^{2}\right)+\frac{1}{2}g^{2}(t)e_{t}(\bm{x})\Delta e_{t}(\bm{x})$ (58) $\displaystyle\quad+e_{t}(\bm{x})r_{t}(\bm{x})-e_{t}(\bm{x})\nabla_{\bm{x}}e_{t}(\bm{x})\cdot\bm{f}(\bm{x},t)$ $\displaystyle=$ $\displaystyle\ \frac{1}{2}g^{2}(t)e_{t}(\bm{x})\nabla_{\bm{x}}e_{t}(\bm{x})\cdot\left(\nabla_{\bm{x}}u_{\theta}(\bm{x},t)+\nabla_{\bm{x}}u^{*}_{t}(\bm{x})\right)+\frac{1}{2}g^{2}(t)e_{t}(\bm{x})\Delta e_{t}(\bm{x})$ $\displaystyle\quad+e_{t}(\bm{x})r_{t}(\bm{x})-e_{t}(\bm{x})\nabla_{\bm{x}}e_{t}(\bm{x})\cdot\bm{f}(\bm{x},t)$ $\displaystyle=$ $\displaystyle\ \frac{1}{4}g^{2}(t)\nabla_{\bm{x}}e_{t}^{2}(\bm{x})\cdot\left(\nabla_{\bm{x}}u_{\theta}(\bm{x},t)+\nabla_{\bm{x}}u^{*}_{t}(\bm{x})\right)+\frac{1}{2}g^{2}(t)e_{t}(\bm{x})\Delta e_{t}(\bm{x})$ $\displaystyle\quad+e_{t}(\bm{x})r_{t}(\bm{x})-\frac{1}{2}\nabla_{\bm{x}}e_{t}^{2}(\bm{x})\cdot f(\bm{x},t).$ $$

Note that $\partial_{t}(\nu_{t}(\bm{x})e_{t}^{2}(\bm{x}))=e_{t}^{2}(\bm{x})\partial_{t}\nu_{t}(\bm{x})+\nu_{t}(\bm{x})\partial_{t}e_{t}^{2}(\bm{x})$, then we have

$$ $\displaystyle\partial_{t}(\nu_{t}(\bm{x})e^{2}_{t}(\bm{x}))=$ $\displaystyle\ \frac{1}{2}g^{2}(t)\nu_{t}(\bm{x})\nabla_{\bm{x}}e_{t}^{2}(\bm{x})\cdot\left(\nabla_{\bm{x}}u_{\theta}(\bm{x},t)+\nabla_{\bm{x}}u^{*}_{t}(\bm{x})\right)$ (59) $\displaystyle\quad+g^{2}(t)\nu_{t}(\bm{x})e_{t}(\bm{x})\Delta e_{t}(\bm{x})$ $\displaystyle\quad+2\nu_{t}(\bm{x})e_{t}(\bm{x})r_{t}(\bm{x})-\nu_{t}(\bm{x})\nabla_{\bm{x}}e_{t}^{2}(\bm{x})\cdot\bm{f}(\bm{x},t)$ $\displaystyle\quad+\frac{1}{2}g^{2}(t)e_{t}^{2}(\bm{x})\Delta\nu_{t}(\bm{x})-e_{t}^{2}(\bm{x})\nabla\cdot\left[\bm{f}(\bm{x},t)\nu_{t}(\bm{x})\right].$ $$

We integrate ([59](#A1.E59)) to get an equation for $\|e_{t}(\cdot)\|^{2}_{L^{2}(\Omega;\nu_{t})}$ given by

$$ $\displaystyle\partial_{t}\left\|e_{t}(\cdot)\right\|^{2}_{L^{2}(\Omega;\nu_{t})}=$ $\displaystyle\ \frac{1}{2}g^{2}(t)\int_{\Omega}\nu_{t}(\bm{x})\nabla_{\bm{x}}e_{t}^{2}(\bm{x})\cdot\left(\nabla_{\bm{x}}u_{\theta}(\bm{x},t)+\nabla_{\bm{x}}u^{*}_{t}(\bm{x})\right)\mathrm{\leavevmode\nobreak\ d}\bm{x}$ (60) $\displaystyle\quad+\ g^{2}(t)\int_{\Omega}\nu_{t}(\bm{x})e_{t}(\bm{x})\Delta e_{t}(\bm{x})\mathrm{\leavevmode\nobreak\ d}\bm{x}$ $\displaystyle\quad+2\int_{\Omega}\nu_{t}(\bm{x})e_{t}(\bm{x})r_{t}(\bm{x})\mathrm{\leavevmode\nobreak\ d}\bm{x}-\int_{\Omega}\nu_{t}(\bm{x})\nabla_{\bm{x}}e_{t}^{2}(\bm{x})\cdot\bm{f}(\bm{x},t)\mathrm{\leavevmode\nobreak\ d}\bm{x}$ $\displaystyle\quad+\frac{1}{2}g^{2}(t)\int_{\Omega}e_{t}^{2}(\bm{x})\Delta\nu_{t}(\bm{x})\mathrm{\leavevmode\nobreak\ d}\bm{x}-\int_{\Omega}e_{t}^{2}(\bm{x})\nabla\cdot\left[\bm{f}(\bm{x},t)\nu_{t}(\bm{x})\right]\mathrm{\leavevmode\nobreak\ d}\bm{x}.$ $$

Note that

$$ $\nabla\cdot\left[\nu_{t}(\bm{x})e_{t}^{2}(\bm{x})\bm{f}(\bm{x},t)\right]=\nu_{t}(\bm{x})\nabla_{\bm{x}}e_{t}^{2}(\bm{x})\cdot\bm{f}(\bm{x},t)+e_{t}^{2}(\bm{x})\nabla\cdot\left[\nu_{t}(\bm{x})\bm{f}(\bm{x},t)\right].$ $$

Then using Lemma [1](#Thmlemma1) and $e_{t}(\bm{x})=0$ for any $(\bm{x},t)\in\partial\Omega\times[0,T]$, we have

$$ $\int_{\Omega}\nu_{t}(\bm{x})\nabla_{\bm{x}}e_{t}^{2}(\bm{x})\cdot\bm{f}(\bm{x},t)\mathrm{\leavevmode\nobreak\ d}\bm{x}+\int_{\Omega}e_{t}^{2}(\bm{x})\nabla\cdot\left[\nu_{t}(\bm{x})\bm{f}(\bm{x},t)\right]\mathrm{\leavevmode\nobreak\ d}\bm{x}=0.$ (61) $$

Similarly, we have

$$ $\displaystyle\int_{\Omega}\nu_{t}(\bm{x})\nabla_{\bm{x}}e_{t}^{2}(\bm{x})\cdot\left(\nabla_{\bm{x}}u_{\theta}(\bm{x},t)+\nabla_{\bm{x}}u^{*}_{t}(\bm{x})\right)\mathrm{\leavevmode\nobreak\ d}\bm{x}$ (62) $\displaystyle=$ $\displaystyle-\int_{\Omega}\nu_{t}(\bm{x})e_{t}^{2}(\bm{x})\left(\Delta u_{\theta}(\bm{x},t)+\Delta u^{*}_{t}(\bm{x})\right)\mathrm{\leavevmode\nobreak\ d}\bm{x}$ $\displaystyle-\int_{\Omega}e_{t}^{2}(\bm{x})\nabla_{\bm{x}}\nu_{t}(\bm{x})\cdot\left(\nabla_{\bm{x}}u_{\theta}(\bm{x},t)+\nabla_{\bm{x}}u^{*}_{t}(\bm{x})\right)\mathrm{\leavevmode\nobreak\ d}\bm{x},$ $$

and

$$ $\int_{\Omega}\nu_{t}(\bm{x})e_{t}(\bm{x})\Delta e_{t}(\bm{x})\mathrm{\leavevmode\nobreak\ d}\bm{x}=-\frac{1}{2}\int_{\Omega}\nabla_{\bm{x}}\nu_{t}(\bm{x})\cdot\nabla_{\bm{x}}e_{t}^{2}(\bm{x})\mathrm{\leavevmode\nobreak\ d}\bm{x}-\int_{\Omega}\nu_{t}(\bm{x})\left\|\nabla_{\bm{x}}e_{t}(\bm{x})\right\|^{2}\mathrm{\leavevmode\nobreak\ d}\bm{x}.$ (63) $$

Plugging ([61](#A1.E61)), ([62](#A1.E62)), and ([63](#A1.E63)) into ([60](#A1.E60)), and using Lemma [2](#Thmlemma2), we have

$$ $\displaystyle\partial_{t}\|e_{t}(\cdot)\|^{2}_{L^{2}(\Omega;\nu_{t})}=$ $\displaystyle-\frac{1}{2}g^{2}(t)\int_{\Omega}(\Delta u_{\theta}(\bm{x},t)+\Delta u^{*}_{t}(\bm{x}))e_{t}^{2}(\bm{x})\nu_{t}(\bm{x})\mathrm{\leavevmode\nobreak\ d}\bm{x}$ (64) $\displaystyle\ -\frac{1}{2}g^{2}(t)\int_{\Omega}\nabla_{\bm{x}}\nu_{t}(\bm{x})\cdot\left(\nabla_{\bm{x}}u_{\theta}(\bm{x},t)+\nabla_{\bm{x}}u_{t}^{*}(\bm{x})\right)e_{t}^{2}(\bm{x})\mathrm{\leavevmode\nobreak\ d}\bm{x}$ $\displaystyle\ -g^{2}(t)\int_{\Omega}\nu_{t}(\bm{x})\|\nabla_{\bm{x}}e_{t}(\bm{x})\|^{2}\mathrm{\leavevmode\nobreak\ d}\bm{x}+2\int_{\Omega}\nu_{t}(\bm{x})e_{t}(\bm{x})r_{t}(\bm{x})\mathrm{\leavevmode\nobreak\ d}\bm{x}$ $\displaystyle\ +g^{2}(t)\int_{\Omega}e_{t}^{2}(\bm{x})\Delta\nu_{t}(\bm{x})\mathrm{\leavevmode\nobreak\ d}\bm{x}.$ $$

Using $\nabla_{\bm{x}}\nu_{t}(\bm{x})=\nu_{t}(\bm{x})\nabla_{\bm{x}}\log\nu_{t}(\bm{x})$ and $\Delta\nu_{t}(\bm{x})=(\Delta\log\nu_{t}(\bm{x})+\|\nabla_{\bm{x}}\log\nu_{t}(\bm{x})\|^{2})\nu_{t}(\bm{x})$,

$$ $\displaystyle\partial_{t}\|e_{t}(\cdot)\|^{2}_{L^{2}(\Omega;\nu_{t})}=$ $\displaystyle-\frac{1}{2}g^{2}(t)\int_{\Omega}(\Delta u_{\theta}(\bm{x},t)+\Delta u^{*}_{t}(\bm{x}))e_{t}^{2}(\bm{x})\nu_{t}(\bm{x})\mathrm{\leavevmode\nobreak\ d}\bm{x}$ (65) $\displaystyle\ -\frac{1}{2}g^{2}(t)\int_{\Omega}\nabla_{\bm{x}}\log\nu_{t}(\bm{x})\cdot\left(\nabla_{\bm{x}}u_{\theta}(\bm{x},t)+\nabla_{\bm{x}}u_{t}^{*}(\bm{x})\right)e_{t}^{2}(\bm{x})\nu_{t}(\bm{x})\mathrm{\leavevmode\nobreak\ d}\bm{x}$ $\displaystyle\ -g^{2}(t)\int_{\Omega}\nu_{t}(\bm{x})\|\nabla_{\bm{x}}e_{t}(\bm{x})\|^{2}\mathrm{\leavevmode\nobreak\ d}\bm{x}+2\int_{\Omega}\nu_{t}(\bm{x})e_{t}(\bm{x})r_{t}(\bm{x})\mathrm{\leavevmode\nobreak\ d}\bm{x}$ $\displaystyle\ +g^{2}(t)\int_{\Omega}e_{t}^{2}(\bm{x})\nu_{t}(\bm{x})(\Delta\log\nu_{t}(\bm{x})+\|\nabla_{\bm{x}}\log\nu_{t}(\bm{x})\|^{2})\mathrm{\leavevmode\nobreak\ d}\bm{x}.$ $$

By Assumption [2](#Thmassumption2) and [3](#Thmassumption3), then we have

$$ $\displaystyle\partial_{t}\|e_{t}(\cdot)\|_{L^{2}(\Omega;\nu_{t})}^{2}\leqslant$ $\displaystyle\ \varepsilon\|r_{t}(\cdot)\|^{2}_{L^{2}(\Omega;\nu_{t})}+C_{1}(\varepsilon)\|e_{t}(\cdot)\|_{L^{2}(\Omega;\nu_{t})}^{2}-m_{1}\|\nabla_{\bm{x}}e_{t}(\cdot)\|^{2}_{L^{2}(\Omega;\nu_{t})}$ (66) $\displaystyle\leqslant$ $\displaystyle\ \varepsilon\|r_{t}(\cdot)\|^{2}_{L^{2}(\Omega;\nu_{t})}+C_{1}(\varepsilon)\|e_{t}(\cdot)\|_{L^{2}(\Omega;\nu_{t})}^{2},$ $$

which follows from applying Young’s inequality and holds for any $\varepsilon>0$. Note that $e_{0}(\bm{x})=0$ for any $\bm{x}\in\Omega$, then using Lemma [3](#Thmlemma3), we have $\forall\ 0\leqslant t\leqslant T$,

$$ $\|e_{t}(\cdot)\|^{2}_{L^{2}(\Omega;\nu_{t})}\leqslant\varepsilon\int_{0}^{t}e^{C_{1}(\varepsilon)(t-s)}\|r_{s}(\cdot)\|^{2}_{L^{2}(\Omega;\nu_{s})}\mathrm{\leavevmode\nobreak\ d}s:=\varepsilon L_{\textrm{PINN}}(t;C_{1}(\varepsilon)).$ (67) $$

Hence, we have proved ([50](#A1.E50)).
In addition, if ([52](#A1.E52)) holds, plugging ([52](#A1.E52)) into ([66](#A1.E66)), we have

$$ $\partial_{t}\|e_{t}(\cdot)\|_{L^{2}(\Omega;\nu_{t})}^{2}\leqslant\varepsilon\|r_{t}(\cdot)\|^{2}_{L^{2}(\Omega;\nu_{t})}+C_{4}(\varepsilon)\|e_{t}(\cdot)\|_{L^{2}(\Omega;\nu_{t})}^{2}.$ (68) $$

Similarly, using Lemma [3](#Thmlemma3), we obtain ([53](#A1.E53)).
Then we prove ([51](#A1.E51)) and ([54](#A1.E54)).
From ([66](#A1.E66)), we have

$$ $m_{1}\|\nabla_{\bm{x}}e_{t}(\cdot)\|^{2}_{L^{2}(\Omega;\nu_{t})}\leqslant\varepsilon\|r_{t}(\cdot)\|^{2}_{L^{2}(\Omega;\nu_{t})}+C_{1}(\varepsilon)\|e_{t}(\cdot)\|^{2}_{L^{2}(\Omega;\nu_{t})}-\partial_{t}\|e_{t}(\cdot)\|^{2}_{L^{2}(\Omega;\nu_{t})}.$ (69) $$

By Assumption [3](#Thmassumption3), we bound $\partial_{t}\|e_{t}(\cdot)\|^{2}_{L^{2}(\Omega;\nu_{t})}$ as follows

$$ $\displaystyle\Big{|}\partial_{t}\|e_{t}(\cdot)\|^{2}_{L^{2}(\Omega;\nu_{t})}\Big{|}=\Big{|}\partial_{t}\left(\int_{\Omega}e_{t}^{2}(\bm{x})\nu_{t}(\bm{x})\mathrm{\leavevmode\nobreak\ d}\bm{x}\right)\Big{|}$ (70) $\displaystyle=$ $\displaystyle\Big{|}\int_{\Omega}e_{t}^{2}(\bm{x})\partial_{t}\nu_{t}(\bm{x})+2\nu_{t}(\bm{x})e_{t}(\bm{x})\partial_{t}e_{t}(\bm{x})\mathrm{\leavevmode\nobreak\ d}\bm{x}\Big{|}$ $\displaystyle\leqslant$ $\displaystyle\int_{\Omega}e_{t}^{2}(\bm{x})\nu_{t}(\bm{x})|\partial_{t}\log\nu_{t}(\bm{x})|\mathrm{\leavevmode\nobreak\ d}\bm{x}+2\Big{|}\int_{\Omega}\nu_{t}(\bm{x})e_{t}(\bm{x})\partial_{t}e_{t}(\bm{x})\mathrm{\leavevmode\nobreak\ d}\bm{x}\Big{|}$ $\displaystyle\leqslant$ $\displaystyle\ B^{\nu}_{0}\|e_{t}(\cdot)\|^{2}_{L^{2}(\Omega;\nu_{t})}+2\left(\int_{\Omega}\nu_{t}(\bm{x})e_{t}^{2}(\bm{x})\mathrm{\leavevmode\nobreak\ d}\bm{x}\right)^{1/2}\left(\int_{\Omega}\nu_{t}(\bm{x})\big{|}\partial_{t}e_{t}(\bm{x})\big{|}^{2}\mathrm{\leavevmode\nobreak\ d}\bm{x}\right)^{1/2}$ $\displaystyle\leqslant$ $\displaystyle\ B^{\nu}_{0}\|e_{t}(\cdot)\|^{2}_{L^{2}(\Omega;\nu_{t})}+2\sqrt{2}\left(\widehat{B}_{0}^{2}+B_{0}^{*2}\right)^{1/2}\|e_{t}(\cdot)\|_{L^{2}(\Omega;\nu_{t})},$ $$

which follows from applying $|\partial_{t}e_{t}(\bm{x})|^{2}=|\partial_{t}u_{\theta}(\bm{x},t)-\partial_{t}u^{*}_{t}(\bm{x})|^{2}\leqslant 2\widehat{B}_{0}^{2}+2B_{0}^{*2}$.
Then plugging ([70](#A1.E70)) into ([69](#A1.E69)), we have

$$ $m_{1}\|\nabla_{\bm{x}}e_{t}(\cdot)\|^{2}_{L^{2}(\Omega;\nu_{t})}\leqslant\varepsilon\|r_{t}(\cdot)\|^{2}_{L^{2}(\Omega;\nu_{t})}+(C_{1}(\varepsilon)+B^{\nu}_{0})\|e_{t}(\cdot)\|^{2}_{L^{2}(\Omega;\nu_{t})}+C_{2}\|e_{t}(\cdot)\|_{L^{2}(\Omega;\nu_{t})}.$ (71) $$

Plugging ([50](#A1.E50)) and ([53](#A1.E53)) into ([71](#A1.E71)) gives ([51](#A1.E51)) and ([54](#A1.E54)) respectively.
∎

###### Lemma 1 (divergence theorem) .

###### Lemma 2 (Green’s first identity) .

###### Lemma 3 (Grönwall’s inequality) .

###### Proof of Lemma 3 .

###### Theorem 7 .

###### Proof of Theorem 7 .

### A.4 Proof of Theorem 3

Given the $L^{2}$ error of the score approximation, Chen et al. () provides an upper bound of KL divergence between the data distribution $\pi$ and the distribution $\widehat{\pi}_{T}$ of approximate samples drawn from the sampling dynamics ([19](#S5.E19)).
We first summarize the results from Chen et al. () in Proposition [1](#Thmproposition1).
Then we prove Theorem [3](#Thmtheorem3) based on Proposition [1](#Thmproposition1).

###### Proposition 1 (Theorem 2.5 in Chen et al. ( 2023a ) ) .

Suppose that $T\geqslant 1$, $K\geqslant 2$, and the $L^{2}$ error of the score approximation is bounded by

$$ $\sum_{k=1}^{N}h_{k}\mathbb{E}_{\bm{x}_{t_{k}}\sim\pi_{t_{k}}}\left\|\nabla_{\bm{x}}\log\pi_{t_{k}}(\bm{x}_{t_{k}})-\bm{s}_{t_{k}}(\bm{x}_{t_{k}})\right\|^{2}\leqslant T\varepsilon_{0}^{2}.$ (72) $$

Then there is a universal constant $\alpha\geqslant 2$ such that the following holds.
Under Assumption [4](#Thmassumption4), by using the exponentially decreasing (then constant) step size
$h_{k}=h\min\{\max\{t_{k},1/(4K)\},1\}$, $0<h\leqslant 1/(\alpha d)$, the sampling dynamic ([19](#S5.E19)) results in a distribution $\widehat{\pi}_{T}$ such that

$$ $\operatorname{KL}(\pi\|\widehat{\pi}_{T})\lesssim(d+M_{2})\cdot e^{-T}+T\varepsilon_{0}^{2}+d^{2}h(\log K+T),$ (73) $$

where the number of sampling steps satisfies that $N\lesssim\frac{1}{h}(\log K+T)$. Choosing $T=\log\left(\frac{M_{2}+d}{\varepsilon_{0}^{2}}\right)$ and $h=\Theta\left(\frac{\varepsilon_{0}^{2}}{d^{2}(\log K+T)}\right)$, we have $N=\mathcal{O}\left(\frac{d^{2}(\log K+T)^{2}}{\varepsilon_{0}^{2}}\right)$ and make the KL divergence $\widetilde{\mathcal{O}}\left(\varepsilon_{0}^{2}\right)$.

###### Proof of Theorem 3 .

As $\bm{s}_{t}(\bm{x})=\nabla_{\bm{x}}u_{\theta}(\bm{x},t)\cdot\mathbbm{1}\{\bm{x}\in\Omega\}$, we have

$$ $\displaystyle\sum_{k=1}^{N}h_{k}\mathbb{E}_{\bm{x}_{t_{k}}\sim\pi_{t_{k}}}\left\|\nabla_{\bm{x}}\log\pi_{t_{k}}(\bm{x}_{t_{k}})-\bm{s}_{t_{k}}(\bm{x}_{t_{k}})\right\|^{2}$ (74) $\displaystyle=$ $\displaystyle\sum_{k=1}^{N}h_{k}\int_{\Omega^{c}}\pi_{t_{k}}(\bm{x})\|\nabla_{\bm{x}}\log\pi_{t_{k}}(\bm{x})\|^{2}\mathrm{\leavevmode\nobreak\ d}\bm{x}+\sum_{k=1}^{N}h_{k}\int_{\Omega}\pi_{t_{k}}(\bm{x})\|\nabla_{\bm{x}}e_{t_{k}}(\bm{x})\|^{2}\mathrm{\leavevmode\nobreak\ d}\bm{x}$ $\displaystyle\leqslant$ $\displaystyle\sum_{k=1}^{N}h_{k}\delta+\sum_{k=1}^{N}h_{k}R_{t_{k}}\|\nabla_{\bm{x}}e_{t_{k}}(\cdot)\|^{2}_{L^{2}(\Omega;\nu_{t_{k}})}\quad\text{(using Assumption \ref{Sampler assumption: bounded region} and \ref{Sampler assumption: base target ratio})}$ $\displaystyle\leqslant$ $\displaystyle\ T\delta+\delta_{1}+C_{5}(\varepsilon)\delta_{2}+C_{2}\sqrt{\sum_{k=1}^{N}h_{k}R_{t_{k}}\delta_{2}},$ $$

where the last inequality follows from Theorem [2](#Thmtheorem2) ($m_{1}=M_{1}=1$) and

$$ $\displaystyle\sum_{k=1}^{N}h_{k}R_{t_{k}}\sqrt{\varepsilon L_{\textrm{PINN}}(t_{k};C_{1}(\varepsilon))}\leqslant$ $\displaystyle\left(\varepsilon\sum_{k=1}^{N}h_{k}R_{t_{k}}L_{\textrm{PINN}}(t_{k};C_{1}(\varepsilon))\right)^{1/2}\left(\sum_{k=1}^{N}h_{k}R_{t_{k}}\right)^{1/2}$ (75) $\displaystyle\leqslant$ $\displaystyle\ \sqrt{\sum_{k=1}^{N}h_{k}R_{t_{k}}\delta_{2}}.$ $$

Then combining ([74](#A1.E74)) and Proposition [1](#Thmproposition1) together gives the results in Theorem [3](#Thmtheorem3).
∎

###### Proposition 1 (Theorem 2.5 in Chen et al. ( 2023a ) ) .

###### Proof of Theorem 3 .

### A.5 Proof of Theorem 4

Note that $u_{t}^{*}(\bm{x})$ satisfies

$$ $\partial_{t}u^{*}_{t}(\bm{x})+\nabla_{\bm{x}}u^{*}_{t}(\bm{x})\cdot\bm{f}(\bm{x},t)+\nabla\cdot\bm{f}(\bm{x},t)-\frac{1}{2}g^{2}(t)\Delta u_{t}^{*}(\bm{x})-\frac{1}{2}g^{2}(t)\left\|\nabla_{\bm{x}}u_{t}^{*}(\bm{x})\right\|^{2}=0.$ (76) $$

And $u_{\theta}(\bm{x},t)$ satisfies

$$ $\partial_{t}u_{\theta}(\bm{x},t)+\nabla_{\bm{x}}u_{\theta}(\bm{x},t)\cdot\bm{f}(\bm{x},t)+\nabla\cdot\bm{f}(\bm{x},t)-\frac{1}{2}g^{2}(t)\Delta u_{\theta}(\bm{x},t)-\frac{1}{2}g^{2}(t)\left\|\nabla_{\bm{x}}u_{\theta}(\bm{x},t)\right\|^{2}=r_{t}(\bm{x}).$ (77) $$

Subtracting ([76](#A1.E76)) for $u^{*}$ from ([77](#A1.E77)) for $u_{\theta}$, we have

$$ $\partial_{t}e_{t}(\bm{x})+\nabla_{\bm{x}}e_{t}(\bm{x})\cdot\bm{f}(\bm{x},t)-\frac{1}{2}g^{2}(t)\left(\left\|\nabla_{\bm{x}}u_{\theta}(\bm{x},t)\right\|^{2}-\left\|\nabla_{\bm{x}}u^{*}_{t}(\bm{x})\right\|^{2}\right)-\frac{1}{2}g^{2}(t)\Delta e_{t}(\bm{x})=r_{t}(\bm{x}).$ (78) $$

Note that $\frac{1}{2}\partial_{t}e_{t}^{2}(\bm{x})=e_{t}(\bm{x})\partial_{t}e_{t}(\bm{x})$ and $\frac{1}{2}\nabla_{\bm{x}}e_{t}^{2}(\bm{x})=e_{t}(\bm{x})\nabla_{\bm{x}}e_{t}(\bm{x})$, then we obtain

$$ $\displaystyle\frac{1}{2}\partial_{t}e_{t}^{2}(\bm{x})=$ $\displaystyle\ \frac{1}{2}g^{2}(t)e_{t}(\bm{x})\left(\left\|\nabla_{\bm{x}}u_{\theta}(\bm{x},t)\right\|^{2}-\left\|\nabla_{\bm{x}}u^{*}_{t}(\bm{x})\right\|^{2}\right)+\frac{1}{2}g^{2}(t)e_{t}(\bm{x})\Delta e_{t}(\bm{x})$ (79) $\displaystyle\quad+e_{t}(\bm{x})r_{t}(\bm{x})-e_{t}(\bm{x})\nabla_{\bm{x}}e_{t}(\bm{x})\cdot\bm{f}(\bm{x},t)$ $\displaystyle=$ $\displaystyle\ \frac{1}{2}g^{2}(t)e_{t}(\bm{x})\nabla_{\bm{x}}e_{t}(\bm{x})\cdot\left(\nabla_{\bm{x}}u_{\theta}(\bm{x},t)+\nabla_{\bm{x}}u^{*}_{t}(\bm{x})\right)+\frac{1}{2}g^{2}(t)e_{t}(\bm{x})\Delta e_{t}(\bm{x})$ $\displaystyle\quad+e_{t}(\bm{x})r_{t}(\bm{x})-e_{t}(\bm{x})\nabla_{\bm{x}}e_{t}(\bm{x})\cdot\bm{f}(\bm{x},t)$ $\displaystyle=$ $\displaystyle\ \frac{1}{4}g^{2}(t)\nabla_{\bm{x}}e_{t}^{2}(\bm{x})\cdot\left(\nabla_{\bm{x}}u_{\theta}(\bm{x},t)+\nabla_{\bm{x}}u^{*}_{t}(\bm{x})\right)+\frac{1}{2}g^{2}(t)e_{t}(\bm{x})\Delta e_{t}(\bm{x})$ $\displaystyle\quad+e_{t}(\bm{x})r_{t}(\bm{x})-\frac{1}{2}\nabla_{\bm{x}}e_{t}^{2}(\bm{x})\cdot f(\bm{x},t).$ $$

We integrate ([79](#A1.E79)) to get an equation for $\|e_{t}(\cdot)\|^{2}_{L^{2}(\Omega)}$ given by

$$ $\displaystyle\partial_{t}\left\|e_{t}(\cdot)\right\|^{2}_{L^{2}(\Omega)}=$ $\displaystyle\ \frac{1}{2}g^{2}(t)\int_{\Omega}\nabla_{\bm{x}}e_{t}^{2}(\bm{x})\cdot\left(\nabla_{\bm{x}}u_{\theta}(\bm{x},t)+\nabla_{\bm{x}}u^{*}_{t}(\bm{x})\right){\textnormal{d}}\bm{x}+\ g^{2}(t)\int_{\Omega}e_{t}(\bm{x})\Delta e_{t}(\bm{x}){\textnormal{d}}\bm{x}$ (80) $\displaystyle\quad+2\int_{\Omega}e_{t}(\bm{x})r_{t}(\bm{x}){\textnormal{d}}\bm{x}-\int_{\Omega}\nabla_{\bm{x}}e_{t}^{2}(\bm{x})\cdot\bm{f}(\bm{x},t){\textnormal{d}}\bm{x}$ $\displaystyle=$ $\displaystyle\ -\frac{1}{2}g^{2}(t)\int_{\Omega}e_{t}^{2}(\bm{x})\cdot\left(\Delta u_{\theta}(\bm{x},t)+\Delta u^{*}_{t}(\bm{x})\right){\textnormal{d}}\bm{x}-g^{2}(t)\int_{\Omega}\left\|\nabla_{\bm{x}}e_{t}(\bm{x})\right\|^{2}{\textnormal{d}}\bm{x}$ $\displaystyle\quad+2\int_{\Omega}e_{t}(\bm{x})r_{t}(\bm{x}){\textnormal{d}}\bm{x}+\int_{\Omega}e_{t}^{2}(\bm{x})\cdot\left[\nabla\cdot\bm{f}(\bm{x},t)\right]{\textnormal{d}}\bm{x}$ $\displaystyle\leqslant$ $\displaystyle\ -\frac{m_{1}}{2}\left(B_{2}^{*}+\widehat{B}_{2}\right)\left\|e_{t}(\cdot)\right\|^{2}_{L^{2}(\Omega)}-m_{1}\left\|\nabla_{\bm{x}}e_{t}(\cdot)\right\|_{L^{2}(\Omega)}^{2}+{\varepsilon}\left\|r_{t}(\cdot)\right\|^{2}_{L^{2}(\Omega)}$ $\displaystyle\ \quad+\frac{1}{{\varepsilon}}\left\|e_{t}(\cdot)\right\|^{2}_{L^{2}(\Omega)}+\mathrm{m}_{2}\left\|e_{t}(\cdot)\right\|^{2}_{L^{2}(\Omega)}$ $\displaystyle=$ $\displaystyle\ C_{1}^{\textrm{U}}({\varepsilon})\left\|e_{t}(\cdot)\right\|_{L^{2}(\Omega)}^{2}+{\varepsilon}\left\|r_{t}(\cdot)\right\|_{L^{2}(\Omega)}^{2}-m_{1}\left\|\nabla_{\bm{x}}e_{t}(\cdot)\right\|^{2}_{L^{2}(\Omega)}$ $\displaystyle\leqslant$ $\displaystyle\ C_{1}^{\textrm{U}}({\varepsilon})\left\|e_{t}(\cdot)\right\|_{L^{2}(\Omega)}^{2}+{\varepsilon}\left\|r_{t}(\cdot)\right\|_{L^{2}(\Omega)}^{2}.$ $$

Note that $e_{0}(\bm{x})=0$ for any $\bm{x}\in\Omega$, then using the Grönwall inequality, we have for any $t\in[0,T]$,

$$ $\|e_{t}(\cdot)\|^{2}_{L^{2}(\Omega)}\leqslant{\varepsilon}\int_{0}^{t}e^{C_{1}^{\textrm{U}}({\varepsilon})(t-s)}\left\|r_{s}(\cdot)\right\|^{2}_{L^{2}(\Omega)}{\textnormal{d}}s:={\varepsilon}L_{\textrm{PINN}}^{\textrm{Unif}}(t;C_{1}^{\textrm{U}}({\varepsilon})).$ (81) $$

Note that from ([80](#A1.E80)),

$$ $m_{1}\|\nabla_{\bm{x}}e_{t}(\cdot)\|^{2}_{L^{2}(\Omega)}\leqslant{\varepsilon}\|r_{t}(\cdot)\|^{2}_{L^{2}(\Omega)}+C_{1}^{\textrm{U}}({\varepsilon})\|e_{t}(\cdot)\|^{2}_{L^{2}(\Omega)}-\partial_{t}\|e_{t}(\cdot)\|^{2}_{L^{2}(\Omega)}.$ (82) $$

We can bound $\partial_{t}\|e_{t}(\cdot)\|^{2}_{L^{2}(\Omega)}$ as follows

$$ $\displaystyle\left|\partial_{t}\|e_{t}(\cdot)\|^{2}_{L^{2}(\Omega)}\right|=\left|\partial_{t}\left(\int_{\Omega}e_{t}^{2}(\bm{x}){\textnormal{d}}\bm{x}\right)\right|=2\left|\int_{\Omega}e_{t}(\bm{x})\partial_{t}e_{t}(\bm{x}){\textnormal{d}}\bm{x}\right|$ (83) $\displaystyle\leqslant$ $\displaystyle\ 2\left(\int_{\Omega}e_{t}^{2}(\bm{x}){\textnormal{d}}\bm{x}\right)^{1/2}\left(\int_{\Omega}\left|\partial_{t}e_{t}(\bm{x})\right|^{2}{\textnormal{d}}\bm{x}\right)^{1/2}\leqslant 2\sqrt{2}\left(\widehat{B}_{0}^{2}+B_{0}^{*2}\right)^{1/2}\|e_{t}(\cdot)\|_{L^{2}(\Omega)},$ $$

which follows from applying $|\partial_{t}e_{t}(\bm{x})|^{2}=|\partial_{t}u_{\theta}(\bm{x},t)-\partial_{t}u^{*}_{t}(\bm{x})|^{2}\leqslant 2\widehat{B}_{0}^{2}+2B_{0}^{*2}$.
Then plugging ([83](#A1.E83)) into ([82](#A1.E82)), we have

$$ $m_{1}\|\nabla_{\bm{x}}e_{t}(\cdot)\|^{2}_{L^{2}(\Omega)}\leqslant{\varepsilon}\|r_{t}(\cdot)\|^{2}_{L^{2}(\Omega)}+C_{1}^{\textrm{U}}({\varepsilon})\|e_{t}(\cdot)\|^{2}_{L^{2}(\Omega)}+C_{2}\|e_{t}(\cdot)\|_{L^{2}(\Omega)}.$ (84) $$

Plugging ([82](#A1.E82)) into ([84](#A1.E84)), we complete the proof.

### A.6 Proof of Theorem 5

As $\bm{s}_{t}(\bm{x})=\nabla_{\bm{x}}u_{\theta}(\bm{x},t)\cdot\bm{1}\{\bm{x}\in\Omega\}$, we have

$$ $\displaystyle\sum_{k=1}^{N}h_{k}\mathbb{E}_{\bm{x}_{t_{k}}\sim\pi_{t_{k}}}\left\|\nabla_{\bm{x}}\log\pi_{t_{k}}(\bm{x}_{t_{k}})-\bm{s}_{t_{k}}(\bm{x}_{t_{k}})\right\|^{2}$ (85) $\displaystyle=$ $\displaystyle\sum_{k=1}^{N}h_{k}\int_{\Omega^{c}}\pi_{t_{k}}(\bm{x})\|\nabla_{\bm{x}}\log\pi_{t_{k}}(\bm{x})\|^{2}{\textnormal{d}}\bm{x}+\sum_{k=1}^{N}h_{k}\int_{\Omega}\pi_{t_{k}}(\bm{x})\|\nabla_{\bm{x}}e_{t_{k}}(\bm{x})\|^{2}{\textnormal{d}}\bm{x}$ $\displaystyle\leqslant$ $\displaystyle\sum_{k=1}^{N}h_{k}\delta+\sum_{k=1}^{N}h_{k}\max_{\bm{x}\in\Omega}\left\{\pi_{t_{k}}(\bm{x})\right\}\cdot\|\nabla_{\bm{x}}e_{t_{k}}(\cdot)\|^{2}_{L^{2}(\Omega)}$ $\displaystyle\leqslant$ $\displaystyle\ T\delta+\delta_{1}\cdot\textrm{Vol}(\Omega)+C_{1}^{\textrm{U}}({\varepsilon})\delta_{2}\cdot\textrm{Vol}(\Omega)+C_{2}\sqrt{\sum_{k=1}^{N}h_{k}\max_{\bm{x}\in\Omega}\left\{\pi_{t_{k}}(\bm{x})\right\}\cdot\delta_{2}\cdot\textrm{Vol}(\Omega)},$ $$

where the last inequality follows from the result in Theorem [4](#Thmtheorem4) and

$$ $\displaystyle\ \sum_{k=1}^{N}h_{k}\max_{\bm{x}\in\Omega}\left\{\pi_{t_{k}}(\bm{x})\right\}\cdot\sqrt{{\varepsilon}L_{\textrm{PINN}}^{\textrm{Unif}}(t_{k};C_{1}^{\textrm{U}}({\varepsilon}))}$ (86) $\displaystyle\leqslant$ $\displaystyle\ \left({\varepsilon}\sum_{k=1}^{N}h_{k}\max_{\bm{x}\in\Omega}\left\{\pi_{t_{k}}(\bm{x})\right\}\cdot L_{\textrm{PINN}}^{\textrm{Unif}}(t_{k};C_{1}^{\textrm{U}}({\varepsilon}))\right)^{1/2}\left(\sum_{k=1}^{N}h_{k}\max_{\bm{x}\in\Omega}\left\{\pi_{t_{k}}(\bm{x})\right\}\right)^{1/2}$ $\displaystyle\leqslant$ $\displaystyle\ \sqrt{\sum_{k=1}^{N}h_{k}\max_{\bm{x}\in\Omega}\left\{\pi_{t_{k}}(\bm{x})\right\}\cdot\delta_{2}\cdot\textrm{Vol}(\Omega)}.$ $$

Combining ([85](#A1.E85)) and Proposition [1](#Thmproposition1), we complete our proof.

## Appendix B Limitations

As we use LMC for collocation generation in DPS, there is a risk of missing modes if short LMC runs do not adequately cover the high-density domain.
In such cases, running LMC for an annealed path of target distributions or adopting the adversarial training method in Wang et al. () for collocation points maybe helpful.
Also, solving high dimensional PDEs via PINN can be challenging, and we may use techniques such as stochastic dimension gradient descent or the Hutchinson trick to scale DPS to high dimensional problems (Hu et al., , ).

## Appendix C Additional Experimental Details and Results

### C.1 Baselines

We benchmark DPS performance against a wide range of strong baseline methods.
For MCMC methods, we consider the Langevin Monte Carlo (LMC).
For LMC, we run 100,000 iterations with step sizes 0.02, 0.002, 0.0002.
Then we choose the samples with the best performance.
As for sampling methods using reverse diffusion, we include RDMC (Huang et al., ), and SLIPS (Grenioux et al., ).
We use the implementation of SLIPS and RDMC from Grenioux et al. () and choose Geom($1,1$) as the SL scheme for SLIPS.
For each algorithm, we search its hyper-parameters within a predetermined grid, similar to Grenioux et al. ().
We also compare with VI-based PIS (Zhang & Chen, ) and DIS (Berner et al., ).
We use the implementation of PIS and DIS from Berner et al. ().

### C.2 Targets

9-Gaussians is a 2-dimensional Mixture of Gaussians where there are 9 modes designed to be well-separated from each other.
The modes share the same variance of 0.3 and the means are located in the grid of $\{-5,0,5\}\times\{-5,0,5\}$.
We set challenging mixing proportions between different modes as shown in Table [3](#A3.T3).

Rings is the inverse polar reparameterization of a $2$-dimensional distribution $p_{z}$ which has itself a decomposition
into two univariate marginals $p_{r}$ and $p_{\theta}$: $p_{r}$ is a mixture of $4$ Gaussian distributions $\mathcal{N}(i,0.2^{2})$ with $i=2,4,6,8$ describing the radial position and $p_{\theta}$ is a uniform distribution over $[0,2\pi)$, which describes the angular position of the samples.
We also set challenging mixing proportions between different modes of $p_{r}$ as shown in Table [4](#A3.T4).

**Table 3: Mixing proportions between 9 modes in 9-Gaussians.**
| Modes | $(-5,-5)^{\prime}$ | $(-5,0)^{\prime}$ | $(-5,5)^{\prime}$ | $(0,-5)^{\prime}$ | $(0,0)^{\prime}$ | $(0,5)^{\prime}$ | $(5,-5)^{\prime}$ | $(5,0)^{\prime}$ | $(5,5)^{\prime}$ |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Weight | $0.2$ | $0.04$ | $0.2$ | $0.04$ | $0.04$ | $0.04$ | $0.2$ | $0.04$ | $0.2$ |

**Table 4: Mixing proportions between 4 modes in rings.**
| Modes | $r=2$ | $r=4$ | $r=6$ | $r=8$ |
| --- | --- | --- | --- | --- |
| Weight | $0.05$ | $0.45$ | $0.05$ | $0.45$ |

Funnel is a classical sampling benchmark problem from Neal (); Hoffman et al. ().
This 10-dimensional density is defined by

$$ $\mu(\bm{x}):=\mathcal{N}(x_{0};0,9)\mathcal{N}(\bm{x}_{1:9};\bm{0},\exp(x_{0})\bm{I}_{9}).$ $$

Double-well is a high-dimensional distribution which share the unnormalized density:

$$ $\mu(\bm{x}):=\exp\left(\sum_{i=0}^{w-1}-x_{i}^{4}+6x_{i}^{2}+0.5x_{i}-\sum_{i=w}^{d-1}0.5x_{i}^{2}\right).$ $$

We choose $w=3$ and $d=30$, leading to a $30$-dimensional distribution contained $8$ modes with challenging mixing proportions between different modes.

### C.3 Diffusion-PINN Sampler

##### Model.

The model architecture of $\mathrm{NN}_{\theta}(\bm{x},t):\mathbb{R}^{d}\times[0,T]\to\mathbb{R}$ in $u_{\theta}(\bm{x},t)$ is

$$ $\mathrm{NN}_{\theta}(\bm{x},t)=\mathrm{MLP}^{\textrm{dec}}\left(\mathrm{MLP}^{\text{embx}}(\bm{x})+\mathrm{MLP}^{\text{embt}}(\mathrm{emb}(t))\right),$ $$

where $\mathrm{MLP}^{\textrm{dec}}$ represents a decoder implemented as MLPs with layer widths $[128,128,128,1]$. The component $\mathrm{MLP}^{\text{embx}}$ serves as a data embedding block and is implemented as MLPs with layer widths $[2,128]$.
$\mathrm{MLP}^{\text{embt}}$ functions as a time embedding block, implemented as MLPs with layer widths $[256,128,128]$. The input to $\mathrm{MLP}^{\text{embt}}$ is derived from the sinusoidal positional embedding (Vaswani et al., ) of $t$.
All these three MLPs utilize the GELU activation function.

##### Training.

In our implementation, we choose $\bm{f}(\bm{x},t)=-\frac{\bm{x}}{2(1-t)}$ and $g(t)=\sqrt{\frac{1}{1-t}}$, leading to the following forward process

$$ $\mathrm{\leavevmode\nobreak\ d}\bm{x}_{t}=-\frac{\bm{x}_{t}}{2(1-t)}\mathrm{\leavevmode\nobreak\ d}t+\sqrt{\frac{1}{1-t}}\mathrm{\leavevmode\nobreak\ d}\bm{B}_{t},\quad\bm{x}_{0}\sim\pi,\quad T_{\textrm{min}}\leqslant t\leqslant T_{\textrm{max}}.$ (87) $$

This admits the explicit conditional distribution $\pi_{t\mid 0}(\bm{x}_{t}|\bm{x}_{0})=\mathcal{N}(\bm{x}_{t};\sqrt{1-t}\cdot\bm{x}_{0},t\cdot\bm{I}_{d})$.
We choose $T_{\textrm{min}}=0.001$ and $T_{\textrm{max}}=0.999$ in practice.
The corresponding log-density FPE becomes

$$ $\partial_{t}u_{t}(\bm{x})=\frac{1}{2(1-t)}\left[\Delta u_{t}(\bm{x})+\left\|\nabla_{\bm{x}}u_{t}(\bm{x})\right\|^{2}+\bm{x}\cdot\nabla_{\bm{x}}u_{t}(\bm{x})+d\right]:=\frac{1}{2(1-t)}\mathcal{L}_{\textrm{L-FPE}}^{\textrm{prac}}u_{t}(\bm{x}).$ (88) $$

We choose $\beta(t)=2(1-t)$ to make training more stable, leading to the following training objective

$$ $\displaystyle L_{\textrm{train}}^{\textrm{prac}}(u_{\theta}):=$ $\displaystyle\mathbb{E}_{t\sim\mathcal{U}[T_{\textrm{min}},T_{\textrm{max}}]}\mathbb{E}_{\bm{x}_{t}\sim\nu_{t}}\left[\left\|2(1-t)\cdot\partial_{t}u_{\theta}(\bm{x}_{t},t)-\mathcal{L}_{\textrm{L-FPE}}^{\textrm{prac}}u_{\theta}(\bm{x}_{t},t)\right\|^{2}\right]$ (89) $\displaystyle\quad+\lambda\cdot\mathbb{E}_{\bm{z}\sim\mathcal{N}(\bm{0},\bm{I}_{d})}\left[\left\|\nabla_{\bm{z}}u_{\theta}(\bm{z},T_{\textrm{max}})+\bm{z}\right\|^{2}\right],$ $$

where $\lambda$ is the regularization coefficient.
It is enough for us to use PINN residual loss without regularization except for Funnel where the regularization is quite useful and we use $\lambda=1$.
To generate collocation points for PINN, we run a short chain of LMC with a large step size.
The hyper-parameters used in LMC for different targets are reported in Table [5](#A3.T5).
We generate fresh collocation points per iteration except for Funnel where we resample new collocation points per $10,000$ iterations.

**Table 5: Hyper-parameters used in LMC for generating collocation points.**
|  | 9-Gaussians | Rings | Funnel | Double-well |
| --- | --- | --- | --- | --- |
| step size | $1.0$ | $0.15$ | $0.02$ | $0.02$ |
| iterations | $60$ | $100$ | $10,000$ | $100$ |
| batch size | $128$ | $200$ | $200$ | $700$ |
| refresh samples per iteration | ✓ | ✓ | ✗ | ✓ |

We train all models with Adam optimizer (Kingma & Ba, ).
The hyper-parameters used in training are summarized in Table [6](#A3.T6).
We use a linear decay schedule for the learning rate in all experiments.

**Table 6: Hyper-parameters for training PINN.**
|  | 9-Gaussians | Rings | Funnel | Double-well |
| --- | --- | --- | --- | --- |
| learning rate | $0.0005$ | $0.0005$ | $0.0001$ | $0.0005$ |
| max norm of gradient clipping | $1.0$ | $1.0$ | $1000.0$ | $1.0$ |
| regularization coefficient $\lambda$ | 0 0 | 0 0 | $1$ | 0 0 |
| total training iterations | $400$k | $1,000$k | $800$k | $1,500$k |

##### Sampling.

The corresponding reverse process is given by

$$ $\mathrm{\leavevmode\nobreak\ d}\bm{x}_{t}=\left(\frac{\bm{x}_{t}}{2t}+\frac{\nabla_{\bm{x}}\log\pi_{1-t}(\bm{x}_{t})}{t}\right)\mathrm{\leavevmode\nobreak\ d}t+\sqrt{\frac{1}{t}}\mathrm{\leavevmode\nobreak\ d}\bm{B}_{t},\quad\bm{x}_{0}\sim\pi_{T_{\textrm{max}}},\quad T_{\textrm{min}}\leqslant t\leqslant T_{\textrm{max}}.$ (90) $$

To simulate ([90](#A3.E90)), we approximate $\pi_{T_{\textrm{max}}}\approx\mathcal{N}(\bm{0},\bm{I}_{d})$ and use the exponential integrator scheme with the score approximation $\bm{s}_{t}(\bm{x})\approx\nabla_{\bm{x}}\log\pi_{t}(\bm{x})$.
In practice, we use $\bm{s}_{t}(\bm{x}):=\nabla_{\bm{x}}u_{\theta}(\bm{x},t)\cdot\mathbbm{1}\{\bm{x}\in\Omega\}$ where $u_{\theta}(\bm{x},t)$ is the approximated log-density provided by the PINN approach (trained by Algorithm [1](#alg1)) and $\Omega$ is a chosen bounded region that covers the high density domain of $\pi_{t}$ for any $t\in[T_{\textrm{min}},T_{\textrm{max}}]$.
We use $\Omega:=\{\bm{x}\in\mathbb{R}^{d}:\|\bm{x}\|\leqslant R\}$ in all experiments, the choice of $R$ is reported in Table [7](#A3.T7).
Our sampling process is summarized in Algorithm [2](#alg2).
We provide more sampling performances of different methods for different targets in Figure [6](#A3.F6) and sample trajectories from DPS in Figure [7](#A3.F7).

**Table 7: The diameter of the truncated region for different targets.**
|  | 9-Gaussians | Rings | Funnel | Double-well |
| --- | --- | --- | --- | --- |
| $R$ | $20$ | $20$ | $2000$ | $30$ |

Figure: Algorithm 2 : Sampling from reverse process

Figure: Figure 6: Sampling performance of different methods for $9$-Gaussians ($d=2$), Rings ($d=2$), Funnel ($d=10$) and Double-well ($d=30$).
Refer to caption: /html/2410.15336/assets/x6.png

Figure: Figure 7: Sample trajectories from DPS for $9$-Gaussians ($d=2$), Rings ($d=2$), Funnel ($d=10$) and Double-well ($d=30$).
Refer to caption: /html/2410.15336/assets/x7.png

## Peer Reviews

confidence: 5
contribution: 2
initial_rating: 3
presentation: 2
questions: - $\ell_reg$ seems ad-hoc. Why regress onto the boundary condition of the score instead of (11b), which seems more aligned with PINN loss?
- Given that x0 comes from LMC already, is the method computationally cheaper compared to standard MCMC methods?
soundness: 2
strengths: - The paper is generally easier to follow. Theoretical results in Sec 5 could be of interested for readers from other domains.
summary: This paper explores diffusion-based models for unnormalized sampling and introduces the Diffusion-PINN Sampler (DPS). Using physics-informed neural networks (PINNs), DPS directly approximates the log-density of SDE marginals, enabling more precise modeling of complex distributions. The authors provide theoretical convergence analysis and validate the method on several synthetic datasets.
weaknesses: - The proposed method is computationally expansive compared to other diffusion-based models to train—due to the Laplacian in PINN loss—and to sample from—due to the evaluation of gradient at every time step.
- Insufficient comparison to modern baselines — in addition to  MC methods like HMC or SMC, there’re many other diffusion/SDE based methods, such as [1,2], just to name a few.
- Experiments were only conducted on rather simple, synthetic, target. I’ll be more convinced to see some higher-dim experiments and/or real-world dataset. 
- Parametrize NN with log mu (12) seems like a strong inductive bias and can be infeasible for many practical applications (e.g., sample Boltzmann distribution for conformation generation) when querying energy functions are expansive. Can the authors provide ablation study without such parametrization? 
- Related works Section is rather short and should be extended.

[1] Particle Denoising Diffusion Sampler (ICML 2024)
[2] Improved sampling via learned diffusions (ICLR 2024)

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 4
questions: See weaknesses.
soundness: 3
strengths: The paper is very well written! It was pleasant to read and easy to follow.

- The authors provide observations of failure of score FPE and propose a motivated and illustrated solution : solving the log-density FPE and plugging it in the reverse process afterwards.

- The authors provide a good amount of experiments and parsimonious illustrations (just what is needed to illustrate their claims).

- The authors assess performances of their methods with good metrics. It could seem naive to say, but a lot of people miss the use of probabilistic measures for probabilistic problems. The choice of KL/Fisher divergence is welcomed along with the $L^2$ of the log-density.

- The use of PINN for such a problem is relevant. It bases the method on well-established (and timely) fields which is nice!

- The overall structure of the work is clean and each part is self-sufficient.

- Analysis, performed on toy problems at first, allows for good interpretation of the results.

- The proposed sampling methods show big improvement over other baselines for considered problems which further motivates it.
summary: The paper presents a new approach for diffusion sampling. Inspired by existing methods that use physics-informed neural networks, they address the problem of solving the reverse diffusion process using the Fokker Plank equation which models the evolution of marginals $p_t(x)$ in the forward process. Specifically, they identify flaws of existing methods and address them. The first one is the failure of the score FPE and the second is the lack of theoretical analysis of proposed methods. They end up with a method called DPS (Diffusion-PINN sampler) that solves for log-density FPE. They show promising results and consistent improvement over other methods, providing, in addition, convergence analysis of the method.
weaknesses: I'm really sorry, but I think you will have to change the name (at least the reduced name) of your method! DPS is already a well known and well established sampling method for posterior problems [1]. I think it would be beneficial to avoid being shadowed by existing method. Moreover, it's more than just a sampling method that you present. It's a conjugate training and sampling algorithm. Your main change is about what your network is trained for (the log density) and not about the FPE which is not new in itself (Although the training and the sampling make a whole). 

Otherwise, I point below a few questions and small concerns.

**Small details**

- Line 60: *to ensure convergence guarantee*, it's nitpicking but isn't it a bit redundant?

**Performance concerns**

- I'm a little bit concerned about the performances of the method both in time and resources. It would be nice to, at least discuss them and at best provide some measures in Appendix for example. I detail below my concerns.

- You model the log-density, then you have to take its gradient in order to compute the score that you plug into the reverse process, right? What is the cost of this compared to direct score matching? Indeed, I presume that taking the grad at every step is not without cost. 

- The PINN loss involves gradient and Laplacian operators. What is the cost of computing such a loss compared to simple denoising score matching (which, at the end, is not more complex than an mse)? Does it slow down training? Do you use optimized tricks like jvp/vjp when possible? Could you comment a bit on that please?

- In table 1, isn't it misleading to report KL only on first dimensions for Funnel and Double-well? Can't you report either on several sub-combinations of dimensions or use another measure? Is it due to the cost of the KL in higher dimension? What about more tractable divergence metrics such as the Sinkhorn divergence?

**Sampling method**

- From your text and as the first requirement in **Algorithm 1** you speak about access to the initial condition which derives from the unnormalized density. It is still not clear to me how you access to such density? 

- You present a sampling method, and even more, for diffusion models. It would be nice to have a higher dimensional problem (I agree that it is not necessary to cover your claims here). This method should be studied in more realistic settings also (with images, high dimensional data, ...) and I'm worried about the cost of computing explicitly the score from the log-density or the cost of the loss in such settings. At least comment on that in the main paper.

- Echoing with the previous comment: What about conditional sampling? Posterior problems. Technically your method can extend easily. You compute the prior score as for now and plug the likelihood to guide the sampling. It would be nice (and I think easy) to show, in Appendix if you lack space, the conditional generation for one toy problem. For example, generating only on one mode in the 9-Gaussians or Rings problem. I'm curious to see the benefit of more accurate prior score to such problems. I'm also curious to see if you still perform much better than other considered methods. Indeed, the likelihood guiding score can sometimes help and simplify a bit the problem.

- Nitpick: In figure 3 you say "Sampling performance". Those figures do not show any performances. Maybe change it to "Samples from different methods ...".


**Additional suggestion**

Figure 5 (center and right) are not really revealing. I know that it's hard to display samples in high dim. Could you put corner plots instead? Even if it's on a subset of dimension of the problem. It would be easier to see the differences with and without regularization if you display marginals along this 2-dimensional joint.


My current score is 6. I would gladly increase my score once my concerns have been addressed and discussed. 
Thanks for your good work!

**[1] Diffusion Posterior Sampling for General Noisy Inverse Problems - ICLR 2023 spotlight**

---

confidence: 4
contribution: 2
initial_rating: 3
presentation: 2
questions: 1. If I understand it correctly, you claim that Example 1 shows that the Fisher divergence can be arbitrarily small, while the KL divergence stays large at the same time. However, clearly, if the Fisher divergence goes to zero, so does the KL divergence (by the definition of a divergence). Therefore, your statement seems to mean that "the KL divergence can be (much) larger than the Fisher divergence", is this correct? This implies the following questions:
    - a) Since the "scale" of any divergence is rather arbitrary, why should the scale of the KL divergence be better suited than the scale of the Fisher divergence? (That said, I think it could be helpful to consider log-scales in Figure 1.)
    - b) Does your observation generalize to distributions that are not Gaussian?
    - c) Are you aware of any bounds that directly relate the KL divergence to the Fisher divergence? (I’m thinking of something similar to the Poincare inequality.)
2. I find the numerical experiments in Figure 4 more enlightening than the statement in Example 1. Do those observations carry over to other (more complicated, high-dimensional) examples? It seems to me that Lai et al. (2023) have found that using the score FPE (as a regularizer) seems to work well. Can you comment on this finding, given your argument against the score FPE? Also, what happens if you use more gradient steps in the example from Figure 4 (both losses seem to not have converged yet)? Does the score FPE eventually converge sufficiently? (Actually, how do you define "score error"? Is it some kind of $L^2$ error?)
3. If one approximates the log-density as you do, one later must take gradients to get a score approximation, which might incur errors. Can you comment on potential numerical implications?
4. Did you use Hutchinson’s gradient estimator that you introduce in lines 262 and following? While it removes a derivative, it is known to increase the variance of gradient estimators. 
5. In (20) one considers $L_\mathrm{PINN}(t_k; C_1(\varepsilon))$ for multiple $t_k$, however, the integral in $L_\mathrm{PINN}$ always starts at $t=0$. Does this mean that the larger $N$ the larger is the sum? (I would not have expected the term to increase with finer time discretization.)

Some comments on notation and typos:
1. The notation seems to be not consistent all the time, e.g. regarding the time index: $g(t), f(x,t)$, but $p_t(x)$ and not $p(x, t)$; also: $u_t(x)$ vs. $u_\theta(x, t)$.
2. Why do you write $\nabla_x$, but not $\Delta_x$ or $\nabla_x \cdot f$? Also, you write $\nabla_x u_\theta(z, T)$ (e.g. in l. 254). Shouldn’t it be $\nabla_z u_\theta(z, T)$? (Alternatively, you should explain your notation. Suggestion: Writing $\nabla$ instead of $\nabla_x$ solves the issue.)
3. Some citations are not correctly formatted (e.g. "monte carlo" instead of "Monte Carlo"), some versions are old.
4. Inconsistent writing: Fokker-Planck vs. Fokker Planck.
5. 83: Definition of norm: "+" is missing.
6. 83: "Let ν denotes".
7. 199: "until the every end".
8. 320: Assumption -> Assumptions (this typo appear multiple times throughout the paper).
9. 326: After (17) replace "." with ",".
10. 333: (Arjovsky et al. (2017)).
11. 355: "Let $\hat{\pi}_T$ denotes".
12. In the appendix, there are often commas or periods missing after equations, see, e.g. (22), (23) etc.
13. 750: "Let $p(x)$ denotes".
soundness: 3
strengths: Sampling via learned diffusion processes is an interesting and active field of research, which is arguably more challenging than the typical generative modeling task since no samples from the target are available. The connection to the underlying PDEs is interesting, however, not novel (see weaknesses). Still, the presentation is sound and offers a detailed and interesting theoretical analysis.

The paper is mathematically well written, however, some (practical) implications and motivations could be made clearer for the convenience of the reader.
summary: The paper studies sampling via time-reversing a diffusion process that starts in the target distribution, i.e. it considers diffusion-based generative modeling, where, however, no data samples from the target distribution are available, but the (unnormalized) target density can be evaluated. The authors suggest to learn the score function that is necessary to sample accurately by approximating the underlying Fokker-Planck PDE (in log-space) with Physics-informed neural networks (PINNs). In particular, it is argued that solving the Fokker-Planck PDE and taking derivatives afterwards is advantageous over solving a corresponding PDE for the score directly. The paper presents theoretical analyses that bound the error of the approximations of the log-density and the score, respectively, by the (weighted) PINN objective. Leveraging previous work, this in consequence allows to bound the sampling error by the PINN objective. Some numerical evaluation on synthetic sampling problems in moderate dimensions is provided, demonstrating a proof of concept of the method.
weaknesses: In my opinion, the paper has the following main weaknesses.

1. **Novelty.** The idea of employing PINNs for diffusion-based sampling has already been suggested in [1], [2], [3] and [4]. While [1] and [2] have been cited, only a preliminary version of [3] seems to be mentioned. In fact, in [3] (almost) the same algorithm has been presented (without relying on sampling from an MCMC algorithm) and multiple numerical experiments (including also variants of the suggested method, such as e.g. deterministic evolutions, general evolutions, annealed/prescribed evolutions) have been conducted.

2. **Numerical evaluation.** While the numerical evaluation in the paper showcases good performance of the suggested algorithm, the considered examples are rather synthetic and of moderate dimension. The question of how the algorithm scales to real-world problems and higher-dimensional examples (say, $d > 100$, which might be challenging for PINNs due to computation of higher-order derivatives) remains open.

3. **Practical relevance of theoretical results.** While the theoretical analyses look sound and are nice from a mathematical perspective, I am not certain about their practical relevance. It is clear that in principle a zero PINN loss implies zero score error and (up to discretization and up to the prior error) perfect sampling. At the same time, the provided bound seems to not allow for quantitative statements since constants cannot be computed and it is (due to the Grönwall inequality) probably far from being tight. Furthermore, it seems that the theoretical results heavily rely on earlier works by [5] and [6].

Some additional comments are the following:

4. It is unclear to me how generating collocation points with MCMC-like algorithms (e.g. LMC) is sufficient. It is known that those algorithms converge exponentially slowly for multimodal targets where the modes are sufficiently separated. Therefore, some exploration strategy seems to be necessary in order to get collocation points that cover regions containing distant modes.
5. You write that trajectory-based alternatives lead to "computational complexity associated with differentiating through SDE solvers". However, this is not necessarily true since off-policy training can be chosen also for trajectory-based methods, e.g. by using divergences other than the KL-divergence, see, e.g., [7], [8].
6. The related work section barely lists any of the many alternative methods for diffusion-based sampling that have been suggested in the last 1-2 years.
7. Concurrently to the work of [9], [1] has also derived the PDE for the log-density (and in fact, also approached it via PINNs, see above).
8. You write about "a linear interpolation (i.e., annealing) path between the target distribution and a simple prior", however it’s linear only in log-space.
9. It would be helpful to define the Fisher divergence in the main text (and not only in Appendix A.2).
10. Sampling quality can usually not be measured well by only one metric and therefore it would be even better to have more metrics (e.g. Wasserstein, ELBO, ESS, comparison to expectation reference values etc.).

I apologize in case I misunderstood certain aspects of your paper and am looking forward to potential corrections.

[1] Julius Berner, Lorenz Richter, and Karen Ullrich. An optimal control perspective on diffusion-based generative modeling. Transactions on Machine Learning Research, 2024.

[2] Bálint Máté and François Fleuret. Learning interpolations between Boltzmann densities. Transactions on Machine Learning Research, 2023.

[3] Sun, Jingtong, et al. "Dynamical measure transport and neural PDE solvers for sampling." arXiv preprint arXiv:2407.07873 (2024).

[4] Albergo, Michael S., and Eric Vanden-Eijnden. "NETS: A Non-Equilibrium Transport Sampler." arXiv preprint arXiv:2410.02711 (2024).

[5] Sitan Chen, Sinho Chewi, Jerry Li, Yuanzhi Li, Adil Salim, and Anru Zhang. Sampling is as easy as learning the score: theory for diffusion models with minimal data assumptions. In The Eleventh International Conference on Learning Representations, 2023b. 

[6] Teo Deveney, Jan Stanczuk, Lisa Maria Kreusser, Chris Budd, and Carola-Bibiane Schönlieb. Closing the ode-sde gap in score-based diffusion models through the fokker-planck equation. arXiv preprint arXiv:2311.15996, 2023. 

[7] Lorenz Richter and Julius Berner. Improved sampling via learned diffusions. In International Conference on Learning Representations, 2024.

[8] Sendera, Marcin, et al. "On diffusion models for amortized inference: Benchmarking and improving stochastic control and sampling." arXiv preprint arXiv:2402.05098 (2024).

[9] Chieh-Hsin Lai, Yuhta Takida, Naoki Murata, Toshimitsu Uesaka, Yuki Mitsufuji, and Stefano Ermon. Fp-diffusion: Improving score-based diffusion models by enforcing the underlying score fokker-planck equation. In International Conference on Machine Learning, pp. 18365–18398. PMLR, 2023.

---

confidence: 4
contribution: 2
initial_rating: 5
presentation: 3
questions: see Weaknesses.
soundness: 3
strengths: 1. The paper is well-composed and the concepts are clearly articulated, making it easy to comprehend. 
2. The innovative approach of learning a diffused log-density via the log-density-PF equation is novel.
3. The math derivation is rigorous.
summary: This paper introduces a novel diffusion-based sampling algorithm called the Diffusion-PINN Sampler (DPS). The DPS estimates the drift term in the reverse diffusion process by solving the governing partial differential equation of the log-density of the underlying stochastic differential equation (SDE) marginals using physics-informed neural networks (PINN). The authors prove that the error of log-density approximation can be controlled by the PINN residual loss, which allows them to establish convergence guarantees for the DPS. Experiments on sampling tasks demonstrate the effectiveness of this method.
weaknesses: 1. When the score of the target distribution is accessible, we can perform sampling using MCMC[1] and ParVI[2][3][4]. Renowned algorithms such as HMC and SVGD should be considered for comparison. It would be interesting to see if PINN-Sampling can surpass these methods.

2. The training objective in equation 14 requires gradients of the neural network. This could potentially be time-intensive as the dimension increases. I am skeptical about the scalability of this method. Could you provide a complexity analysis of the training process in relation to dimensionality?

3. Assumption 1 is not plausible, as the suggested scenario is unlikely to occur. Consequently, it's not feasible to bound the error within a constrained domain. It's probable that the subsequent theorems will also hold when $\Omega = \mathbb{R}^d$ and $\nu_t$ are light-tailed distributions.

4. In Example 1, as $\tau$ approaches zero, the variance of $\pi^M$ and $\hat{\pi}^M$ becomes unbounded. However, in real-world applications, the variance of our data distribution is always bounded. Could you give an example where the data distribution has a bounded variance but still has an arbitrarily small $\tau$?

5. The scope of the experiments appears to be rather limited, being mostly synthetic and low-dimensional. It would be beneficial to see applications to higher-dimensional cases, Gaussian processes, Bayesian Neural Networks, and real-world data experiments.

References:

[1] Neal R M. MCMC using Hamiltonian dynamics[J]. arXiv preprint arXiv:1206.1901, 2012. 

[2] Liu Q, Wang D. Stein variational gradient descent: A general-purpose Bayesian inference algorithm[J]. Advances in neural information processing systems, 2016, 29. 

[3] Liu C, Zhuo J, Cheng P, et al. Understanding and accelerating particle-based variational inference[C]//International Conference on Machine Learning. PMLR, 2019: 4082-4092. 

[4] Wang F, Zhu H, Zhang C, et al. GAD-PVI: A General Accelerated Dynamic-Weight Particle-Based Variational Inference Framework[C]//Proceedings of the AAAI Conference on Artificial Intelligence. 2024, 38(14): 15466-15473.

Meta-review: This paper uses physics-informed neural networks to establish a diffusion-based sampling algorithm. There are two salient concerns from the reviewers about this submission. One is that physics informed neural networks (PINN) have already been used in diffusion-based samplers. It is not yet clear that this paper's modification / improvement is significant enough (measured in terms of effectiveness against existing benchmarks) to grant it a publication at ICLR. Another concern is about the scalability of PINN. I would suggest that the authors take another pass to address the above mentioned problems.

Decision: Reject