Functional Gradient Flows for Constrained Sampling
Authors: Shiyue Zhang
Longlin Yu
Ziheng Cheng
Cheng Zhang
Venue: NeurIPS 2024 poster (2024)
Keywords: particle-based variational inference
constrained sampling
functional gradient flow
boundary integral
Primary area: probabilistic_methods
TL;DR: A new functional gradient particle-based variational inference method for sampling on constrained domains.

Abstract: Recently, through a unified gradient flow perspective of Markov chain Monte Carlo (MCMC) and variational inference (VI), particle-based variational inference methods (ParVIs) have been proposed that tend to combine the best of both worlds. While typical ParVIs such as Stein Variational Gradient Descent (SVGD) approximate the gradient flow within a reproducing kernel Hilbert space (RKHS), many attempts have been made recently to replace RKHS with more expressive function spaces, such as neural networks. While successful, these methods are mainly designed for sampling from unconstrained domains. In this paper, we offer a general solution to constrained sampling by introducing a boundary condition for the gradient flow which would confine the particles within the specific domain. This allows us to propose a new functional gradient ParVI method for constrained sampling, called *constrained functional gradient flow* (CFG), with provable continuous-time convergence in total variation (TV). We also present novel numerical strategies to handle the boundary integral term arising from the domain constraints. Our theory and experiments demonstrate the effectiveness of the proposed framework.

Functional Gradient Flows for Constrained Sampling
Shiyue Zhang"
Longlin Yu" School Mathematical Sciences Peking University 11 ylyu@pku.ed eduep.uu.dd..com edu. cn
School of Mathematical Sciences
Peking University zhangshiyuedstu, pktpisstt.pkuuu.... edu.
Ziheng Cheng"
Department of Industrial Engineering and Operations Research University of California, Berkeley
beaheng_chenggberkeeeeeee edu edu
Cheng Zhang!
School of Mathematical Sciences and Center for Statistical Science Peking University
chengzhang@math, pku. edu. cn
Abstract
Recently, through a unified gradient flow perspective of Markov chain Monte Carlo (MCMC) and variational inference (VI), particle-hased variational inference methods (ParVIs) have been proposed that tend to combine the best of both worlds. While typical ParVIs such Stein Variational Gradient Descent (SVGD) approxi- mate the gradient flow within a reproducing kernel Hilbert space (RKHS), many attempts have been made recently to replace RKHS with more expressive function spaces, such as neural networks. While successful, these methods are mainly de- signed for sampling from unconstrained domains. In this paper, we offer a general solution to constrained sampling by introducing a boundary condition for the gradi- ent flow which would confine the particles within the specific domain. This allows us to propose a new functional gradient Par method for constrained sampling, called constrained functional gradient flow (CFG), with provable continuous-tine convergence in total variation (TV). We also present novel numerical strategies to handle the boundary integral term arising from the domain constraints. Our theory and experiments demonstrate the effectiveness of the proposed framework.
Introduction
Efficiently approximating and sampling from unnormalized distributions is a fundamental and challenging task in probabilistic machine learning, especially in Bayesian inference. Various methods, including Markov Chain Monte Carlo (MCMC) and Variational Inference (VI), have been developed to address the intractability of the target distribution. VI, the inference problem reformulated as optimization task that aims to find approximation within specific distribution family that minimizes the Kullback-Leibler (KL) divergence to the posterior (Jordan et al. 1999 Wainwright
& Jordan 2008 Blei al. 2016). Leveraging efficient optimization algorithms,, VI is often fast during training and more scalable to large datasets. However, its approximation power may be limited depending the chosen family of variational distributions. In contrast, MCMC methods generate
"Equal contribution.
t Corresponding author.
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

