Thank you for your careful review and helpful feedback! We address your concerns as follows.

1. When the score of the target distribution is accessible, we can perform sampling using MCMC[1] and ParVI. ...

Additional experimental results for more baselines, including HMC and SVGD, are provided in the global response.

2. The training objective in equation 14 requires gradients of the neural network. This could potentially be time-intensive as the dimension increases. ...

We propose an unbiased Hutchinson gradient estimator to reduce the computational cost of evaluating the Laplacian in the training objective. This approach effectively addresses the issue of time-intensive computations that arise as the dimensionality increases, which would otherwise hinder the scalability of our method, as demonstrated in the numerical complexity analysis presented in the global response.

3. Assumption 1 is not plausible, as the suggested scenario is unlikely to occur. ...

First, making or using this assumption is common in the analysis of PINN, as seen in [1] and [2]. Second, this assumption is essential for two reasons: (1) it ensures the uniqueness of the solution within the bounded region, and (2) it guarantees the existence of the continuous constants required for our analysis. Moreover, it is still worthwhile to explore the possibility of extending our analysis to the domain $\Omega = \mathbb{R}^d$ and certain light-tailed distributions $\nu_t$. Specifically, we could construct a continuous modifier function $\delta_\theta^*(\boldsymbol{x}, t)$ such that $\delta_\theta^*(\boldsymbol{x}, t) = u_t^*(\boldsymbol{x}) - u_\theta(\boldsymbol{x}, t)$ for any $(\boldsymbol{x}, t) \in \partial \Omega \times [0, T]$. This would allow us to extend our analysis by leveraging the light-tailed property of $\nu_t$ and the construction of $\delta_\theta^*(\boldsymbol{x}, t)$.

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. ...

We provide an example that exhibits the same properties as in Example 1, with bounded variance. Specifically, we consider the mixture distributions:

then we obtain the following lower bound for the KL divergence:

Additionally, we provide an upper bound for the Fisher divergence:

Finally, for given $\boldsymbol{a}_1, \boldsymbol{a}_2 \in \mathbb{R}^d$ and as $\sigma^2 \to 0$, we observe that $\tau$ approaches zero in Example 1, while the variance of both $\pi^M$ and $\hat{\pi}^M$ remains bounded.

5. The scope of the experiments appears to be rather limited, being mostly synthetic and low-dimensional. ...

We present additional high-dimensional experiments on a 50-dimensional Double-Well with 32 separated modes, incorporating challenging weights in the global response.

References

[1] Deveney T, Stanczuk J, Kreusser L M, et al. Closing the ODE-SDE gap in score-based diffusion models through the Fokker-Planck equation[J]. arXiv preprint arXiv:2311.15996, 2023.

[2] Wang C, Li S, He D, et al. Is $L^ 2$ Physics Informed Loss Always Suitable for Training Physics Informed Neural Network?[J]. Advances in Neural Information Processing Systems, 2022, 35: 8278-8290.

Thanks for your clarification!

Do you plan to give a revised manuscript?

Thanks for your reply. I also have a further concern about the newly provided experiments.

1. Your algorithm can achieve a 90% distance reduction on your provided experiments compared to HMC or SVGD. How many particles you use for SVGD, this is a key hyperparameter.

2. Why does HMC perform worse than LMC? From my knowledge, HMC is a momentum-augmented accelerated version of LMC; their convergence behavior after enough iterations should be similar.

3. These MCMC or ParVI algorithms need mix time, I want to see the convergence curve and detailed hyperparameters setting of these baselines in the revised paper, or I cannot be convinced. Or can you provide your code with scripts?

4. I also want to know the visual time-costs comparison, how much time do you need to train on the 9-gaussian level task?

Thanks for your response! We have addressed your concerns as follows.

1. How many particles used for SVGD ...

The 9-Gaussians, Rings, and Double-Well targets in our experiments are all multi-modal with challenging (unequal) mixing proportions, making them particularly difficult for traditional particle-based variational inference methods like SVGD to sample from. Additionally, the time and memory costs for SVGD increase significantly with the number of particles. Therefore, we chose to use 1,000 particles for SVGD in our experiments, a standard and practical setting for this method. This hyper-parameter choice has been explicitly included in our revised paper.

2. HMC ...

First, we want to clarify that we utilize a parallel version of LMC in our paper. Specifically, we run multiple independent LMC chains, one for each initial sample. For example, to generate 10,000 samples from the target distribution, we begin with 10,000 initial points (e.g., sampled from a standard Gaussian) and run Langevin dynamics independently for these 10,000 points over several iterations (e.g., 100,000). The final iteration outputs of the LMC chains are then collected as samples from the target distribution. This parallel approach leverages the diversity of initializations and independent chain updates to produce significantly higher-quality samples than standard LMC.

Previously, we used a non-parallel (standard) version of HMC as an additional baseline, which resulted in poorer performance compared to parallel LMC. We have updated the HMC implementation to a parallel version and fine-tuned the hyper-parameters (step size = 0.1, number of Leapfrog iterations = 10, total iterations = 10,000), achieving results that are comparable to or better than those of LMC. Additionally, we have clarified our use of the parallel version of LMC in the revision to avoid any misunderstanding.

3. convergence curve ...

The convergence curves for well-trained DPS, SVGD, and HMC for 9-Gaussians are presented in the table below. Notably, leveraging the reverse process of diffusion models, our sampling algorithm achieves significantly faster convergence compared to MCMC and ParVI methods. In contrast, the MCMC method (HMC) converges to the target distribution at an extremely slow rate.

Table: Convergence curves of well-trained DPS, SVGD, and HMC for 9-Gaussians (KL divergence vs. number of steps).

4. time-costs ...

The following table compares the training time (in iterations) and sampling accuracy (measured by KL divergence) achieved by our DPS. On the 9-Gaussians experiment, DPS requires only 20,000 training iterations—approximately 12 minutes—to achieve nearly perfect sampling, demonstrating its high efficiency.

Table: Training time (iterations) vs. sampling accuracy (KL divergence) for DPS on 9-Gaussians.

Thank you for your thoughtful review and helpful feedback! We address your concerns as follows.

1. Novelty ...

The application of PINN to diffusion-based sampling has gained significant attention in recent research. However, it is important to highlight the unique contributions of our work compared to related studies. First, we utilize LMC for generating collocation points, which proves to be more efficient, particularly in high-dimensional settings, compared to uniform sampling. This claim is supported by our theoretical results in Section 5.3. Additionally, our work identifies the limitations of using PINN to solve the score Fokker–Planck equation (FPE) for mixing proportion identification in sampling tasks. Furthermore, we provide a comprehensive convergence analysis for our proposed algorithm. These points collectively demonstrate the novelty of our approach. A more detailed discussion of related work will be included in our revision.

2. Numerical Evaluation ...

First, we extend our experiments on the Double-Well target to 50 dimensions, consistent with the setup in [3]. The results are shown in the global response. Notably, our target distribution is characterized by challenging mixing proportions, making it particularly difficult to sample from. While we recognize that solving higher-dimensional PDEs using PINN can be challenging, recent works have introduced techniques to scale PINN to higher-dimensional scenarios. Incorporating some of these methods could help address this limitation in our approach.

3. Practical relevance of theoretical results ...

It is intuitively clear that achieving zero PINN loss would lead to zero score error and perfect sampling via diffusion models. However, in practice, attaining zero PINN loss is impossible. Our theoretical results aim to provide a quantitative error analysis for PINN and a convergence guarantee for our sampling algorithm. Specifically, we bound the score error in terms of the PINN residual loss and upper-bound the sampling error in KL divergence by the PINN loss.

While the constants in our analysis cannot be explicitly computed, the convergence results remain valid. Moreover, some constants can be controlled through techniques such as weight clipping.