samples from the posterior through Markov chain that satisfies the detailed balance condition (Duane et al. 1987 Robert & Stramer, 2002, Neal 201T, Welling & Tch 2011 Chen et al. 2014). While MCMC is asymptotically unbiased, it may be slow to converge, and assessing convergence can be challenging.
Recently, there has been growing interest in the gradient flow formulation of both MCMC and VI, leading to the development of particle based variational inference methods (ParVIs) that tend to combine the best of both worlds (Liu & Wang 2016, Chen et al. 2018 Liu et al. 2019 di Langosco
2022). From the variational perspective, ParVIs take a non parametric approach where the approximating distribution is represented as set of particles.
at. 2021 Fan et al. 2022 Alvarez-Melis
These particles are iteratively updated towards the steepest direction to reduce the KL divergence to the posterior, following the gradient flow in the space of distributions with certain geometries. This non-parametric nature ParVIs significantly enhances its flexibility compared to classical parametric VIs, and the interaction between particles also makes VarVIs more particle-efficiet than MCMCs.
One prominent particle-based variational inference technique Stein Variational Gradient Descent (SVGD) (Liu & Wang (016). It calculates the update directions for particles by approximating the gradient flows of the Kulbback-Leibler (KL) divergence within a reproducing kernel Hilbert space (RKHS), where the approximation takes a tractable form (Liu 2017 Chewi et al. (020). However, the performance of SVGD heavily depends on the choice of the kernel function and the quadratic computational complexity of the kernel matrix also makes it impractical to use a large number of particles. As kernel methods are known to have limited expressive power, many attempts have been made recently to expand the function class for gradient flow approximation (Hu et al., (0188) Grathwohl et al. 2020 di Langosco al., 2021 Dong et al., 2223 Cheng et al., (023). By embracing a more expressive set of functions, such as neural networks, these functional gradient approaches have shown improved performance over vanilla SVGD while not requiring expensive kernel computation.
While MCMC and VI methods have shown great success in sampling from unconstrained domains, they often struggle when dealing with target distributions supported on constrained domains. Sam- pling from constrained domains is an important and challenging problem that appears in various fields, such as topic modeling (Blei al., 2003), computational statistics and biology (Morris. 2002) Lewis et al. 2012 Thiele et al. 2013). Recently, several attempts have been made to extend classical sampling methods like Hamiltonian Monte Carlo (HMC) or SVGD to constrained domains (Brubaker et al. 2012 Byrne& & Girolami 2013 Liu & Zhu 2018 Shi et al. (022). However, these extensions either involve computationally expensive numerical subroutines such as solving nonlinear systems of equations, or rely on intricate implicit and symplectic schemes or mirror maps that require case-by-case design effort tailored to specific constraint domains.
In this paper, we propose a functional gradient ParVI method for sampling from probability distribu- tions subject to constrained domains with general shapes. We demonstrate that functional gradient approaches for Var can be seamlessly adapted to constrained domains by learning the gradient flows with a vector field that adheres boundary condition. Intuitively, this boundary condition would confine particles within the specified domain, and it is indeed a sufficient condition for gradient flows of probability measures confined to this domain. Following previous works (di Langosco et al., 2021; Dong et al. 2023) Cheng et al. 2023). we employ the regularized Stein discrepancy objective for functional gradient estimation where we directly incorporate the boundary condition for the vector field into the design of its neural network approximation. Due to the domain constraints, integration by parts now would lead to an additional boundary integral term in the training objective. We, there- fore, derived an effective approach for properly evaluating this term for general boundaries. Extensive numerical experiments across different constrained machine learning problems are conducted to demonstrate the effectiveness and efficiency of our method.
2 Related Work
To better capture the gradient flow, number of functional gradient methods for ParVIs have been proposed recently that employ a larger and more expressive class of functions than reproducing kernel Hilbert spaces (RKHS). In particular, di Langosco et al. (202T) used neural networks to learn the Stein Discrepancy with an L2 regularization term, and updated the particles based on the learned witness functions. Dong et al. (2023) and [heng et al. (2023) provided extensions that accommodates broader class of regularizers.

---

Sampling from constrained domains is generally more challenging compared to the unconstrained ones. [Brubaker et al. (2012) proposed constrained version of HMC for sampling on implicit manifolds by applying Lagrangian mechanics to Hamiltonian dynamics, which requires solving
a nonlinear system of equations for each numerical integration step. Lan al. (2014) focused on constraints that can be transformed into hyper-spheres, which can be viewed special case of Brubaker et al. (2012) that has close-form update formulas. Zhang al. (2020); Ahn Chewi (202I);
Shi et al. (2022) extended Langevin algorithms and SVGD to constrained domains via mirror maps. However, these methods require explicit forms of transformations (e.g., spherical augmentation and mirror maps) that capture the constraints, which would limit their applications simple domains. Bubeck et al. 2018); Brosse et al. (017);; y; Salim & Richtárik (2020) considered projected Langevin algorithms through projection oracle, which of high computational cost for complex constrained domains. Zhang et al. (2022) proposed an orthogonal-space gradient flow approach for sampling in manifold domains with equality constraints, which employed similar strategy ours. Our method different from theirs in that it is designed for domains with inequality constraints using ParVIs with neural network gradient flow approximations, and we also proposed novel numerical strategies to address boundary-related challenges.
3 Background
Notations We use denote particles in Rd and (xlg(x) S 0} denote the constrained domain. Notation 1n is the indicator function  Let || denote the standard Euclidean norm of vector. Let ((Rd) denote the set of probability distributions on Rd that are absolute continuous with respect the Lebesgue measure, and we do not distinguish the probability measure with its density function. Let DKL. be the Kullback-Leiibe divergence and TV be the total variation distance, We use C(X, y) denote the space of continuous mappings from N' to V, and the shorthand C(X) for C).
3.1 Particle-based Variational Inference
Let p" E P(RR)) be the target distribution we wish to sample from. We can frame the problem of sampling into a KL divergence minimization problem
arg min (kl(p) pEo
(I)
where Q C P(RR)) is the space of probability measures. Particle-based variational inference methods (la)VIs) can described within this framework where Q represented as a set of particles. Starting from an initial distribution po and an initial particle zo Po. we update the particle xz following dæ, = Ue (2, (dt where 0 : IR + Rd is the velocity field at time t. The density pt of x, follows the continuity equation dp:/ /dt = - (trPt), and the KL divergence decreases with the following rate:
d (PL(Pelle") = - Ep, ( log vt). "",)))
P"
di
VarVIs aim to find the optimal velocity field that minimizes (2) in Hilbert space H with the squared norm regularizer as follows
- Ep. (V log 2, wx) + 11 vz) pt.
min NEW
(3)
3.2 Functional Wasserstein Gradient
Different choices of the space H in (3) lead to different gradient flow algorithms. SVGD Liu & Wang (2016) chooses to be Reproducing Kernel Hilbert Space with kernel k(.,)). This way, the optimal velocity field has close form solution
2 (-) Ep, V logr" (x) + v,k(;))],
albeit the design of kernels can be restrictive for the flexibility of the method.
When 7/ = _ (p:), then v.- = V log V. is also known as the Wasserstein gradient of KL divergence (Jordan al. 1998). Since the score function of particle distribution V log Pt is generally inaccessible.

---

recent works propose to parameterize as a neural network and train it through (3 di Langosco et al. (202I) Dong et al. (2023); Cheng et (2023). Due to Stein's 's identity, ( has tractable form
2/ = arg min Bp. #EF
(V logp", "eep","" .)) V
||ll
(5)
2
where F is the neural network family. This method is called functional Wasserstein gradient in contrast to the kernelized version.
4 Main Method
Consider sampling from a target distribution (*(x) supported on Ω {g(x) S 0}, where -
R is continuously differentiable function. Throughout the work, we assume that p* admits differentiable positive density function on Ω
We first reformulate the problem into constrained optimization in the space of probability measures:
min PEP(R4)
(KI,()lp"),
((2) =
However, note that p" is only supported on Ω2 the problem is ill-posed if p(?) < 1, since in this case p will be absolute continuous with respect to p" and (pll""))) cannot be defined. Besides, note that the initial particle distribution po is generally assumed be fully supported on Idd and thus does not meet the constraint.
To fix this issue, natural idea is to drive the particle distribution pı towards Ω to satisfy Pc(2) = 1 and train the functional gradient flow on constrained domains by minimizing the regularized Stein discrepancy (RSD)
p" 1 Pt (- log 2 2 + [111113) dx.
min LRSD TEF
These stringent requirements make the construction of velocity field for constrained sampling far more volatile than conventional functional gradient for unconstrained domains.
In the rest of this section, we first present the idea of velocity design in Section 4. and then derive
the tractable training objective in Section 4.2 The practical algorithm is proposed in Section 4.3
4.1 Necessity of Piece-wise Velocity Field
We first show that globally continuous velocity field may fail achieve the exact minimum of RSD under the constraint. To ensure that particles inside will never escape, the velocity should satisfy the boundary condition
04 0, on an
(8)
In fact, if (8) holds, by Stoke's formula and the continuity equation @ (x) -V (p)(z))(())). we have
d di (())
d Ep, In [ σt (())dr -- f v (po()))((z))de = p((z)v((z) AdS.
dt
Jnn
Thus conclude
Proposition 4.1. vf 7 < 0 on Oh, then ((( will not decrease.
However, it is possible that there does not exist an optimal solution of 7 that meets the boundary condition (8) and continuous on 0Ω, as illustrated the following ID examplee..
Example 4.2. Let p" x ) (2221) and exp( exp(_!"=tt) (m=4)²) We have V log p- &/ 2 in
Ω and LASD = 1 plu(x) (1)+ 2 s. The constraint on boundary is e(1) $ 0, (11) 2 0. Hence there is optimal v in C(R).
1

---

Therefore, it is necessary to design velocities for particles inside and outside separately. For particles located outside S (contributing the probability p(g(z) > 0)), it prudent to make the velocity drive the particles into the domain. On the other hand, for those already inside Q the velocity field should be able to fit the target distribution. Overall, reasonable to use the following piece-wise construction:
Vg |V
hnet (0<00 A
(G220--
(10)
Here > a constant. The gradient direction Vg will push the particles into Q. And htne is continuous neural net to learn the optimal direction in Ω. We observed that the magnitudes of the gradients of different particles vary significantly in the numerical experiments and thus use the normalized gradient Vg instead. Vall
4.2 Training Objective
Based on the velocity design in 10 hnet to leamn the optimal direction in Z through 77.
now ready to derive a tractable training objective for
Note that in principle, lnet should be trained after (22) 1. However, in order to enhance efficiency in practice, we use inner particles minimize RSD hefore the particles are driven the constrained domain. At this phase, we replace pt in 7 with the conditional measure ft = (=(-
For any u E ((S2,,R.) expand and we get
pv T v log P' dx + 2.1o Pellu||drr
LASD
Pe.
2
(11)
1 - logp" +V vjdr 2 pr 11/1//ddx
/ dev ödS. Jan
Jn
Substitute with the continuous extension of hnet Q2 and we obtain the training objective for Anet-
Note that the first two terms in [1] are aligned with (5) by substituting in the integration domain, while the last term is a boundary integral specific to inequality constraints. This is one of the essential differences between training functional gradients on constrained and unconstrained domains.
4.3 Estimation of the Boundary Integral
The boundary integral Lo po- -id is computationally intractable in general. Inspired by the co-area formula (Federer 2014 we can estimate it using heuristic band-wise approximation as follows
/ pc 7dS 2c h S.
pu Ādæ = $ p(z)o(x)
p((x)) dx p(a E Sh) n p(x)
/ p(z)v(x) ndr,
80
h
S,
(12)
where p(x) p(x) such that Jss ((x)dx 1, and Sn (x Ω d(x, 02) < h} is band-like area the boundary with small bandwidth 0, which further has proxy
s(x) < 0, o (
Vg(x) [Vo())] 2 0.
(13)
Please refer to Appendix Afro more explanations. This way, area to obtain a Monte Carlo estimate of the boundary integral. Suppose there are m particles Q of which (x, 1y are the band-like area. Then p(# E Sx)) m.
use the particles in the bake-like
= a and thus
f pu Pu- 80
AdS idS R o(z,) V
(14)
mh
We summarize all the techniques above propose the Constrained Functional Gradient (CFG) in Algorithm I We use the neural network structure designed in Appendix D for Anet-
5

---

Algorithm CFG: Constrained Functional Gradient
Require: Unnormalized target distribution y", initial particles { initial neural network param-
eters We {Nu. )u weight parameter X, iteration number L, L', particle step size a, parameter step size 7). bandwidth h
for k 0.... L 1 do
Assign wo = wkk
Identify {- 1==1
f Shh, (+h [xffat e 12}
for = 0, .0 TT' do
Set Andg = Jnt Et Vg
2
Compute
ERSD(w)
logp' (x
Nnw (x/() (afe))))
m
h (z)) Vg
2
meh
VV((22)
Update w/ 1+1 = w - 1/ Cn end for Assign Wh+11 = w
Compute Uwi+1 (frx-1 Eati Vg) 1 [gcoo - A -V (g20))
Update particles x1/11 2 + 000 (xt) for = =1,.. N end for return Particles (ri)
5 Theoretical Analysis
In this section, we present the convergence guarantee of CFG in terms of TV distance. Consider the continuous dynamic dæ, vy (2, dt where v is defined in 110). Let Pt be the law of 2,- We first investigate when particles can get into the domain.
Assumption 5.1. [V(())] 2 C > 0 for x E Q and Mo max(g(20)) : xo E supp(po)) < co.
For simplicity, we do not delve into the case that the particles outside may get stuck at local extremal point of y before getting in Ω. Based on this assumption, the particles will enter the constrained domain in finite time:
Theorem 5.2. There exists a finite to
Ma. such that g(xx) S O, VI 2 lo- XC
The proof is deferred Appendix B. Since domains with inequality constraints do necessarily require the velocity field to gradually decay near the boundary as in Zhang et al. (2022), directly using A TV (g20) allows for swift penetration of the particles into the constrained domain.
Vy
Once the particles are in the constrained domain, an important technique is to extend the target distribution be non-trivially supported on the whole space by convolution with Gaussian kernel and use it as a proxy. Denote the convolution the target distribution and the Gaussian distribution as 0- p= * (O, o²1), where o can be arbitrarily small. Now the KL divergence DKL(pz | is well defined on IRd. Moreover, by the convolution properties, TV Ipp)) + 0 as 0 Stein & Shakarchi 2005).
Assumption 5.3 (Poincaré inequality). The target distribution p" satisfies var,- [f] $ EEE | for any smooth function f.
This is common assumption in convergence analysis that describes the degree of convexity of the target distribution (Chewi et al. 2022). contrast to Zhang et al. (2022). our approach relies solely on the Poincaré property of the original target distribution, rather than requiring it for the conditional measure on zero measure set. Following Cheng al, (2023), we make the following assumption on the approximation error of neural networks.

---

Assumption 5,4. For any to. Pell/nec(t) V log en 1|
Now, we are ready to present the main result.
Proposition 5.5. Suppose that p" satisfies K-PI and the data distribution Pt is smooth. Denote
F(pr. 10) prl VallVlog log p" 12
then the following bound holds for any t > to
TV(pellp') I 8 (2)F((p))) O(a)
(15)
Relating the gradient of KI. divergence and the Fisher information and letting o - O, we have Theorem 5.6. Under assumption 5.4 and the assumptions in Proposition 5.5 suppose also that the KL divergence (Pt (pt,ll""") < cc. The following bound holds
min TV(p: s SO(T co). t
(16)
Please refer to Appendix
B.2 and B.3
for detailed proofs.
6 Experiments
In this section, we compare CFG to other bascline methods for constrained domain sampling, including MSVGD (Shi et al. 2022), MIED (Li et al. 22022), Spherical HMC (Lan et al. (0014), PD-SVGD and Control SVGD Liu et al. 202T). Throughout the section, we use neural networks with particles with Gaussian distributions unless otherwise stated. We All the experiments were implemented with Pytorch. More details can be found in Appendix [ The code is available at hhttts::/gggtee
2 hidden layers and initialize
use the neural network structure designed as in Appendix [
com wwwpeeehhhanonoonnnaaaiieee
6.1 Toy Experiments
We first conduct 2-D experiments to test the effectiveness of CFG sampling from truncated Gaussian distributions within various constrained domains, including ring-shaped, cardioid-shaped and double-moon-shapee domains. We covered non-convex domains, including unconnected domains. use 1000 particles for all domains. From the left of Figure [ we can see that the particles quickly converge to the target distribution without escaping the constrained domain. Comparison to MIED (via various distributional metrics) on the first three domains are deferred to Appendix D.I where CFG provides better approximation in most cases.
To compare with MSVGD and MIED, we additionally conduct the experiment of sampling from truncated gaussian mixture distribution within the block-shaped domain. The initial distribution is the uniform distribution. The right of Figure shows that our method outperformed MSVGD and achieved comparable results to MIED in terms of Wasserstein-2 distance and energy distance.
6.2 Bayesian Lasso
The Lasso method is broadly used in model selection to avoid overfitting by imposing a penalty on model parameters. Park Casella (2008) introduced a Bayesian alternative, named Bayesian Lasso, that replaces the plenty term with the double exponential prior distribution Pprior (3) oc Lan et al. (2014) then proposed Spherical HMC method, which introduced more flexibility choosing priors with explicit l -nom constraints: Poeior (3) (11)11 For
1, it corresponds to the Lasso method. For more general case when > 1, it is called Bayesian Bridge regression.
Following Lan et al. 2014), we choose the prior to be the truncated Gaussian distribution, where p(3) (0, II). This leads to the Bayesian regularized linear regression model (Lan et al. 2014):
02 ~ 310² o
1)1 <))
(17)
7

---

2200
Wasserstin Distancce (Sinkhorn)
MSVGD
CFG
0.35
Ean
K40 220
Far 2200
Ground Iats
0.30
0.25
far
0.045
0.040
CFG
0.030
0.025
Ben.
K40 200
Iar 2000
Grunnd Fate
Ereegg
0.020
0.020
0.005
100
200
400
Iter
Figure 1: Left: CFG sampled particles at different numbers of iterations constrained domains (ring, cardioid, double-moon, block). Right: The convergence curves of MSVGD, CFG and MIED on the block constraint.
The posterior distribution for 8 is Aly,X,ot o2 + +1)-1)1(1111, S r), where 3 = (I" This posterior has an inequality constraint and our method can apply.
6.2.1 Synthetic Dataset
first generate a 20-dimensional dataset (X, y), where X E x000 y E 1000 and X, + e, where 0, o....... 0), where the first half of components is tens and the second half is zeros. Let GOLs denote the estimates obtained by ordinary least squares (OLS) regression. We follow Park & Casella 2008) and set = [| ||1-
(0, 25I). We set the true regression coefficients to be Struc
We compare method Wasserstein-2 distance and energy distance with Spherical HMC and MIED using different numher of particles. The ground truth is obtained via rejection sampling using 100,000 posterior samples. We use the hodd) scheme to adapt the bandwidth to the number of particles N (see Appendix C. where we choose hgd 3 0.1. From Figure 2 we see that CFG performs the best in both metrics when N is small. This indicates the sample efficiency of particle-base brtdcle-based variational inference for constrained sampling, which aligns with the findings in Liu & large.
Wang (2016) for unconstrained sampling. All methods provide similar results when
It is worth noting that as kernel- based method, the time complexity of MIED scales quadratically as N increases, while functional gradient methods like CFG scale linearly. Meanwhile, we observed in our experiments that MIED tends to require more iterations to converge when the number of particles increases. These issues add to the overall time cost for MIED to achieve accurate approximation. especially when large N is used (see Figure Appendix D.2). Similar issues of SVGD are also stated in Dong et (2023). Please refer Appendix [D.2 for detailed information.
6.2.2 Real Dataset
Following Lan et al. 2014). we also evaluate our method using the diabetes dataset discussed in Park Casella (2008). We compare the posterior median estimates given by Spherical HMC, CFG and
MIED for the Lasso regression model (q 1) and the Bridge regression model (g 1.2). as the
shrinkage factor r/||BoLs||, varies from 0 to We 5000 particles for CFG and MIED.

---

1.0
CFG
0.8
Wassersteins Distance
MIED
0.6
0.6
Energy
0.4
0.2
0.0
10
10a
10a
10
10a
Sample Size (N)
Sample Size (N)
Figure Left: W distance of SPH, CFG and MIED versus the number of particles, Right:
Energy distance of SPH, CFG and MIED versus the number of particles. Both synthetic dataset.
Aia
0.0200
0.0755
0.0150
Brigg Conffcinnts
00075
0.050
0.0025
MED 0.0000
5PH
C7C
Figure 3: Left: Bayesian Lasso (g 1) using Spherical HMC (upper left), CFG (upper middle) and MIED (upper right). Bayesian Bridge Regression (g = 1.2) using Spherical HMC (lower left) CFG (upper middle) and MIED (upper right). Right: Results of monotonic Bayesian neural network with e
0.01. Only the portion below 0.02 shown on the y-axis to better display the performance of models satisfying constraint.
From Figure 3 we can see that our method aligns well with Spherical HMC and MIED, which shows the effectiveness of our method.
6.3 Monotonic Bayesian Neural Network
We use the COMPAS dataset following the setting in Liu et al. (0021) Given a Bayesian neural network 0). define the constraint function g( = m (8) with (m)
EgwDill (en-D[l(--0...... 100.00999)); 0)) 1l1]. Here, Tmeso denotes the subset of features which the output should be monotonic. Note that in Liu et (202)), the constraint is only in the sense of expectation, i.e.,
E, [0)] < 0. the context of monotonie neural network (Karpf 1991] Liu et al. 2020), however, the monotonicity constraint is (0) < 0, which defines domain with inequality constrains and thus our method can be applied.
We use two-layer ReLU neural network with 50 hidden units. With different threshold E {0.005, 0.01, 0.05). we compare our CFG with MIED and SVGD variants in Liu al. (2021): PD SVGD and Control SVGD (C-SVGD). The number of particles is 200 and the results are shown Table We observe that in all the tasks, our method outperforms the other two baselines in terms of both test accuracy and test likelihood. "Ratio Out" denotes the proportion of particles outside Q= (0( S 0} to the total number particles. All methods can successfully force almost all particles into the domain. We further plot in the right of Figure 3 the averaged test monotonic loss Lmono against test accuracy during the training process for 0.01. For MIED, the test accuracy higher when more particles outside the domain, while lower when forcing the particles into the

---

domain during training. Notably, Vanilla SVGD exhibits lower accuracy and higher monotonic loss, indicating the importance of constraint samplingg methods.
Furthermore, our method can scale up to even higher dimensional real problems. We experimented on the COMPAS dataset using larger monotonic BNNs and on the 276-dimensional larger dataset Blog Feedback Liu et al. 22200). Please refer to Appendix D.3 detailed information.
Table Results of monotonic Bayesian neural network under different monotonicity threshold. The results are averaged from the last 10 checkpoints for robustness.
TEST ACC
TEST NLL
RATIO OUT (5)
PD-SVGD
C-SVDD
MIED
PD-SVGD
C-SVGD
MIED
PD-SVGD C-SVGD MIED
CPG
0.0
J104
635
678
0.
0.005
005
650
595 800 .850 000
636
.33
676 001
.632 .000
00
6.0
7
Conclusion
We proposed new functional gradient ParVI method for constrained domain sampling, named CFG, which uses neural networks to design piece-wise velocity field that satisfies non-escaping boundary condition. We presented novel numerical strategies to deal with boundary integrals arising from the domain constraints. We showed that our method has TV distance convergence guarantee. Empirically, we demonstrated the effectiveness of our method on various machine learning tasks with inequality constraints.
Acknowledgments and Disclosure of Funding
This work was supported by National Natural Science Foundation of China (grant no. 12201014 and grant no. 12292983). The research of Cheng Zhang was supported 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. The authors are grateful for the computational resources provided by the High-performance Computing Platform of Peking University. The authors appreciate the anonymous NeurIPS reviewers for their constructive feedback.
References
Kwangjun Ahn and Sinho Chewi. Efficient constrained sampling via the mirror-langevin algorithm. In A. Beygelzimer, Y. Dauphin, P. Liang, and J. Wortman Vaughan (eds.), Advances in Neural Infor- mation Processing Systems, 2021.
David Alvarez-Melis, Yair Schiff, and Youssef Mroueh. Optimizing functionals on the space of probabilities with input convex neural networks. Trunsactions on Machine Learning Research, 2022. ISSN 2835-8856.
D. M. Blei, A. Y. Ng, and M. I. Jordan. Latent dirichlet allocation. allocaationn Journal of Machine Learning Re- search, 3(am)):933-1022 2003. hattttts::///ml> csail http:llmmitttueueeee mmt..eeu/ppaaerss/3//lee...
html
David M. Blei, Alp Kucukelbir, and Jon D. McAuliffe. Variational inference: A review for statisticians. Journal of the American Statistical Association, 1 112:859 877, 2016.
Nicolas Brosse, Alain Durmus, Éric Moulines, and Marcelo Pereyra. Sampling from a log-concave distribution with compact support with proximal langevin monte carlo. In Conference on learning theory, pp. 319-342. PMLR, 2017.
Mareus Brubaker, Mathieu Salzmann, and Raquel Urtasun. A family of meme methods on implicitly defined manifolds. In Artificial Intelligence and Statistics, pp. 161-172. PMLR, 2012.
Sébastien Bubeck, Ronen Eldan, and Joseph I Lehec. Sampling from a log-concave distribution with projected langevin monte carlo. Discrete & Computational Geometry, 59:757-783, 2018.
10

---

Simon Bymne and Mark Girolami. Geodesic monte carlo on embedded manifolds. Scandinavian Journal of Statistics, 40(4):825-845, 2013.
Changyou Chen, Ruiyi Zhang, Wenlin Wang, Bai Li, and Liqun Chen, unified particle-optimization framework for scalable bayesian sampling. ArXiv, abs/ 805. 1659, 2018.
Tianqi Chen, Emily Fox, and Carlos Guestrin. Stochastic gradient hamiltonian monte carlo. In International Conference on Machine Learning, pp. (683-1691, 2014.
Z. Cheng, S. Zhang, L. Yu, and C. Zhang. Particle-based variational inference with generalized wasserstein gradient flow. Advances in Neural Information Processing Systems, 2023.
S. Chewi, M.A. Erdogdu, M. Li, R. Shen, and S. Zhang, Analysis of langevin monte carlo from poincaré to log-sobolev. Proceedings of Thirty Fifth Conference on Learning Theory. 178:1-2, 2022.
Sinho Chewi. Thibaut Le Gouic, Chen Lu, Tyler Maunu, and Philippe Rigollet. Svgd as a kernel- ized wasserstein gradient flow of the chi-squared divergence. Advances Neural Information Processing Systems, 33, 2020.
Thomas A. Courtade. Bounds on the poincaré constant for convolution measures. Annales de l'Institut Henri Poincaré Probabilités et Statistiques, 56:566-579, 2020.
Lauro Langosco di Langosco. Vincent Fortuin, and Heiko Strathmann. Neural variational gradient descent. arXiv preprint arXiv:2107.10733 2021.
Hanze Dong. Xi Wang, LIN Yong, and Tong Zhang. Particle-based variational inference with preconditioned functional gradient flow. The Eleventh International Conference on Learning Representations, 2023, URL htttps://pperre
S. Duane, A. D. Kennedy, B J. Pendleton, and D. Roweth. Hybrid Monte Carlo. Physics Letters B, 195(2):216 - 222, 1987.
Jiaojiao Fan, Qinsheng Zhang, Amirhossein Taghvaei, and Yongxin Chen. Variational wasserstein gradient flow. In International Conference on Machine Learning, pp. 6185-6215. PMLR, 2022.
H. Federer. Geometric measure theory. Springer, 2014.
Will Grathwohl, Kuan-Chieh Wang, Jörn-Henrik Jacobsen, David Duvenaud, and Richard Zemel. Learning the stein discrepancy for training and evaluating energy-based models without sampling. In International Conference on Machine Learning, pp. 3732-3747. PMLR, 2020.
Tianyang Hu, Zixiang Chen, Hanxi Sun, Jincheng Bai, Mao Ye, and Guang Cheng. Stein neural sampler. arXiv preprint 1rXiv:1810.03544 2018.
Michael I. Jordan, Zoubin Ghahramani, T. Jaakkola, and Lawrence K. Saul. An introduction to variational methods for graphical models. Machine Learning, 37:183-233, 1999.
Richard Jordan, David Kinderlehrer, and Felix Outo The variational formulation of the fokker-planck equation. SIAM journal mathematical analysis, 29(1):1-17, 1998.
Jørgen Karpf. Inductive modelling in law: example based expert systems in administrative law. In Proceedings of the 3rd international conference on Artificial intelligence and law, pp. 297-306, 1991.
S. Lan, B. Zhou, and Shahbaba. Spherical hamiltonian monte carlo for constrained target distributions. The 31st International Conference Machine Learning, 32:629-637, 2014.
N. F. Lewis, H. Nagarajan, and B. O. Palsson. Constraining the metabolic genotype-phenotype relationship using a phylogeny of in silico methods. Nature Reviews Microbiology, 10(4):291-305 2012.
Lingxiao Li, Qiang Liu, Anna Korba, Mikhail Yurochkin, and Justin Solomon. Sampling with mollified interaction energy descent. arXiv preprint arXiv:2210.13-000 2022.

---

Chang Liu and Jun Zhu. Riemannian stein variational gradient descent for Bayesian inference. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 32, 2018.
Chang Liu, Jingwei Zhuo, Pengyu Cheng, Ruiyi Zhang. Jun Zhu, and Lawrence Carin. Understanding and accelerating particle-based variational inference. International Conference on Machine Learning. PMLR, pp. 4082-4092, 2019.
Qiang Liu. Stein variational gradient descent as gradient flow. Advances in neural information processing systems, 30, 2017.
Qiang Liu and Dilin Wang. Stein variational gradient descent: general purpose bayesian inference algorithm. Advances in neural information processing systems, 29, 2016.
Xingchao Liu, Xing Han, Na Zhang, and Qiang Liu. Certified monotonic neural networks. Advances in Neural Information Processing Systems, 33:15477-15438, 2020.
Xingchao Liu, Xin Tong, and Qiang Liu. Sampling with trusthworthy constraints: A variational gradient framework. Advances in Neural Information Processing Systems, 34:23557-235663 2021.
B. M.rris. Improved bounds for sampling contingency tables. Random Structures Algorithms, 21(2):135 146, 2002.
Radford Neal. MCMC using hamiltonian dynamics, In S Brooks, A Gelman, G Jones, and XL Meng (eds.), Handbook of Markov Chain Monte Carlo, Chapman & Hall/CRC Handbooks of Modern Statistical Methods. Taylor & Francis, 2011. ISBN 9781420079425. URL http://books
google. cot../ocsooocoooooooooooooooco dfsaoiiIA
Park and G. Casella. The bayesian lasso. Journal of the American Statistical Association, 103(482):
681-686, 2008.
Andrey Pilipenko. An introduction stochastic differential equations with reflection, volume 1. Universitätsverlag Potsdam, 2014.
G. O. Robert and O. Stramer. Langevin diffusions and metropolis-hastings algorithms. Methodology and Computing in Applied Probability, 4:337-357, 2002.
Adil Salim and Peter Richtárik. Primal dual interpretation of the proximal stochastic gradient langevin algorithm. Advances in Neural Information Processing Systems, $3:3786-3705 2020.
Jiaxin Shi, Chang Liu, and Lester Mackey. Sampling with mirrored stein operators. In International Conference on Learning Representations, 2022. UtL.pppp::/pppaaerreee.ee////e
eMudnJsb115
Elias M. Stein and Rami Shakarchi. Real Analysis. Princeton University press, 2005.
I. Thiele, N. Swainston, R. M. Fleming, A. Hoppe. S. Sahoo, M. K. Aurich, H. Haraldsdottir, M. L. Mo, O. Rolfsson, M. D. Stobbe, and et al. community-driven global reconstruction of human metabolism. Nature Biotechnology, 3 31(5):419-423) 2013.
M. Wainwright and M. I. Jordan. Graphical models, exponential families, and variational inference. Foundations and Trends in Maching Learning, 1(1-2):1-305, 2008.
Max Welling and Yee W Teh. Bayesian learning via stochastic gradient langevin dynamics. In International Conference on Machine Learning, pp. 681-688, 2011.
Kelvin Shuangjian Zhang, Gabriel Peyré, Jalal Fadili, and Marcelo Pereyra. Wasserstein control of mirror langevin monte carlo. In Conference Learning Theory, pp. 3814-3844 PMLR, 2020.
Ruqi Zhang, Qiang Liu, and Xin T. Tong. Sampling in constrained domains with orthogonal-space variational gradient descent. Advances in Neural Information Processing Systems, 2022.
12

---

A Estimation of Boundary Integral
The boundary integral Ja pc- rid.S a crucial factor that distinguishes ParVI on constrained domain from unconstrained case, due to its computational intractability in general. We derive heuristic band-wise approximation (12) and thus the remaining issue is identify which particle is in the band-wise area {d(ar 02) - h
Proposition A.1. For any h > O, the sufficient condition for {a Q d(x 000 S h2 is
{ g(x)= S 0.
(
9 I + h [V 2 O.
Vo(x)
(18)
Vg(a) Proof. For any Ωt define f(r) 9 [[ Then f(-) is continuous with f(0)
g(x) & O. If f(h) 2 O, by intermediate value theorem, there exists hg E [0, hi] such that f(h)) 0. And we have either g(x)
g(ze) where zo
ho Vg(z) Therefore æ E 80 or [aa(a]]l
2o E 8N, implying dxx, a9) < h since || zoll = hn I h.
In fact, if g(x) S then d(x, 80) S equivalent max ((y) 0. When h is small and Vg(a)) * 0, the maximum point of LHS is approximately = x + h [[Vg(()]] steepest ascent direction. Hence (18) is reasonable proxy when is small.
Vg(a) since Vg(x) is the
B Proofs
B.1 Proof of Theorem 5.2
Theorem B.1. Let Mo max(g( < oo, there exists finte to S
Ma XO
such that g(xx)) S
0, Vt to.
Proof. It suffices to prove that for any Zp, there exists a finite to <
2(z XC
such that g(2,) S 0, Vt 2 to.
Since
Je g(()) = (Vg)²
=Vg)² hnne [g<0)
Al|Vg|| [g20)
Using assumption || C, we can state that there exists finite to S s(zo) 0, vt 2 to-
such that g(xx) s
XC
Firstly prove that there is a lo
g(zn) XC
such that g(2()) S 0.
By contradiction, if such to does not exist, we have g(2)) > 0 for all t [o,
o(zo) AC then dg
-/ VglV integrate both sides yields g(xx) g(20) = v(eg) contradiction.
A||Vgl|dt s -((z)). then g(2)) S 0.
Then prove that g(2)) S 0, Vt 2 to
not, there exists a minimum t > to. such that g(x) > 0. then d = -Al[Vgll < 0, which
means for sufficiently small 5 > 0. (2(-5) > O, contradiction.
B.2 Proof of Proposition 5.5
Proposition B.2, Suppose that pd satisfies R-PI and the data distribution Pt is smooth. Denote F(p., f*) -o Sp, dr, then the following bound holds for 1, In
TV(p:ll")) s V8(K + (F)(
(19)
13

---

Proof. By equivalent expression of TV distance, it suffices to prove that for any function ff satisfying
[fl 1.
gg [f] Ep. |/]l +
Note that for > t0. pr (92) =
(o50))"(2) 1 (100spp((t)((()dd + 11020000 (x)
BBp [5]
Ep. [f] =
/ (o p((z))f(z)d
(0000))(0))
/ 1000 (x) (x)))) +
where e2 = [ 1920) (xr
(w) pr (z)))((r)dd O(a) by the L, convergence of Gaussian convo-
lution Stein & Shakarch (2005).
The first term can be bounded by
( (ALosoo (f*(x) p((x))
-
-( ((1os50)0( (z)
(x)((())
Pc(x)
1)f(x))d
(x)x
/ (ososE (1)(
P((x)
1)2 / (1))
P((x)
1) (y))
(x(x)
(x(xx)
2 / 11050)0 (0)
pc(x)) 1)2 / [Los0)00())) (z) Pc(a))
(x) + +1ddr
(x(x)
$ 4 / 10000) (1)(
p+(
1)2
(x)
$ 88(1 - / |1000) (P(()))a)(
8(1 / ((n(o)))
8(1 b) with f 10 (0))
pr(x) dx 1.
(x(
Then note that
pp(x)) varg" V (x))
PP((a) for (x)
Eg-
82 b2.
From Courtade, (2020), we have the Poincaré coefficient of N' o²1) coefficient of pa p" [0 o² is no larger than x + o².
of(0,0²1) !) is o², and the Poincaré
Therefore, by the Poincaré inequality, we have
Pr(x) {T(()) < (k (m) / | f
Pe(x) 85
1 1 b S (1 - b2) varg. V
dr
= =m= / ((000pp((=))ss>p 8p- (2)
=
And
Ex. [f] -1p, [fll S V8(K + F2)F E2
14

---

with - Ol)).
B.3 Proof of Theorem 5.6
Theorem B.3. Following the same assumptions in Proposition 5.5 Suppose also that the KL
divergence (KI(Pto:PP oo, the following bound holds
(20)
Proof. By equivalent expression of TV distance, it suffices to prove that for any function f satisfying I/l 1,
min EE - Ep, [f]l OT cb)
rST
now only consider the time after lo. which means that pr. (?2) =
First we state the relation between 6r = (O I² and p=(x). When o + 0, by the pointwise
convergence of the convolution Stein & Shakarchi (2005), we have Vf* + Vp" (z) and p" p" (x). Note that for g(x) S 0, pd p"(x)) and p" (²(x) are strictly positive, thus sp- (x) - Sp (z), as - for any x E Q
Note that 5* (I:) is supported on Rd, we have
d
D DKL (pt dt
=-f pov' pavE(opp) (8p0 8p. Jdz
/ polbnne [
Vg
| (020)) ( Sp, )d
[ Polbbet ((050)) Sp, )dx spx / pil(hnet + A [V ((-=0)) ) (sppp Sp, )d
Vg
- / (650)]l180 8pp 12xx - Pr. 11000) - (sp- - (p))) 8pp )dx
Pt: (050) sp. 2d + e0 / Pr: 1050011880 Spp 12 + 4co / pt (0500 Mn (89- Sp, )l
$ (1 - Gp) / 1 sp, ddr 200 S 1 ||nent - (8p sp.) 2d
+ 2en
P: 11.115001lpp 11sp sp 12
For 60 < 1, the second term is O(e), and the third term is O(a).
Using the notation F(pe.p") f Ps- (g50)) pl Sp. 1l integrating both sides yields the following:
[' F(De.8')dt & (' F(pr.f')dt 1-
So
min F(p.,p") 6SST
D
(21)
(1 (1-
15

---

Plugging the result of P()
21 in Proposition 5.5 and let - O, we have
min Bp [/] Ep, [fll s OCT +
Finally, by triangular inequality:
Epp [f] - Ep. [f]| 1B [f] Ep, [f]] E [f] - Eg. [f]]
Using the result that Gaussian convolution converge original L1 function in L1 (Stein & Shakarchi 2005), we have Egp [f] Ep [f]l + when + 0 for any p'. This means
mist 1 [f] Ep [f]l O(T 1 + e the proof is complete.
C Insights on Bandwidth Selection
Error Analysis of Boundary Integral Estimation We first present the error analysis of band-wise approximation, based on which we derive a heuristic O((dN)=)) scheme for bandwidth selection when adjusting the number of particles.
1 Σ v(os))"V((,,)//V(((e)) - pe.ppd.. rdS. We claim Denote the approximation error e mi
that the mean square error E[e?) = O(h + H)
The idea is leverage bias-variance decomposition.
estimate the bias of the basd-wise approximation, using Lagrange's Mean Value Theorem
1 : T. [ p(z)o vde
/ po pdS
A1
on
1
[/ [ p(z)u adS / p(z)u adS dl
(22)
-0 /-" ldl)
=(()..
Then we turn to bound variance. From 12 we have A Ss, p(z)o-ndr = plaeSa) And
Ss. ((z)u(z)) üdr..
h
u(z,)) nxxx) p(x E Sh))
J p(z)u(z) vdr
A2 = mh
h
Ø=1
mh (x;)) 7 7(xz)) mh S.
p(z)v(x) vdr
(23)
mh p(x)v(x) Adr p(x Sn) / p(x))(()) rdz Sh
Sn
A3 + A4-
The variance of Monte Carlo estimation is |[[A]] = (
By Central Limit Theorem, E[(" "h - p(2 $1 = O and thus E[A O(2).
Finally, through bias- bias-variance decomposition, have
O(2² + 1 1
E E
E[A + A3 = OEEA + A AID O(122 + n m
(24)
n
16

---

-0.5
-28
-1.
-20
-10
2.5
-2.
-10
-10
-35
10
10'
304
30
125
10
307
205
207
19
10"
10a
30"
Semple
Sump
Figure 4: Left: MSE of boundary integral estimation of distribution p1- Middle: MSE of boundary
integral estimation of distribution p2. Right: MSE of boundary integral estimation of distribution p3-
However, n depends on h implicitly. We use a simple example to determine this relation. Consider the target distribution being the uniform distribution, and the constrained domain being the ball in Rd with radius R, then p(a E Sr) x R RR-(R-A)" (dh)). In this case, the total mean square error of order O + ZNT which implies the optimal h is O((dN)-5) in this case the total mean square error is of order ((ddN) 3).
R
In practice, one can set h = (dN)) and tune hy through sparse grid search.
Simulation Verification of Boundary Integral Estimation To verify the above analysis, we conduct toy 2D simulation study estimating boundary integral of 3 different velocities at the boundary using different numbers of samples from different distributions. The constrained domain is the block area f {x |z1|| < 2, |22| 2}.
The number of particles are 10?, 103, 10 105, 10 The types velocity and distribution are listed in Table 2 The corresponding true values of boundary integral Jan pv rdS are listed in the Table
Table The types of velocities and distributions.
Types 2
3
ñ p(r) Unif((?)
v(x)
(1))
(a) x
O)(0,0),1) on
N((0--1 1) on
Table 3: The true values of boundary integral in the verification experiment.
Velocities
Distributions
1/2
V3
Pi
I
0
0
P2
0.226259
P3
0.911333
0
-0.617187
We estimate [E[e"] by 10 trials for difffernt number of particles and plot the curves as in Figure 4 We set ho (d and hod 0.5. The slope of the log-log plot is approximately - for all types of velocities and distributions, which supports the previous analysis.
Additionally, if we do not follow the scheme of h (oN) 1, and for example, fixing 0.5 (102)= or = 0.5 (10%)- - (the starting and ending edgewidth of the adaptive edgewidth scheme h = 0.5 N 4).. Figure 5 shows that the estimated B[e²] will tend to go up in certain cases, which supports the rationality of the h ho(dN)- scheme.
17

---

-0.5
-1.0
(00L((MSE)
-155
-2.
-2.
10a
10a
10
10a
10a
Sample Size (N)
Figure MSE of boundary integral estimation of distribution Pi and velocity ii using fixed edgewidths and adaptive edgewidth.
D Additional Details of Experiments
Structure of Neural Networks Fnet as follows
parameterize Anet (10) using two neural networks faet and
Unne = Snet "net Vg.
(25)
where fnet : Rd R2, Znet Rd + R. Here, the -2 act Vg term relates to reflection process in reflected stochastic differential equation Pilipenko regarded as a reflection term that forces the particle to stay inside Q as it neutralizes the outward Normal--ointing component. Heuristically. Vg serves the reflection direction and = contributes in learning the magnitude. Empirically, we observe that incorporating the gradient of the boundary Vg substantially alleviates the issue of particles adhering the boundary during training, and may slightly improve the sampling quality. Please refer to D.1 below.
equationPPllipenke (2014), and always points inward Q. It can be
D.1 Toy Experiments
Setting Details On toy experiments, we conduct CFG on four 2-D constrained distributions. The experiments are implemented on Intel 2.30GHz CPU with RAM 16384MB and NVIDIA GeForce RTX 3060 Laptop GPU with total memory 14066MB. The density and constrained domain are in Table The truncated gaussian mixture is 9 equally weighted gaussian distribution with all standard variance equal to 0.2, the centers are (-1.7,-1.7) (-1.7, 0), (-1.7, 1.7), (0, -1.7), (0, O), (0, 1.7), (1.7, -1.7), (1.7, 0), (1.7, 1.7).
Table 4: Four 2-D constrained distributions implemented in the toy experiments.
Name
Density
Constrained domain
Ring Cardioid
p" (x) xc N0, , I) on Ω p² (x) x (0, /) on Ω
Ω x S |
< 4]
Ω
{x
x + (5 Φ
2/3,2 < <4
Double-moon
p (x) x on
Ω
{x]]
logg(1))
2,
20g 3) 2(11 18)2
where ((x)
-(||)
Block
Truncated gaussian mixture on Ω
Ω
x
[x
<
and 322
S 2}
For CFG, fnet and Znet are three-layer neural networks with LeakyReLU activation (negative slope=0.1). The number of hidden units 128, except for the RING on which 256 is used. Both neural nets are trained by Adam optimizer with learning rate 0.002, except for RING on which 0.005 used. The number inner-loop gradient updates for fnet and znot is set 10, except for the RING on which 3 is used. The total number of iterations is 2000 and the step size of particle is 0.005 except for the RING on which 0.01 is used. The band width is set 0.05 except for BLOCK on which 0.001 is used. À in the picce-war picce-wise construction of the velocity field is chosen to be I.
18

---

On the BLOCK experiment, the total number of iterations is 2000. The ground truth is obtained by rejection sampling using 10t posterior samples. We run the accuracy experiments on 3 random seeds. The average results and the standard errors of the means are represented in the figure using lines and shades. For MSVGD, we follow Shi et al. (2022) and set the learning rate to 0.05 (selected from {0.01, 0.05, 0.1)). The kernel is the IMQ kernel and the kernel width is 0.1 (selected from {0.05, 0.1, 0.2)). For MIED, we use the default hyperparameters in Li et al. (2022) and set the particle stepsize to be 0.01 (selected from {0.005, 0.01, 0.02)).
Additional Comparison MIED Table shows the additional comparison of CFG and MIED on the first three constrained domains. CFG achieved better results in most of the cases.
Table 5: Wasserstein-2 distance and energy distance between the target distribution and the variational approximation on the three toy datasets.
Wasserstein-2 distance (Sinkhorn)
Energy distance
Name
MIED
CFG
MIED
CFG
Ring
0.1074
0,1087
0.0004
0.0003
Cardioid
0.1240
0.1141
0.0016
0.0005
Double-moon
0,4724
0.1660
0.0629
0.0022
Additional Ablation Studies From Fig Table 6 we can see that, without estimating the bound-
ary integral leads to failure sampling in most cases. This indicates the essential difference between constrained and unconstrained domain sampling. Adding Enet can achieve slight improvement in avoiding clustering near the boundary and better KL convergence.
Table 6: Ablation results of not estimating boundary integral, with and without 2net-
Wasserstein-2 distance (Sinkhorm)
Energy distance
Ring
Cardioid
Double-moon
Block
Ring
Cardioid
Double-moon
Block
w/o boundary integral
0.2138
0.2321
0,4866
0,2438
0.0097
0.1147
0.0068
0.0073
w/o Znet
0.1248
0,2234
0.1217
0.2422
0.0013
0.0009
0.0049
0.0073
w/ Znet
0.1087
0.1660
0.1141
0.2416
0.0003
0.0005
0.0022
0.0072
Geometric Generalization of Accommodating Multiple Constraints Our proposed method can generalize and accommodate multiple constraints (including more equality and inequality constraints) and more complicated geometries.
For additional equality constraints, our proposed framework still apply by adopting the idea of Zhang et al. (2022). We could use velocity fields like (y(x) to guided the particles to satisfy the equality constraints, and use our method to propose velocity fields substituting (x) to satisfy the inequality constraints. Adding these two types of velocity presents the desired velocity field. To illustrate the idea, demonstrate the training process of 3D toy ring distribution example. Suppose the coordinate of the particle is x = (x, y, z). The target distribution is a truncated standard gaussian located in the ring shaped domain in xOy plane. This corresponds to the equality constraint
and the inequality constraint 1 < x² + 2² < 4. The initial distribution is the 3D standard gaussian distribution. We choose Tn (x)
-iign (z) 1.5, From figure 7 we can see that the particles collapse to the xOy plane and converge inside the ring domain.
For additional inequality constraints, need multiple velocity field outside the constrained domain for particles to enter the constrained domain, and the boundary integral term can be similarly estimated using band-wise approximation.
D.2 Bayesian Lasso
Setting Details The experiments are implemented on Intel 2.30GHz CPU with RAM 16384MB and NVIDIA GeForce RTX 3060 Laptop GPU with total memory 14066MB.
19

---

w/o boundary integral
w/o Znot
Znet
Ground Truth
-1
-2
2
0
Ground Truth
-3
2
2
-2
0
22
2
-2
0
Ground Truth
-2
0
2
4
2
0
2
4
0
4
2
0
2
4
Ground Truth
-1
-2
-2
2
2
Figure 6: Ablation sampling results of not estimating boundary integral (left), not using Znet (middle left), using znet (middle right) and the ground truth (right).
On synthetic dataset, we run the experiments on 5 random seeds. The average results and the standard enrors of the means are represented in the figure using lines and shades. The number of particles N are chosen from {5, 15, 30, 70, 180, 400, 900, 2000).
For CFG, Jnct and Znet are three-layer neural networks with LeakyReL.U activation (negative slope=0.1). The number of hidden units is 256. Both neural are trained by Adam optimizer with learning rate 0.0005 for 10 iterations in the inner loop. The total number of iterations is 1000 and the stepsize of partiele 0.004. in the piece- construction of the velocity field is chosen be
For Spherical HMC, we follow the same main setting as Lan et al. (2014). The number of burn-in epochs is 4000, and the number of leap-frog is selected from (50, 100
For MIED, we use the default hyperparameters in Li al. (2022) and set the particle stepsize to gradually decay when the particle number N goes up. We set particle stepsize to 0.5 for N E (5, 15, 30, 70), 0.35 for E (180, 400), 0.3 for N 900, and 0.2 for N 2000. The number of iterations is 4000 for convergence except for N 2000, which needs 10000 iterations converge.
20

---

Iter:
Iter: 100
Iter: 1000
Ground Truth
2010
2210
2210
3210 ub4o
Figure 7: Illustration of generalizing method accommodating equality and inequality con- straints.
020
3119
22125
00/0
0.7%
12
Figure 8: First row: Wasserstein-2 distance of SPH, CFG and MIED versus the number of particles on dimensions 5, 10, 15, 20, Second row: Energy distance of SPH, CFG and MIED versus the number of particles on dimensions 5, 10, 15, 20. Both on synthetic dataset.
On real diabetes dataset, for CFG, fnet and Znet are three-layer neural networks with LeakyReLU activation (negative slope=0.1). The number of hidden units is 50. Both neural nets are trained by Adam optimizer with learning rate 0.005 for 10 iterations in the inner loop. The total number of iterations is 300 and the bandwidth h = 1. The number of particles is 5000 and the stepsize of particle is selected from (0.9, 1.05, 1.2).
For Spherical HMC, we follow the same setting as Lan et al. 2014). The number of total epochs is 1 11000 and the number of burn-in epochs is 1000.
For MIED, we use the default hyperparameters in Li et al. 2022 and the particle stepsize is selected from (0.1, 0.2, 0.5). The number of iterations is 2000.
Additional Comparisons on the Synthetic Dataset We further compare the sampling results of SPH. CFG and MIED on different dimensions. We choose dimensions 5, 10, and 20. From the design of Strue. the samples of first two dimensions are close to 10, and the last two dimensions are close O. Figure & shows that CFG achieves comparable results to MIED on Wasserstein-2 distances, while slightly better result on energy distance in most of the dimensions. SPH eventually catches up when the particle number grows up. This alignes with the results stated Section 6.2.1
Larger number of particles is important for higher sample qualities. Similar to Dong et al. (203)) CFG is more scalable than MIED in terms of particle numbers N. This is because CFG obtained with iterative approximation, the complexity is O(N). On the other hand, MIED incurs complexity ( due to the calculation of the weight denominator Eig ett..
The Nonlinear Time Complexity of MIED
Figure 9 plot the Wasserstein-2 Distance versus the computation time using at most 1000 iterations for CFG and 10000 iterations for MIED. The number of particles are E (900, 2000, 4000). clear that MIED entails more computation cost with the increase of particle numbers for larger N.
21

---

1.2
1.0
Wasserstein-2
0.8
0.6
0.4
0
200
400
600
800
Figure 9: Wasserstein-2 distance of CFG and MIED versus the training time using different number of particles.
CFG(h=0.1)
CFG/h=0.01)
logualEnergy Distance)
CFGth-0.1N -)
10-1
10a
10a
10a
Sample Size (N)
Figure 10: Choosing the adaptive bandwidth (red) against fixed bandwidths for the Bayesian Lasso experiment on synthetic dataset.
From the ending points of each curves (denotes the computation time for 1000 iterations of CFG and 10000 iterations of MIED), time cost of CFG only grows linearly (from 130s to 280s to 600s approximately), while the time cost of MIED grows quadratically (from 100s to 220s to 830s approximately), which supports the above claim of time complexity.
Additionally, figure g indicates that for large N, the number of iterations for MIED should also increase. 4000 iterations is sufficient for N 900, while 10000 iterations is needed for N 2000, and more than 10000 iterations is needed for MIED to completely converge to achieve better sampling results.
The Effect of Choosing the Adaptive Bandwidth Figure compareed the energy distance between the adaptive bandwidth scheme = 0.1N- and fixed bandwidth = 0. 0.01 of CFG. We can see that the adaptive scheme achieved slightly better energy distance results than the fixed bandwidth.
D.3 Monotonic Bayesian Neural Network
Setting Details The experiments are implemented Nvidia GeForce RTX 4080 Laptop GPU with memory 12 GB. We use 200 particles and 200 data batch size for stochastic gradient in all methods.
For CFG, both fnet and žnet are three layers with LeakyReLU activation (negative slope 0.1). foet have 300 hidden units while =net have 200. Both neural nets are trained by Adam optimizer with learning rate 0)001 for 10 iterations in the inner loop. The bandwidth h 0.02 and the stepsize of particle is 5e-5. The total number of iterations is 200. the piece-wise construction of the velocity field is chosen to be 100.
22

---

For PD SVGD and Control SVGD, we run 3000 iterations with particle stepsize selected from {0.001, 0.002, 0.005} ensure convergence.
For MIED, we run 1500 iterations with default hyperparameters in Li et al. 2022) ) and particle stepsize selected from (2,3,5) 5) 10-1
Additional Experiments on Scaling Up to Higher Dimensions To test our method on higher dimension, first experiment the COMPAS dataset using larger monotonic BNNs. The number of neurons the layer of BNN increased from 50 to 100, making the particle dimension increase from 903 to 1502. From Table our method still achieved the best accuracy results compared to other methods, while achieving competitive log-likelihood results.
For even higher dimension, we additionally experiment on the 276 16--dimensionl larger dataset Blog Feedback Liu et al. (2020) using monotonic BNN. The particle dimension 13903. From Table 8 our method still achieved the best result.
Table 7: Results of a larger monotonic Bayesian neural network on COMPAS dataset under different monotonicity threshold. The results are averaged from the last 10 checkpoints for robustness. For each monotonicity threshold, the best result is marked in black bold font and the second best result is marked in brown bold font. Positive proportion of particles outside the constrained domain is marked in red.
TEST ACe
TEST NLL
RATIO OUT (5)
PD-SVGD
C-SVGD
MIED
CPG
PD-SVGD
C-SVG
MIED
CRG
PD-SVGD
C-SVGD MIED
CPG
0.05
818 L. 008
630
004
568 0000
849 1. 001
639 1 000
.630 1: .001
.665 000
857 1. 000
0.0
0.0
0.0
0.0
0.01
618 = 100
688
004
309
000
.651 + 000
630
1: 0022
.681 + .001
666 000
860 #: 000
13.0
0.0
0.0
0.0
0.005
818
000
638
001
571
001
.653
000
630 + 00)
.631 4 .001
664 +
857 + 000
24.0
0.0
0.0
0.0
Table 8: Results of monotonic BNN on a higher 776-dimensional dataset Blog Feedback. The particle dimension is 13903. The results are averaged from the last 10 checkpoints for robustness.
TEST RMSE
RATIO OUT (%)
PD-SVGD C-SVGD
MIED
CFG
PD-SVGD C-SVGD MIED
CFG
.205 + .011 .27 .000 .212 0001
.204 .00 .000
0.0
0.0
0.0
0.0
E
Limitations and Future Work
Extensions to SVGD. The main idea of our work is enforcing the particles into the constrained domain, and then using neural networks capture the target distribution. The latter part can also be accomplished by using kernels like in SVGD. The boundary integral may be treated in the similar way. We believe that our framework can be adapted to larger family of variational sampling methods.
More Efficiently Choosing Bandwidth h. this paper, we only consider the simple case of uniform target distribution in choosing better adaptive scheme for bandwidth h. Proposing adaptive scheme based on the current particles is an interesting topic, which will be left for future work.
Extensions to GWG. Cheng et al. (2023) proposed a variational framework including general geometries by using 4p norm regularization. In this case, the treatment of boundary integral is the same as in the RSD loss expansion. The effect of using other regularization terms is left for future study.
Relations with the Reflected Stochastic Differential Equation. Contrast to SDE-based sampling methods, our work is based on ODE. Both SDE and ODE can simulate trajectories including reflection. We proposed heuristic understanding of the relation between SDE and ODE reflection in this work, more thorough theoretical analysis may be presented. We leave this to future research.
23

---

Broader Impact
This paper presents work whose goal is advance the field of machine learning. There are many potential societal consequences of our work, none which we feel must be specifically highlighted here. As far as we are concerned, our paper has no potential negative societal impacts.
24

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: The last four sentences of the abstract and the last four sentences of the introduction reflect the paper's contributions and scope.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification: Please refer the appendix E
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create separate "Limitations" section their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't 't acknowledged judgment and recognize that individual actions in favor of transparency play impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
the paper. The authors should use their best
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
25

---

Justification:
Please refer to appendix B
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: Please refer to appendix D and the codes and toy datasets the supplemental materials.
Guidelines:
The answer NA means that the paper does not include experiments. If the paper includes experiments, No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how reproduce that algorithm. (b) If the contribution is primarily
model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility, In the case of closed-source models, it may be that access to the model is limited in some way (e.g., registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
26

---

Answer: [Yes]
Justification:
Please refer to the codes and toy datasets in the supplemental materials.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://niis.. cc/
for more details.
While we encourage the release of code and data, we understand that this might not be possible, "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (c.g., for a new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https guwde//ddeeeeee Cooeeeoodeesssseeeeeee //nips.com htte://bb///uudddd///eeueeiessss ssioislicy for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible supplemental material (appended the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (c.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary understand the results?
Answer: [Yes]
Justification: Please refer to appendix D and the codes and toy datasets in the supplementary materials.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment Statistical Significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification: Please refer to the right of figure I figure and figure
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the enror bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
27

---

It should be clear whether the error bar is the standard deviation or the standard error of the mean.
is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out range (e.g. negative еπхо rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: Please refer to appendix D and figure 9
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper).
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics
Answer: [Yes]
Justification:
We follow the NeurIPS Code of Ethics.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., there a special consid- eration due to laws or regulations in their jurisdiction).
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification: Please refer to appendix F As far as we are concerned, our paper has no potential negative societal impacts. We do not engage in negative fairness, privacy or security issues.
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
28

---

The conference expects that many papers will be foundational research and not tied particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms monitor how system learns from feedback over time, improving the efficiency and accessibility of MI..
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: Our paper poses no such risks for misuse.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Lieenses for existing assets
Are the creators or original owners of assets (e.g., code, data, models), used in the paper. properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification: We cite the original papers in the experiments. Please refer to Section 6
Guidelines:
The answer NA means that the paper does not existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (c.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (c.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi thcode. com com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
29

---

If this information is not available online, the authors are encouraged to reach out to the asset's 's creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification: Please refer the supplemental materials.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/moodl as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification: Our paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According the NeurIPS Code of Ethies, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: Our paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurlPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
30


Reviews:
[{"review_id": "15xwU4hhWb", "reviewer": "Reviewer_s6xw", "summary": "The paper introduces an extension of the ParVI method for sampling on measures supported on constrained subset $\\Omega$ of $\\mathbb{R}^d$, where $\\Omega$ is the sublevel set of a phase-field function $g$.     The paper presents a construction in terms of a discontinuous particle velocity field, with the discontinuity occurring on the boundary of the domain.   Outside the domain,   the velocity field pulls the particles towards the boundary of the domain, while within the domain the normal parVI approach is used.\n\nTo leverage the usual Stein formulation through integration by parts, a new boundary term occurs which must be calculated.   This is approximated by a spatial integral along a thin band within the domain.\n\nThe authors prove that (i) particles eventually arrive at the domain boundary (using the assumption that there are no stationary points in the level set function $g$ outside the domain);  and (ii) under a universal approximation assumption, the empirical distribution converges in TV, assuming a Poincare inequality holds.\n\nThe author demonstrate the method on a number of synthetic and real examples.", "questions": "(1) How limiting is the assumption that $\\lVert \\nabla g\\rVert \\geq C$?   I can see this being a reasonable assumption for simple constraints, but what about others where the boundary is complex, e.g. corrugated, or is non-convex?\n\n(2) Can the authors clearly explain the bias-variance trade-off in choosing $h$, and the implications of that for the performance of the method.\n\n(3) Can the authors explain why the particles do not 'stick' at the boundary when the step-size gets very small?  Similarly, can they explore the effects of the large discontinuity.   \n\n(4) For the problems considered in this paper, could the authors not simply instead use rejection sampling for a uniform distribution within the constrained domain and then use that as an initial condition for the constrained parvi?    Potentially then you wouldn't even need to consider the boundary condition?", "limitations": "Limitations have been adequately addressed.", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The paper presents an alternative approach to the mirror descent strategies to handling sampling from constrained domains.   This is potentially very useful in settings where there is not an obvious parametrisation of the constrained measure which lends itself to mirror-descent approaches -- to my knowledge this is novel.     The methodology is generally applicable, and the authors have been very clear about their assumptions and provided some good theoretical results.   The experiments are sufficiently challenging.", "weaknesses": "The main weakness to me lies in how the boundary term is handled, for two reasons.   My main concern is what happens when the step-size is very small.   In this case,  I can easily see scenarios where:\n\n(1) The particles reach the boundary and remain stuck there.   What is pushing the particles across the boundary.   When the step-size is not small this is ok, because particles will 'overshoot' and then get picked up by the other side of the velocity field.    This isn't really addressed in the text, nor is the effect of the discontinuity on the particle movements, which I expect to be non-trivial. \n\n(2) The handling of the boundary integral is challenging, because there is a clear bias-variance trade-off.  The smaller $h$ is, the less likely it is to find particles in the domain, and thus we would expect the variance of the integral estimator to be huge.    This isn't really addressed in the text.\n\n(3) I believe (though happy to be corrected on this), that most of the examples could be handled with the mirror approach?  It would have been nice to see examples where the alternative is simply impossible.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your insightful review and valuable questions! We address your comments and questions as below.\n\n### Weaknesses\n\n`W1`: \nMy main concern is what happens when the step-size is very small. In this case, I can easily see scenarios where:\n(1) The particles reach (...) non-trivial.\n\nA1: Thanks for the question! First, we want to emphasize that in practice the step-size needs not to be very small (the algorithm would still convergence inside the constrained domain as the functional gradient approaches zero). Assumption 5.1 provides positive velocity field on the boundary. Also, when particles get stuck near the boundary, the boundary integral part would provide extra \"force\" to push these particles towards the inner side of the domain if the local density of the particles is higher than the target. With proper step size, empirically the discontinuity on the particle movements does not affect the overall convergence. Please refer to Figure 6 and Table 6 (listed below) in appendix D.1 for empirical results.\n\n`W2`: The handling of the boundary integral is challenging, (...) This isn't really addressed in the text. \n\nA2: Thanks for the question! The bandwidth $h$ needs to be adjusted according to the number of particles $N$ to get good performance. To put it simple, $h$ should be large if $N$ is small and vice versa. The choice of $h$ needs to strike a good balance between the band-wise numerical approximation error of the integral and the variance due to Monte Carlo estimation. Please refer to appendix C for detailed error analysis and simulation verification. We will make this reference clearer in our revision.\n\n\n`W3`: I believe (though happy to be corrected on this), (...) simply impossible.\n\nA3: Thanks for the question! For non-convex and complex boundaries where the mirror map is infeasible, it is hardly possible to use the mirror method. We have experimented on these cases in our paper, including the ring, cardioid, double-moon experiments in section 6.1 and the monotonic BNN experiments. Additionally, when mirror method is applicable such as the block experiment, we have shown in the right of Figure 1 that we still achieved better results than the mirror method MSVGD.\n\n### Questions\n`Q1`:  How limiting is the assumption that $||\\nabla g||\\ge C$? I can see this being a reasonable assumption for simple constraints, but what about others where the boundary is complex, e.g. corrugated, or is non-convex? \n\nA1: Thanks for the question! Firstly, we would like to be clear that the main focus of our theoretical analysis is on the convergence behavior when the particles are inside the constrained domain. The assumption that $||\\nabla g||\\ge C$ is mainly to guarantee that the particles would enter the constrained domain. In fact, any assumption would work here as long as it is sufficient enough that the particles will not get stuck at local modes outside the domain. Secondly, this assumption is not very limiting, which is adopted from [1]. In general, if $g$ has local extrema outside the constrained domain which hinder the entering, the particles may stuck or collapse at these local modes and affect the training. This would be left for future work.\n\n[1] Zhang, R., Liu, Q., & Tong, X. T. Sampling in constrained domains with orthogonal-space variational gradient descent. Advances in Neural Information Processing Systems, 2022.\n\n`Q2`:  Can the authors clearly explain the bias-variance trade-off in choosing $h$, and the implications of that for the performance of the method.\n\nA2: Please refer to appendix C for detailed error analysis and simulation verification. We will refer to this discussion earlier in our revision.\n\n`Q3`:  Can the authors explain (...) large discontinuity.\n\nA3: Thanks for the question! Generally speaking, the additional boundary integral term in RSD will endow the trained velocity field with a repulsive factor that prevents clustering near the boundary unless the target distribution encourages the particles to do so. This is different from the unconstrained cases and demonstrates the importance of properly estimating the boundary integral. Empirically, for positive step size, we can see from Figure 6 in appendix D that the 'sticking' phenomenon does not occur.\n\n`Q4`: For the problems considered in this paper, could the authors not simply instead use rejection sampling for a uniform distribution within the constrained domain and then use that as an initial condition for the constrained parvi? Potentially then you wouldn't even need to consider the boundary condition?\n\nA4: Thanks for the questions! For your first question, for complex real problems such as the monotonic BNN, it may be sample inefficient to use rejection sampling for a uniform distribution simply to guarantee that the initial particles are inside the constrained domain. The constrained domain of some problems may be difficult to locate and may have a low acceptance rate for rejection sampling. For your second question, even if the initial particles are all inside the domain, the boundary condition is still essential for correct sampling of the constrained distribution. Please refer to appendix D.1 for detailed ablation study on with and without estimating the boundary integral (Table 6 in appendix D.1 is listed below). The latter suffers from the 'sticking' to the boundary problem as you stated before while the former does not.\n\nTable: Ablation results of not estimating boundary integral, with and without $z_{net}$.\n\n|Name|Ring(W2 distance)|Cardioid(W2 distance)|Double-moon(W2 distance)|Block(W2 distance)|Ring(Energy distance)|Cardioid(Energy distance)|Double-moon(Energy distance)|Block(Energy distance) |\n|----|----|----|----|----|----|----|----|----|\n|w/o boundary integral|0.2138|0.2321|0.4866|0.2438|0.0097|0.1147|0.0068|0.0073|\n|w/o $z_{net}$|0.1248|0.2234|0.1217|0.2422|0.0013|0.0009|0.0049|0.0073|\n|w/ $z_{net}$|**0.1087**|**0.1660**|**0.1141**|**0.2416**|**0.0003**|**0.0005**| **0.0022**|**0.0072**|", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Cf116AwrBf", "reviewer": "Reviewer_YhS1", "summary": "The authors study functional Wasserstein gradient descent methods with constraints in support of target distribution. They study the variational problem of vector fields using the least square formula between the current step and the gradient of score functions. To maintain the constrained set, the authors restrict the choices of vector fields to particular linear combinations of vector fields and use the neural network functions to approximate these vector fields. The authors prove the convergence of the proposed method in total variation. Numerical examples demonstrate the effectiveness of the proposed method.", "questions": "Can authors illustrate how to pick lambda? How does lambda affect the algorithm, at least for numerical reasons?", "limitations": "There are no limitations.", "rating": 6, "confidence": 5, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "1. The authors carefully analyze the algorithm's importance in handling the sampling domain's boundary. In particular, example 4.2 is a good example. \n\n2. The algorithm variational formulation is well-designed with nice numerical examples.", "weaknesses": "1. The authors construct (10) as the projected vector field. The choices of lambda seem to be lacking. \n \n2. The authors miss essential literature in the field. One needs to discuss the related references carefully. \n\n1. Wang, Y., Li, W. Accelerated Information Gradient Flow. J Sci Comput 90, 11 (2022). \n2. A. Lin et.al. Wasserstein Proximal of GANs. GSI, 2022. \n3. Wang, Y. et al. Optimal Neural Network Approximation of Wasserstein Gradient Direction via Convex Optimization. 2022.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your careful review and valuable questions! We address your comments and questions as below.\n\n### Weaknesses\n\n`W1`: The authors construct (10) as the projected vector field. The choices of lambda seem to be lacking.\n\nA1: Thanks for the suggestion! We choose lambda to be 1 in toy experiments and Bayesian Lasso experiments. We choose lambda to be 100 in monotonic BNN experiments. \nLambda is used to control how fast the particles can be pushed inside the constrained domain. According to the problem scale and dimensionality, we select lambda in a reasonable range. We will make this point clearer in our revision.\n\n\n`W2`: The authors miss essential literature in the field. One needs to discuss the related references carefully. \n\nWang, Y., Li, W. Accelerated Information Gradient Flow. J Sci Comput 90, 11 (2022).\n\nA. Lin et.al. Wasserstein Proximal of GANs. GSI, 2022.\n\nWang, Y. et al. Optimal Neural Network Approximation of Wasserstein Gradient Direction via Convex Optimization. 2022.\n\nA2: Thanks for the suggestion! These references implement Wasserstein gradient flows in multiple regions. Wang \\& Li (2022) considered the Accelerated Information Gradient flow, which is the generalization of the Hamiltonian flow integrating multiple metrics. They provided the particle version for the flow with detailed theoretical analysis, and they also experimented on common sampling tasks including BLR and BNN. Lin et.al. (2022) applied the Wasserstein proximal method to GAN, the underlying theoretical guarantee is about Wasserstein natural gradient flow. Wang et.al. (2022) considered a similar problem of using neural networks to approximate the Wasserstein gradient flow as [1], the main difference is that Wang et.al. (2022) considered parameterizing the velocity field using the gradient of a neural network, rather than directly using a neural network. All the above references do not consider the problem of constrained domain sampling, for which we designed functional gradient flows. We will add these references in our revision.\n\n[1] Zhang, R., Liu, Q., & Tong, X. T. Sampling in constrained domains with orthogonal-space variational gradient descent. Advances in Neural Information Processing Systems, 2022.\n\n\n### Questions\n\n`Q1`: Can authors illustrate how to pick lambda? How does lambda affect the algorithm, at least for numerical reasons? \n\nA1: Thanks for the question! Empirically, we pick lambda to be 1 for toy experiments and Bayesian Lasso experiments, and we pick lambda to be 100 for monotonic BNN experiments. Lambda is used to control how fast the particles can be pushed inside the constrained domain. According to scale of the problem, relatively large lambda may provide swift entering for some constrained domains, but may incur some instability if lambda is too large that particles overshoot the constrained domains.\nIn practice, for small scale experiments where the step size can be large, one may choose a small lambda. For large scale problems (e.g., BNN) where the step size is small, one may want to use a large lambda to drive the particles to the constrained domain faster.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "jLXFIADDex", "reviewer": "Reviewer_PQmA", "summary": "This paper addresses the challenging problem of constrained sampling from an unnormalized distribution. Building on the principles of Stein variational gradient descent (SVGD), which merges the strengths of variational inference and Markov Chain Monte Carlo (MCMC), this work innovates by applying SVGD in a constrained setting. Specifically, it introduces boundary conditions for the gradient flow to confine particles within a designated domain, demonstrating convergence in total variation. Numerical experiments highlight the method's effectiveness.", "questions": "If we have additional equality constraints, does the method still work?", "limitations": "The paper has discussed the limitations.", "rating": 7, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The paper is well-written.\n2. The proposed method is technically sound and demonstrates convergence in total variation.\n3. Experiments on both synthetic and real-world datasets showcase the method's potential.", "weaknesses": "In Bayesian neural network experiments, the paper overlooks numerous baselines of existing Bayesian Neural Networks (BNNs) using variational inference, MCMC, and Laplace approximation.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your careful review and valuable questions! We address your comments and questions as below.\n\n### Weaknesses\n\n`W1`: In Bayesian neural network experiments, the paper overlooks numerous baselines of existing Bayesian Neural Networks (BNNs) using variational inference, MCMC, and Laplace approximation.\n\nA1: Thanks for the advice! Existing methods are useful in common unconstrained problems, but these methods generally do not apply to problems with constrained domains. \n\n\n### Questions\n\n`Q1`: If we have additional equality constraints, does the method still work?\n\nA1: Thanks for the question! Our proposed framework can still apply to these types of problems by adopting the idea of [1]. We could use velocity fields like $v_{\\sharp}(x)$ to guided the particles to satisfy the equality constraints, and use our method to propose velocity fields substituting $v_{\\perp}(x)$ to satisfy the inequality constraints. We add these two types of velocity to obtain the desired velocity field. To illustrate the idea, we demonstrate the training process of a 3D toy example in Figure 2 of the attached pdf file. Suppose the coordinate of the particle is $(x,y,z)$. The target distribution is a truncated standard gaussian located in the ring shaped domain in xOy plane. This corresponds to the equality constraint $z=0$ and the inequality constraint $1\\le x^2+y^2\\le 4$. The initial distribution is the 3D standard gaussian distribution. We can see that the particles collapse to the xOy plane and converge inside the ring domain.\n\n[1] Zhang, R., Liu, Q., & Tong, X. T. Sampling in constrained domains with orthogonal-space variational gradient descent. Advances in Neural Information Processing Systems, 2022.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "oY0pBnCqXv", "reviewer": "Reviewer_2Jzp", "summary": "The authors develop a solution to constrained sampling by introducing a boundary condition for the gradient flow which would confine the particles within the specific domain. This gives a new functional gradient ParVI method for constrained sampling, called constrained functional gradient flow (CFG), with provable continuous-time convergence in total variation (TV). They also present numerical strategies to handle the boundary integral term arising from the domain constraints. They provide theoretical and experimental study to show the effectiveness of the proposed framework.", "questions": "There is no question from the reviewer.", "limitations": "N/A", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "The authors provide extensive theoretical analysis to support the proposed method. The authors provide solid theoretical results along with the proposed method. The authors also provide experiments on different datasets to validate the approach.  The paper is well written.", "weaknesses": "Experiments on a real-world application could strengthen the paper.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your careful review! We address your comments as below.\n\n### Weaknesses\n\n`W1`: Experiments on a real-world application could strengthen the paper.\n\nA1: Thanks for the suggestion! In our experiments, the monotonic BNN is applied to real-world problems. The setting is originally derived from [1] about modelling in law. Additionally, it is also widely applied in social problems, including predicting admission decisions ([2]). For higher dimensional real-world applications, please refer to the third part of the global response.\n\n\n[1] Karpf, J. (1991, May). Inductive modelling in law: example based expert systems in administrative law. In Proceedings of the 3rd international conference on Artificial intelligence and law (pp. 297-306).\n\n[2] Liu, X., Tong, X., & Liu, Q. (2021). Sampling with trusthworthy constraints: A variational gradient framework. Advances in Neural Information Processing Systems, 34, 23557-23568.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "fVrOzCZDWE", "reviewer": "Reviewer_SeKQ", "summary": "The paper proposed a particle based method (using neural networks) to sample probability densities  which are supported on subdomain of $\\mathbb{R}^d$.  This is done in the spirit of Stein variational gradient descent but using neural network instead of kernels and by incorporating the constraint into the functional which leads to a boundary integral which is handled via a band approximation. A complete convergence to equilbrium analysis is given, under a suitable assumption on the neural network approximation. The algorithm is tested on a variety of synthetic and real data sets and the results are competitive compared to other methods on the market (e.g MIED).", "questions": "1) Can the authors accommodate  multiple constraints and more complicated geometries? Are those generalizations possible? A discussion of these issues would be welcome.\n\n2) How could the theoretical analysis be strengthened? The rate of convergence TV norm is slow and one would expect (?!) much faster convergence in suitable norm.", "limitations": "The limitations are addressed adequately.", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "This a very well written paper which propose a new approach to constrained sampling. The approach is very well motivated and is a natural extension of SVGD. The reviewer enjoyed very much the clarity of the exposition and the conceptual clarity of the approach.\n\nA nice convergence analysis in continuous time is provided which is relatively straightforward and does not provide strong guarantees.There is lots more to be done on the theoretical side but this is a nice first step.\n\nThe model is tested on a variety of examples with satisfactory results and the model compares quite well (similar accuracy) with other algorithms on the market.", "weaknesses": "Relative weakness is that the proposed method does not seem to outperform existing algorithms.  Maybe a more detailed comparison with existing method interms of costs would be useful.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your careful review and valuable questions! We address your comments and questions as below.\n\n### Weaknesses\n\n`W1`: Relative weakness is that the proposed method does not seem to outperform existing algorithms. Maybe a more detailed comparison with existing method in terms of costs would be useful.\n\nA1: We achieved better results than MSVGD in the block experiment, and achieved better results mostly than MIED in the ring, cardioid and the double moon experiment. The latter results are in appendix D.1, which are listed below as Table 1. For large-scale experiments monotonic BNN, our method also outperforms other methods, including MIED. We list the monotonic BNN results in our main text (Table 1 in section 6.3) below as Table 2.\n\nTable 1:  Wasserstein-2 distance and energy distance between the target distribution and the variational approximation on the three toy datasets.\n\n|Name |MIED(W2 distance)| CFG(W2 distance)| MIED(Energy distance)| CFG(Energy distance) |\n|----|----|----|----|----|\n|Ring|**0.1074**| 0.1087| 0.0004|**0.0003**|\n|Cardioid|0.1240|**0.1141**| 0.0016|**0.0005**|\n|Double-moon|0.4724|**0.1660**| 0.0629|**0.0022**|\n\nTable 2: Results of monotonic Bayesian neural network under different monotonicity threshold. The results are averaged from the last 10 checkpoints for robustness.\n\n|$\\epsilon$|PD-SVGD(Test acc)|C-SVGD(Test acc)|MIED(Test acc)|CFG(Test acc)|PD-SVGD(Test NLL)|C-SVGD(Test NLL)|MIED(Test NLL)|CFG(Test NLL)|PD-SVGD(Ratio out/%)|C-SVGD(Ratio out/%)|MIED(Ratio out/%)|CFG(Ratio out/%)|\n|----|----|----|----|----|----|----|----|----|----|----|----|----|\n|$0.05$|$.647\\pm .007$|$.649\\pm .002$|$.596\\pm .002$|$\\mathbf{.661\\pm .000}$|$.634\\pm .002$|$.633\\pm .001$|$.684\\pm .000$|$\\mathbf{.632\\pm .000}$|$0.5$|$0.0$|$3.0$|$0.0$|\n|$0.01$|$.645\\pm .004$|$.650\\pm .002$|$.590\\pm .001$| $\\mathbf{.660\\pm .001}$|$.635\\pm .001$|$.634\\pm .001$|$.678\\pm .002$|$\\mathbf{.632\\pm .000}$|$0.0$|$0.0$|$5.0$|$0.0$|\n$0.005$|$.645\\pm .005$|$.650\\pm .002$|$.586\\pm .000$|$\\mathbf{.659\\pm .001}$|$.635\\pm .001$|$.633\\pm .002$|$.676\\pm .001$|$\\mathbf{.632\\pm .000}$|$0.0$|$0.0$|$6.0$|$0.0$|\n\n\n### Questions\n\n`Q1`: Can the authors accommodate multiple constraints and more complicated geometries? Are those generalizations possible? A discussion of these issues would be welcome.\n\nA1: Thanks for your question! The generalizations are possible. Please refer to the first part of the global response.\n\n`Q2`: How could the theoretical analysis be strengthened? The rate of convergence TV norm is slow and one would expect (?!) much faster convergence in suitable norm.\n\nA2: Thanks for the question! Our TV convergence rate is comparable to [1]. TV metric is broadly used in many machine learning problems because it is well defined for any two distributions. KL metric needs the initial distribution being absolute continuous with respect to the target distribution, while TV is not restricted by this requirement. Considering other suitable norm that applies to the constrained domain problem is very interesting, and we leave it for future work.\n\n[1] Zhang, R., Liu, Q., & Tong, X. T. Sampling in constrained domains with orthogonal-space variational gradient descent. Advances in Neural Information Processing Systems, 2022.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "nz822r0EXD", "reviewer": "Reviewer_VuAH", "summary": "The paper proposes to adapt particle based variational infence to the context in which the distribution of interest is supported on a subset of R^d. It is supposed that the subset on which the distribution is supported is characterized as a lower level set of a function, g, and that this function has a gradient whose norm is bounded away from zero on the complement of the supporting set. Consequently, a gradient flow which descends g when outside the set and follows gradient dynamics for which (a mollified modification of) the original distribution of interest is a fixed point within that support might be expected to provide a good approach to sampling.\n\nThe paper then concerns itself with constructing a particular (Wasserstein-type) gradient flow with such properties, establishing some theoretical properties of the proposed method and then exploring its performance empirically in a number of settings. The approach could be viewed as an extension of an existing method for dealing with distributions over manifolds in which the inequality constraint within the present work is replaced with an equality constraint and much of the novel work is concerned with addressing that difference in this particular framework.", "questions": "Time complexity seems important in comparing methods. Other than Appendix D.2 this isn't much discussed in this manuscript and that appendix isn't particularly enlightening. I think most readers would like an answer to the question \"How much computation is required to obtain an answer of a given quality (however that is assessed) using this and competing methods and how does that vary as the required quality varies?\" Are you able to give at least some heuristic numerical answer to that question?\n\nThe choice of $g$ function only seems to be made explicit for the monotonic neural network example where it is not a simple function. What were the choices used in the other numerical examples and can you offer users any guidance on how they should go about specifying the particular choice of $g$ for other problems (while it may be straightforward to specify a $g$ such that it is zero on the boundary of the domain there are clearly a great many possible choices and it isn't obvious a priori how sensitive performance will be to the choice or how to make that choice)?\n\nMore generally, how should users specify the various tuning parameters other than $g$ itself -- number of particles, architecture and parameters for the neural networks f_net and z_net, ...\n\nI don't understand the second part of Assumption 5.1: in line 139 it is stated that $p_0$ is generally assumed to be supported on all of $\\mathbb{R}^d$; in the first part of the assume the gradient of $g$ is assumed to be bounded away from zero outside the support set. How can $g$ have gradient bounded away from zero across the support of $p_0$ (less $\\Omega$) if that support is all of $R^d$ and $g$ be bounded over that support?\n\nDo you have in mind particular settings where the greater flexibility enjoyed by the CFG method in specifying constrained domains is important? I was disappointed that the numerical study didn't show such a setting. This flexibility seems like it should be a major advantage of the method and one thing that was absent from the numerical study was evidence that the proposed method either solves problems which existing methods cannot or performs substantially better in those settings. Such evidence would, I think, make the paper stronger.", "limitations": "I think this is done reasonably well.\n\nI would have liked to see more discussion of the tuning of the algorithmic parameters.", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "The problem of sampling from distributions constrained to subsets of R^d is an important one and extending the reach of modern sampling methods to this context is a valuable contribution.\n\nThe general framework is an appealing one: it allows for fairly arbitary supports to be considered, provided only that a good choice of g can be obtained to characterize the support and the general approach is intuitively appealing; if within the support of the distribution of interest then pursue a standard gradient flow approach and outside that set modify the dynamics such that the flow moves towards the support.\n\nReasonably competitive performance seems to be obtained in the numerical examples considered.", "weaknesses": "I did find the paper a little difficult to read in places and the choice/definition of notations could be improved. For example, in the definition of the set $\\Omega$ in line 97 it is not made explicit it is a subset of $\\mathbb{R}^d$ and $g$ is never identified; something like \n$$\\Omega = \\{ x \\in \\mathbb{R}^d : g(x) \\leq 0\\} \\text{ for some } g:\\mathbb{R}^d \\to \\mathbb{R}$$\nwould be much clearer. And in line 243 (as an example of an odd choice which is made throughout the manuscript to write densities evaluated at a point with that point being entirely implicit):\n$$p(\\beta) = \\mathcal{N}(0,\\sigma^2 I)$$\ndoesn't feature $\\beta$ on the RHS and looking at (17) I at first thought you intended the $p(\\beta)$ term to be uniform; the trivial edit to the more precise statement $p(\\beta) = \\mathcal{N}(\\beta; 0,\\sigma^2 I)$ would undoubtedly save many readers time. These are just some examples and while individually trivial my overall impression was that careful editing for presentation could significantly improve the manuscript.\n\nThe examples all seem rather small by modern standards with \"ground truth\" estimates obtained by rejection sampling in one case (notwithstanding the fact that scale is not the only thing which makes inference challenging and there is always scope for methods which perform well in other important contexts). Is it feasible to scale the method up, for example to larger neural networks?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your careful review and valuable questions! We address your comments and questions as below.\n\n### Weaknesses\n\n`W1`: I did find (...) improve the manuscript. \n\nA1: Thanks for your suggestion! We will modify the notations accordingly in our revision.\n\n`W2`: The examples all (...) larger neural networks?\n\nA2: Thanks for your question! Our method is feasible for learning larger neural networks and on higher dimensional real problems. Please refer to the third part of the global response. \n\n### Questions\n\n`Q1`: Time complexity seems (...) to that question?\n\nA1: Thanks for your question! Figure 8 in appendix D.2 should be able to answer your question with numerical results. Take the sample size $N=4000$ as an example. Given a quality of Wasserstein distance smaller than 0.6, the computation time of CFG is approximately 110s while MIED is approximately 320s. Changing the required quality to Wasserstein distance smaller than 0.4, the computation time of CFG change to approximately 160s while MIED change to approximately 400s. In general, CFG is more particle efficient as a functional gradient method.\n\n`Q2`: The choice of (...) make that choice)?\n\nA2: Thanks for your question! We would like to clarify that we focus on problems for which the constrained condition is provided in advance, so we choose the $g$ function to be the simplest way of representing the required constraint. For toy experiments, please refer to appendix D.1 for the analytical descriptions of the constrained domains $\\Omega$. Corresponding to $\\Omega$, we choose $g(x)=(||x||^2-1)(||x||^2-4)$ for the ring experiment, $g(x)=x_1^2+(\\frac{6}{5}x_2-x_1^{\\frac{2}{3}})^2-4$ for the cardioid experiment, $g(x)=e^{-2}-\\frac{e^{-2(x_1-3)^2}+e^{-2(x_1+3)^2}}{e^{2(||x|| - 3)^2}}$ for the double moon experiment, and $g(x)=|x_1-x_2|+|x_1+x_2|-4$ for the block experiment. For Bayesian Lasso, $\\Omega$ is the $l_q$ ball, and we choose $g(x)=||x||_q-r$ to match the problem scenario. \n\n`Q3`: More generally, (...) $f_{net}$ and $z_{net}$, ...\n\nA3: Thanks for the question! We experimented on choosing different number of particles in the Bayesian Lasso experiment, the results showed that increasing the number of particles can improve the sampling quality. Our experiment results are not very sensitive to the choice of architecture and the parameters of $f_{net}$ and $z_{net}$. Properly selecting the bandwidth $h$ according to the particle number can improve the results. One possible approach is to use the adaptive scheme to select the bandwidth as proposed in section 6.2.1 and appendix C. This scheme is backed by the error analysis in appendix C, and could be verified in the Bayesian Lasso experiment. Figure 1 in the attached pdf file showed that the adaptive scheme achieved slightly better energy distance results than the fixed bandwidth. More details on the parameter setting are provided in appendix D for each experiment. \n \n`Q4`: I don't understand (...) over that support?\n\nA4: Thanks for your question! This assumption is adopted from [1], which is to guarantee that the particles outside the domain can enter the constrained domain. In practice, we could take a sufficiently large domain that covers all the initial particles and the given constrained domain. $g$ could be bounded on this domain and $\\nabla g$ is bounded away from zero, providing feasibility of training. This assumption does not serve for our main focus, which is the convergence analysis when the particles are inside the constrained domain. \n  \n[1] Zhang, R., Liu, Q., & Tong, X. T. Sampling in constrained domains with orthogonal-space variational gradient descent. Advances in Neural Information Processing Systems, 2022.\n\n`Q5`: Do you have (...) the paper stronger.\n \nA5: Thanks for your question! Firstly, from the best of our knowledge, there are not many methods that can directly deal with non-convex and complex hard constrained domains besides our method and the MIED method. As stated in the introduction, the mirror method is restricted to problems where mirror maps can be obtained, while projection methods suffers from high computation cost for complex domains because of the complexity in finding projection points. Additionally, compared to methods that can deal with this setting such as MIED, our method achieved better results on the cardioid, the double moon and the monotonic BNN experiment. We list Table 6 in appendix D.1 as below for the comparisons on the toy experiments. Secondly, using neural network to learn the optimal velocity field is more powerful way to implement ParVI, as it enhances the expressiveness and flexibility compared to classical kernel methods ([1][2][3]). This advantage can be verified by our better results compared to kernel based methods MSVGD, PD-SVGD and C-SVGD. Using the simple and flexible neural network framework, our method is also more scalable than MIED in terms of particle numbers N. Please refer to appendix D.2 for detailed discussion and verification.\n\nTable:  Wasserstein-2 distance and energy distance between the target distribution and the variational approximation on the three toy datasets.\n|Name|MIED(W2 distance)|CFG(W2 distance)|MIED(Energy distance)|CFG(Energy distance)|\n|----|----|----|----|----|\n|Ring|**0.1074**|0.1087|0.0004|**0.0003**|\n|Cardioid|0.1240|**0.1141**|0.0016|**0.0005**|\n|Double-moon|0.4724|**0.1660**|0.0629|**0.0022**|\n\n[1] di Langosco, L. L., Fortuin, V., & Strathmann, H. (2021). Neural variational gradient descent. arXiv preprint arXiv:2107.10731.\n\n[2] Dong, H., Wang, X., Lin, Y., & Zhang, T. (2022). Particle-based variational inference with preconditioned functional gradient flow. arXiv preprint arXiv:2211.13954.\n\n[3] Cheng, Z., Zhang, S., Yu, L., & Zhang, C. (2024). Particle-based Variational Inference with Generalized Wasserstein Gradient Flow. Advances in Neural Information Processing Systems, 36.\n\n### Limitations\n\nI would (...) parameters.\n\nA1: Please refer to A3 of Q3 of our response for more details.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
We thank all reviewers for their constructive feedback, and will modify our paper accordingly in our revision. Here we address some of the common issues raised by the reviewers.

**Geometric generalization of our method**

Some of the reviewers are interested whether our proposed method can generalize and accommodate multiple constraints (including more equality and inequality constraints) and more complicated geometries. For additional equality constraints, our proposed framework can still apply by adopting the idea of [1]. We could use velocity fields like $v_{\sharp}(x)$ to guided the particles to satisfy the equality constraints, and use our method to propose velocity fields substituting $v_{\perp}(x)$ to satisfy the inequality constraints. Adding these two types of velocity presents the desired velocity field. To illustrate the idea, we demonstrate the training process of a 3D toy example in Figure 2 of the attached pdf file. Suppose the coordinate of the particle is $(x,y,z)$. The target distribution is a truncated standard gaussian located in the ring shaped domain in xOy plane. This corresponds to the equality constraint $z=0$ and the inequality constraint $1\le x^2+y^2\le 4$. The initial distribution is the 3D standard gaussian distribution. We can see that the particles collapse to the xOy plane and converge inside the ring domain. For additional inequality constraints, we may need multiple velocity fields outside the constrained domain for particles to enter the constrained domain, and the boundary integral term can be similarly estimated using band-wise approximation.

[1] Zhang, R., Liu, Q., & Tong, X. T. Sampling in constrained domains with orthogonal-space variational gradient descent. Advances in Neural Information Processing Systems, 2022.


**The assumptions and parameters on the domain entering phase**

Some of the reviewers are interested in the effect of lambda in line 158, and some are interested in assumption 5.1. These questions are about the phase when particles enter the constrained domain. We would like to be clear that the main focus of our problem is on the convergence when the particles are inside the constrained domain. The feasibility of entering the constrained domain is the basic requirement of the training. The assumptions on entering the constrained domain could be relaxed as long as sufficient enough that the particles will not stuck at local modes outside the domain. Lambda in line 158 controls the entering velocity of the particles. With respect to the problem scale and dimensionality, we select lambda in a reasonable range to achieve swift entering without overshooting the constrained domain.


**Scaling our method up to real-world application**

Some of the reviewers are interested whether our method could scale up to high dimensional real problems. We experimented on the COMPAS dataset using larger monotonic BNNs, the number of neurons in the layer of BNN increased from 50 to 100, making the particle dimension increase from 903 to 1502. From Table 1 in the attached pdf file, our method still achieved the best accuracy results compared to other methods, while achieving competitive log-likelihood results. To test our method on higher dimension, we additionally experimented on the 276-dimensional larger dataset Blog Feedback ([1]) using monotonic BNN. The particle dimension is 13903. From Table 2 in the attached pdf file, our method achieved the best result.

[1] Liu, X., Han, X., Zhang, N., & Liu, Q. (2020). Certified monotonic neural networks. Advances in Neural Information Processing Systems, 33, 15427-15438.


We hope our response has adequately addressed the reviewers' questions and concerns, and look forward to reading any other additional comments.

Meta-review: This paper investigates the constrained sampling problem in Particle-based Variational Inference methods and proposes a framework called Constrained Functional Gradient Flow (CFG) to restrict samples within a specific domain. This is achieved by introducing a boundary condition for the gradient flow. The paper also provides a convergence guarantee in total variation for the proposed method.

Almost all reviewers agree that the paper merits acceptance, and after carefully reading the rebuttal and discussion, I tend to concur. Please incorporate the reviewers’ feedback into the final version of the paper, with particular attention to the following concerns:

1- Include the geometric generalization of your method and clarify the assumptions about entering the boundary and how to set $\lambda$ and other hyper-parameters. (Response to Rev s6Xw, and YhS1)

2- Include the experimental results of real-world applications. 

3- Include the related baseline and references and their complexity comparison. (Response to Rev YhS1, PQmA, and VuAH)
Decision: Accept (poster)