We agree that the analysis presented in this work may not be very tight. Improving the analysis and tightening the bounds is an avenue for future exploration.

Additionally, while we derive our results using techniques from [6] and theoretical results on diffusion models, it is important to highlight that we extend the analysis in [6]. Specifically, we consider the distribution of collocation points $\nu_t$, whereas [6] only addresses uniform sampling which would become inefficient in higher dimensions.

4. It is unclear to me how generating collocation points with MCMC-like algorithms (e.g. LMC) is sufficient. ...

In general, obtaining perfect samples using MCMC methods is challenging. However, for our purposes, we only require samples that adequately cover the high-density domain, which is more feasible. In this context, MCMC methods are sufficient. Specifically, in practice, a short LMC chain with a large learning rate can generate samples that effectively cover the high-density domain, even though it may struggle to accurately capture the correct mixing proportions for multimodal distributions. Nonetheless, these `rough' samples are sufficient for training within the PINN framework for accurate log-potential estimation, as demonstrated by our theoretical analysis and numerical experiments (our method can capture the correct mixing proportions).

5. You write that trajectory-based alternatives lead to "computational complexity associated with differentiating through SDE solvers" ...

We acknowledge that certain modifications can address this drawback. However, it is worth noting that these techniques have only emerged recently. We will revise our discussion accordingly in the next version.

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.

The related works section primarily summarizes and compares recent diffusion-based sampling techniques via solving PDEs. Other diffusion-based methods, such as RDMC, PIS, DIS, and DDS, are discussed in the introduction. We will revise our discussion of related works on diffusion-based sampling in the next version and sincerely appreciate the valuable comments and suggestions regarding additional references.

11. 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. ...

Our statement demonstrates that the KL divergence can be (much) larger than the Fisher divergence. To build intuition, consider the mixture model

where $\pi_1(x)$ and $\pi_2(x)$ are two distributions with well-separated high-density regions, and $w_1, w_2 > 0$ such that $w_1 + w_2 = 1$. The score function of this mixture model is given by

Let $\Omega_1$ and $\Omega_2$ denote the high-density regions of $\pi_1(x)$ and $\pi_2(x)$, respectively, which are assumed to be well-separated. Suppose $\pi_1(x)$ and $\pi_2(x)$ are light-tailed distributions (e.g., Gaussian) and satisfy

The high-density region of the mixture model $\pi(x)$ is thus $\Omega_1 \cup \Omega_2$. Within this domain, for $x \in \Omega_1$, we approximately have

and for $x \in \Omega_2$,

Consequently, the score function $\nabla_x \log \pi(x)$ provides little information about the weights $w_1$ and $w_2$ on the high-density regions of $\pi(x)$, leading to a nearly zero Fisher divergence. In contrast, the KL divergence between two mixture models with different weights $w_1$ and $w_2$ remains a positive constant. This observation highlights a fundamental difference between the two divergence measures and suggests that it is general for mixture models with light-tailed distributions, including Gaussians. We want to emphasize here that it is not the scale of different divergence that matters here. It is all about the power of distinguishing the weights of isolated components for different divergence, which in turn affect the magnitude of the corresponding signals when used for optimization. This is known as blindness of score matching in the literature ([i], [ii]).

In Figure 1, the Fisher divergence is observed to be very small, approaching zero in our numerical implementation. Hence, changing the scale should not affect this result.

Additionally, recall that if a distribution $\mu$ satisfies the logarithmic Sobolev inequality (LSI) with constant $C_{\rm LSI}$, we have

for any distribution $\nu$, where the Fisher divergence is defined as

Therefore, our observation aligns with the fact that the LSI constant $C_{\rm LSI}$ for the mixture model can be large (even unbounded) when the high-density regions $\Omega_1$ and $\Omega_2$ are well-separated.

[i] L. K. Wenliang. Blindness of score-based methods to isolated components and mixing proportions. arXiv preprint arXiv:2008.10087, 2020.

[ii] Mingtian Zhang, Oscar Key, Peter Hayes, David Barber, Brooks Paige, and Francois-Xavier Briol. Towards healing the blindness of score matching. In NeurIPS 2022 Workshop on Score-Based Methods, 2022.

Thanks for your response! We have addressed your concerns as follows.

1. Novelty ...

Compared to [1], where the combination of PINN and sampling is only briefly mentioned in the appendix, our work provides a comprehensive investigation of diffusion-based sampling using PINN. This includes investigating limitations highlighting the drawbacks of using the score Fokker-Planck equation (FPE) in diffusion-based sampling to identify mixing proportions, which is known as blindness of score matching, as well as offering a complete theoretical analysis of the sampling algorithm combined with the PINN approach.

In contrast to [2], it is important to note that [2] considers a linear interpolation path (in log space) for sampling, which differs from the reverse process of diffusion models employed in our work.

Finally, in comparison to [3], which is a concurrent work to ours (notably, our work was submitted to NeurIPS 2024 and was very close to an acceptance with an average score of $6$!), our study includes additional theoretical analysis and focuses on more challenging situations on identifying mixing proportions for target distributions with isolated modes, setting it apart in terms of depth and contribution.

2. Numerical evaluation ...

The scalability of our approach primarily depends on the effectiveness of PINN in high-dimensional settings. Fortunately, recent advancements ([1,2]) have provided several methods that successfully improve the ability of PINN to solve relatively high-dimensional PDEs. These scaling techniques for PINN can be directly incorporated into our sampling framework, enabling our sampling algorithm to scale effectively as well.

We did not include the results of the log-variance versions of PIS and DIS because their execution time, using the officially released implementation, exceeds two days. Nonetheless, we agree that PIS-LV and DIS-LV would perform comparably to our method.

3. Practical relevance ...

We want to emphasize that our theoretical results establish a connection between the approximation error of both the log-density and the score in terms of the PINN residual loss. Specifically, if the PINN residual loss is sufficiently small, then the approximation errors for both the log-density and the score will also be small.

Intuitively, the PINN residual loss provides a higher-order constraint on the log-density approximation, which is different from a simple $L_2$ loss. In order to get a small PINN residual loss, not only the log-density approximation needs to be accurate, but also its higher-order gradients (including the score) that are involved in the PINN residual evaluation need to be accurate as well. To illustrate this, we provide the following example. In our implementation, the PINN residual is defined as:

For simplicity, let’s consider the one-dimensional case where $u_t(x) = -0.5x^2$ . In this case, the residual becomes:

However, if we introduce a perturbation $\tilde{u}_t(x) = -0.5x^2 + 0.1\sin(100x)$ , the PINN residual becomes:

which can be significantly large. This example is consistent with our theoretical findings, demonstrating that a small PINN residual loss (not the $L_2$ loss) is critical for ensuring accurate (score) approximations.

Furthermore, we conducted numerical experiments to address your concern. Specifically, we computed the PINN residual loss for both the well-trained model and perturbed well-trained model on the 9-Gaussian target. The results are summarized in the following table:

Table: Comparison of PINN residual between well-trained models with and without slight perturbations

From the above results, it is evident that the perturbation causes a significantly large PINN residual loss, despite only a slight variation in the log-density. Due to the large PINN residual loss, achieving a small score approximation error is not feasible.

4. MCMC algorithms. ...

Traditional MCMC algorithms typically require an extremely long time to generate perfect samples, especially when dealing with high-dimensional multi-modal target distributions. This is due to the inherent trade-off between exploration and accuracy, which leads to prolonged mixing times. However, in our work, we use MCMC methods (e.g., LMC) solely to generate collocation points, with the goal of sufficiently covering the entire high-density domain rather than obtaining perfect samples. This relaxed requirement allows us to adopt a larger step size for LMC, enabling sufficient exploration while eliminating the original trade-off. This approach highlights a key advantage of off-policy training in our framework.

5. Comparable results with other papers ...

It is important to highlight that, in our experiments, the multi-modal target distributions are designed with challenging (unequal) mixing proportions, making them more difficult to sample from. This may explain why most baseline methods tend to perform worse.

6. Fisher v.s. KL ...

Yes. Our results are mainly for mixture models. In this paper, our goal is to develop an efficient sampling algorithm for multi-modal distributions, which is challenging for traditional methods. For simpler distributions (i.e., single-mode distributions), traditional MCMC methods such as HMC would be enough to provide good samples.

References:

[1]. Zheyuan Hu, Zekun Shi, George Em Karniadakis, and Kenji Kawaguchi. Hutchinson trace estimation for high-dimensional and high-order physics-informed neural networks. Computer Methods in Applied Mechanics and Engineering, 424:116883, 2024a.

[2]. Zheyuan Hu, Khemraj Shukla, George Em Karniadakis, and Kenji Kawaguchi. Tackling the curse of dimensionality with physics-informed neural networks. Neural Networks, pp. 106369, 2024b.

Thanks for your response! We address your concerns as follows.

1. Novelty ...

We agree that [2] also considered learnable paths other than linear interpolation in the log space. However, learning these paths requires extra computation and may incur instability during training.

3. Practical relevance ...

Apologies for the misunderstanding. To clarify, our original goal is to show that the approximation error for the score function can be effectively controlled as long as the log-density is well-approximated by minimizing the PINN residual.

4. MCMC ...

As mentioned before, in our case, we are not running traditional MCMC to get samples from the target distribution for the collocation points. Due to the off-policy nature of our method, we can run a much shorter and noisier (and parallel) MCMC (e.g., LMC) with large step sizes that are more likely to capture the high-probability regions. This is more practical than uniform sampling from a pre-selected high probability region (adopted in previous PINN-based sampling algorithms) that are often hard to obtain, especially in high dimensions.

Thanks for your suggestion! We will modify the name accordingly in our revision.

1. I'm a little bit concerned about the performances of the method both in time and resources. ...

We provide measurements of training time and sampling time across different dimensional cases in the complexity analysis section of the global response.

2. You model the log-density, then you have to take its gradient in order to compute the score ...

In the complexity analysis section of the global response, we compare the sampling time using direct score estimation versus taking the gradient of the approximated log-density. Our results show that sampling time is halved when score estimation is directly employed. Nonetheless, the sampling time of our method is already highly efficient. For instance, sampling 10,000 points over 1,000 time steps takes less than 2 seconds.

3. The PINN loss involves gradient and Laplacian operators. ...

Our training objective is inherently more complex than the simple denoising score matching objective, as it includes terms like Laplacian operators. Consequently, our training process is slower compared to simple denoising score matching. Notably, we have not yet incorporated optimized techniques such as jvp/vjp for implementation. Instead, we address the computational challenges by proposing an unbiased Hutchinson gradient estimator to avoid the costly evaluation of the Laplacian. This estimator significantly reduces computational costs and accelerates training, as detailed in the numerical complexity analysis section of the global response.

Moreover, it is important to highlight a key distinction: simple denoising score matching relies on access to a large amount of training data sampled from the target distribution, which is available in the context of generative modeling. In contrast, our focus is on sampling from an unnormalized density, where such training data is unavailable. As a result, we have to utilize a more complex objective to accurately approximate the score function in this setting.

4. In table 1, isn't it misleading to report KL only on first dimensions for Funnel and Double-well? ...

This evaluation standard leverages the properties of the target distributions (Funnel and Double-Well). In the Funnel case, given the first dimension, the remaining dimensions are independent and identically distributed. Similarly, in the Double-Well case, the distribution is dimension-independent, with all dimensions following a standard Gaussian distribution except for the first three (in the 30-dimensional case) or the first five (in the 50-dimensional case).

As a result, in both cases, analyzing the KL divergence over the first five dimensions is sufficient to capture the differences between the target and approximate distributions. This approach effectively highlights the key discrepancies while avoiding unnecessary computation across redundant dimensions.

5. 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. ...

In this work, we focus on sampling from an unnormalized density, such as posterior sampling in Bayesian inference. In this setting, we assume access to the exact formula of the unnormalized density, but we do not have access to samples from the target distribution. Therefore, our approach inherently requires access to the unnormalized density, unlike generative modeling, which typically relies on samples from the target distribution.

6. 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 ...

The main purpose of our paper is to leverage the reverse process of diffusion models to design efficient sampling algorithms, instead of generative models. Therefore, we do not expect our method to be studied in real high dimensional data such as images. In those cases, one could use diffusion models directly and with data samples, score approximation can be easily obtained with denoising score matching.

Furthermore, leveraging PyTorch’s vectorized implementation ensures that the cost of score computation during sampling remains efficient. While it is slower than using direct score estimation, it is still highly performant—sampling 10,000 points over 1,000 time steps for 100-dimensional tasks takes less than 2 seconds. In addition, we introduce an unbiased Hutchinson gradient estimator to mitigate the high computational cost associated with evaluating the Laplacian term in the PINN loss. Our numerical results, presented in the complexity analysis section of the global response, substantiate these claims.

For higher-dimensional problems, see the results in our global response.

Thank you for your thoughtful review and valuable feedback. We address your specific questions and comments below.

1. 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 ...

First, we propose an unbiased Hutchinson gradient estimator to compute the PINN loss, significantly reducing the computational cost associated with evaluating the Laplacian in the loss function. Additionally, we compare the sampling time using direct score estimation versus taking the gradient of the approximated log-density. Our results show that sampling time is halved when score estimation is directly employed. Nonetheless, the sampling time of our method is already highly efficient. For instance, sampling 10,000 points over 1,000 time steps takes less than 2 seconds. For further discussion and additional experimental results, please refer to the complexity analysis section of the global response.

2. 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.

We have included additional experimental results for modern baselines, such as HMC and the recent diffusion-based methods introduced in [2], in the global response.

3. 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.

We present additional higher-dimensional experiments on $50$-dimensional Double-Well with $32$ separated modes, following the same setup as in [1] but with more challenging weights in the global response.

4. Parametrize NN with log mu (12) seems like a strong inductive bias ...

First, our parameterization ensures that the parameterized solution satisfies the initial condition, eliminating the need to consider the initial loss in our training objective. However, we agree that this parameterization may incur higher computational costs, especially when evaluating the energy of the target distribution is expensive.

In such cases, we can use a simple neural network and incorporate the initial loss term into the training objective, which can be computed via unbiased stochastic estimation to reduce the cost of evaluating the energy. Namely, we use the following training objective,

Notably, in our experiments, the density function can be computed efficiently, so we chose this parameterization to avoid the additional burden of the initial loss term.

Additionally, we conducted an ablation study on the 9-Gaussians task. The results demonstrate that including the initial loss in the training objective, rather than using this parameterization, still provides a good score approximation (and naturally achieves effective sampling). These experimental results are presented in the ablation study section of the global response.

5. Related works Section is rather short and should be extended.

Thanks for the suggestion! We want to mention that the introduction section also discussed several related methods, including VI methods, MCMC methods, RDMC, PIS, DIS, and DDS. In the related works section, we further discuss other diffusion-based sampling techniques that involve solving ODEs/PDEs, to highlight the unique contribution of our approach. We will add more discussion to related works in our revision. Given that this is a rapidly growing research area, we admit that there may be some missing references and would greatly appreciate it if you could provide more related works as well.

Thanks for raising the score! For further discussions on the scalability and novelty of our work, please refer to the global response!