Title: SURROGATE-BASED CONSTRAINED LANGEVIN SAMPLING WITH APPLICATIONS TO OPTIMAL MATERIAL CONFIGURATION DESIGN

Abstract: We consider the problem of generating configurations that satisfy physical constraints for optimal material nano-pattern design, where multiple (and often conflicting) properties need to be simultaneously satisfied. Consider, for example, the trade-off between thermal resistance, electrical conductivity, and mechanical stability needed to design a nano-porous template with optimal thermoelectric efficiency. To that end, we leverage the posterior regularization framework and show that this constraint satisfaction problem can be formulated as sampling from a Gibbs distribution. The main challenges come from the black-box nature of those physical constraints, since they are obtained via solving highly non-linear PDEs. To overcome those difficulties, we introduce Surrogate-based Constrained Langevin dynamics for black-box sampling. We explore two surrogate approaches. The first approach exploits zero-order approximation of gradients in the Langevin Sampling and we refer to it as Zero-Order Langevin. In practice, this approach can be prohibitive since we still need to often query the expensive PDE solvers. The second approach approximates the gradients in the Langevin dynamics with deep neural networks, allowing us an efficient sampling strategy using the surrogate model. We prove the convergence of those two approaches when the target distribution is log-concave and smooth. We show the effectiveness of both approaches in designing optimal nano-porous material configurations, where the goal is to produce nano-pattern templates with low thermal conductivity and reasonable mechanical stability.

1. INTRODUCTION:

In many real-world design problems, the optimal design needs to simultaneously satisfy multiple constraints, which can be expensive to estimate. For example, in computational material design, the goal is to come up with material configurations, or samples, satisfying a list of physical constraints that are given by black-box numerical Partial Differential Equations (PDE) solvers. Such solvers (for example, the Boltzmann Transport Equation solver) are often complex, expensive to evaluate, and offer no access to their inner variables or their gradients.

We pose this design-under-constraints problem as sampling from a Gibbs distribution defined on some compact support. The problem of sampling from a distribution with unknown likelihood that can only be point-wise evaluated is called black-box sampling (Chen & Schmeiser, 1998; Neal, 2003) . We show in this paper that constrained black-box sampling can be cast as a constrained Langevin dynamics with gradient-free methods. Zero-order optimization via Gaussian smoothing was introduced in Nesterov & Spokoiny (2017) and extended to black-box sampling with Langevin dynamics in Shen et al. (2019) . We extend this approach to the constrained setting from a black-box density with compact support.

However, one shortcoming of this approach is that it is computationally very expensive since it requires repeatedly querying PDE solvers in order to get an estimate of the gradient. To alleviate computational issues, we propose Surrogate Model Based Langevin dynamics, that consists of two steps: (i) Learning (using training data) an approximation of the gradient of the potential of the Gibbs distribution. We show that learning the gradient, rather than the potential itself, is important for the mixing of the Langevin dynamics towards the target Gibbs distribution. We devise several objective functions, as well as deep neural-network architectures for parameterizing the approximating function class, for learning the gradient of the potential function. (ii) We then use the surrogate gradient model in the constrained Langevin dynamics in lieu of the black-box potential. Using the surrogate enables more efficient sampling, since it avoids querying the expensive PDE solvers, and obtaining gradients is as efficient as evaluating the functions themselves using automatic differentiation frameworks such as PyTorch or TensorFlow.

To summarize, our main contributions are as follows:

1. We cast the problem of generating samples under constraints in the black-box setting as sampling from a Gibbs distribution. 2. We introduce Constrained Zero-Order Langevin Monte Carlo, using projection or proximal methods, and provide the proof of its convergence to the target Gibbs distribution. 3. We introduce Surrogate Model Based Projected Langevin Monte Carlo via learning the gradient of the potential of the Gibbs distribution using deep neural networks or reproducing kernel spaces, and prove its convergence to the target distribution when used in conjunction with projection or proximal based methods. We shed the light on the importance of the approximation of the gradient of the potential, and we show how to achieve this using Hermite and Taylor learning. 4. We showcase the usability and effectiveness of the proposed methods for the design of nanoporous configurations with improved thermoelectric efficiency. The design consists of finding new configurations with optimized pore locations, such that the resulting configurations have favorable thermal conductivity (i.e., minimal κ) and desired mechanical stability (von Mises Stress σ ≤ τ , where τ is some preset threshold).1 

2. FROM CONSTRAINTS SATISFACTION TO SAMPLING FROM A GIBBS DISTRIBUTION: POSTERIOR REGULARIZATION:

In black-box optimization problems (such as the material design under consideration), the goal is to find a posterior distribution q of samples satisfying a list of equality and inequality constraints: ψ j (x) = y k , j = 1 . . . C e , and φ k (x) ≤ b k , k = 1 . . . C i where x ∈ Ω and Ω ⊂ R d is a bounded domain. We assume a prior distribution p 0 (whose analytical form is known). The main challenge in black-box optimization is that the functions ψ j and φ k can be only evaluated point-wise, and neither do we have functional forms nor access to their gradients. For example, ψ and φ might be obtained via aggregating some statistics on the solution of a nonlinear PDE given by a complex solver.

To make the problem of learning under constraints tractable, we choose Lagrangian parameters λ j > 0 and obtain the following relaxed objective:

min q, Ω q(x)=1 KL(q, p 0 ) + Ce j=1 λ j E x∼q (ψ j (x) -y k ) 2 + Ci k=1 λ k E x∼q (φ k (x) -b k ) + (1)

The formulation in Eq. 1 is similar in spirit to the posterior regularization framework of Ganchev et al. (2010) ; Hu et al. (2018) . However, we highlight two differences: (i) our focus is on constrained settings (where Ω is bounded), and (ii) we assume a black-box setting. We first obtain: Lemma 1 (Constraint Satisfaction as Sampling from a Gibbs Distribution). The solution to the distribution learning problem given in Eq. 1 is given by:

π(x) = exp(-U (x)) Z 1 x∈Ω (2)

where

U (x) = -log p 0 (x) + Ce j=1 λ j (ψ j (x) -y k ) 2 + Ci k=1 λ k (φ k (x) -b k ) + and Z = x∈Ω exp (-U (x)) dx.

Lemma 1 shows that the constraint satisfaction problem formulated in Eq. 1 amounts to sampling from a Gibbs distribution defined on a compact support given in Eq. 2. Sampling from a Gibbs distribution (also known as Boltzmann distribution) has a long history using Langevin dynamics. In the white-box setting when the functions defining the constraints have explicit analytical forms as well as their gradients, Langevin dynamics for Gibbs distribution sampling defined on a compact domain Ω and their mixing properties were actively studied in Bubeck et al. (2015) ; Brosse et al. (2017) . In the next Section, we provide a more detailed review.

Remark 1 (Relation to Bayesian Optimization). While in Bayesian optimization we are interested in finding a point that satisfies the constraints, in our setting we are interested in finding a distribution of candidate samples that satisfy (black-box) constraints. See (Suzuki et al., 2019) for more details.

Remark 2. For the rest of the paper, we will assume p 0 to be the uniform distribution on Ω, which means that its gradients are zero on the support of the domain Ω. Otherwise, if p 0 is known and belongs to, for instance, an exponential family or a generative model prior (such as normalizing flows), we can sample from π using a mixture of black-box sampling on the constraints (ψ j , φ k ) and white-box sampling on log(p 0 ).

3. WHITE-BOX SAMPLING: CONSTRAINED LANGEVIN DYNAMICS:

We review in this section Langevin dynamics in the unconstrained case (Ω = R d ) and the constrained setting (Ω ⊂ R d ). Below, • denotes the Euclidean norm unless otherwise specified. We are interested in sampling from

EQUATION

Preliminaries. We give here assumptions, definitions and few preliminary known facts that will be useful later. Those assumptions are commonly used in Langevin sampling analysis (Dalalyan, 2017; Bubeck et al., 2015; Brosse et al., 2017; Durmus et al., 2019) .

1. Assumption A: We assume Ω is a convex such that 0 ∈ Ω, Ω contains a Euclidean ball of radius r, and Ω is contained in a Euclidean ball of radius R. (For example, Ω might encode box constraints.)

The projection onto Ω, P Ω (x) is defined as follows: for all x ∈ Ω,

P Ω (x) = arg min z∈Ω x -z 2 . Let R = sup x,x ∈Ω ||x -x || < ∞.

2.. Assumption B::

We assume that U is convex, β-smooth, and with bounded gradients:

∇ x U (x) -∇ y U (y) ≤ β x -y , ∀x, y ∈ Ω (β-smoothness). ∇U (x) ≤ L, ∀x ∈ Ω (Boundedness).

The Total Variation (TV) distance between two measures µ, ν is defined as follows:

T V (µ, ν) = sup A |µ(A) -ν(A)|. Pinsker Inequality relates KL divergence to TV: T V (µ, ν) ≤ 2KL(µ, ν).

Unconstrained Langevin Dynamics. In the unconstrained case, the goal is to sample from a Gibbs distribution π(x) = exp(-U (x))/Z that has unbounded support. This sampling can be done via the Langevin Monte Carlo (LMC) algorithm, which is given by the following iteration:

EQUATION

where ξ k ∼ N (0, I d ), η is the learning rate, and λ > 0 is a variance term.

Constrained Langevin Dynamics. In the constrained case, the goal is to sample from π(x) = exp(-U (x))/Z1 x∈Ω ,. We discuss two variants:

Projected Langevin Dynamics. Similar to projected gradient descent, Bubeck et al. (2015) introduced Projected Langevin Monte Carlo (PLMC) and proved its mixing propreties towards the stationary distribution π. PLMC is given by the following iteration :

EQUATION

In essence, PLMC consists of a single iteration of LMC, followed by a projection on the set Ω using the operator P Ω .

Proximal Langevin Dynamics. Similar to proximal methods in constrained optimization, Brosse et al. (2017) introduced Proximal LMC (ProxLMC) that uses the iteration:

X k+1 = 1 - η γ X k -η∇ x U (X k ) + η γ P Ω (X k ) + 2ληξ k , k = 0 . . . K -1, (ProxLMC) (6)

where η is the step size and γ is a regularization parameter. In essence, ProxLMC (Brosse et al., 2017) performs an ordinary LMC on U γ (x) = U (x) + i γ Ω (x), where i γ Ω (x) is the proximal operator:

i γ Ω (x) = inf y i Ω (x) + (2γ) -1 x -y 2 = (2γ) -1 x -P Ω (x) 2 ,

where i Ω (x) = 0 for x ∈ Ω and i Ω (x) = ∞ for x / ∈ Ω. Therefore, the update in Eq. 6 is a regular Langevin update (as in Eq. 4) with potential gradient

∇ x U γ (x) = ∇ x U (x) + γ -1 (x -P Ω (x)).

We denote by µ PLMC K and µ ProxLMC K the distributions of X K obtained by iterating Eq. 5 and Eq. 6 respectively. Under Assumptions A and B, both these distributions converge to the target Gibbs distribution π in the total variation distance. In particular, Bubeck et al. (2015) showed that for η = Θ(R 2 /K), we obtain:

EQUATION

Likewise, Brosse et al. (2017) showed that for 0 < η ≤ γ(1 + β 2 γ 2 ) -1 , we obtain:

EQUATION

where the notation

α n = Ω(β n ) means that there exists c ∈ R, C > 0 such that α n ≥ Cβ n log c (β n ).

4. CONSTRAINED LANGEVIN DYNAMICS IN THE BLACK-BOX SETTING:

We now introduce our variants of constrained LMC for the black-box setting where explicit potential gradients are unavailable. We explore in this paper two strategies for approximating the gradient of U in the black-box setting. In the first strategy, we borrow ideas from derivative-free optimization (in particular, evolutionary search). In the second strategy we learn a surrogate deep model that approximates the gradient of the potential. Below, let G : Ω → R d be a vector valued function that approximates the gradient of the potential, ∇ x U . We make:

Assumption C. The surrogate gradient G satisfies E G(Y k ) 2 < ∞, ∀k.

Surrogate Projected Langevin Dynamics. Given Y 0 , the Surrogate Projected LMC (S-PLMC) replaces the potential gradient ∇ x U in Eq. 5 with the surrogate gradient G:

EQUATION

Surrogate Proximal Langevin Dynamics. Similarly, the Surrogate Proximal LMC (S-ProxLMC) replaces the unknown potential gradient ∇ x U in Eq. 6 with the gradient surrogate G:

Y k+1 = 1 - η γ Y k -ηG(Y k ) + η γ P Ω (Y k ) + 2ληξ k , k = 0 . . . K -1 (S-ProxLMC) (10)

We now present our main theorems on the approximation properties of surrogate LMC (S-PLMC, and S-ProxLMC). We do so by bounding the total variation distance between the trajectories of the surrogate Langevin dynamics (S-PLMC, and S-ProxLMC) and the true LMC dynamics (PLMC and ProxLMC). Theorem 1 is an application of techniques in Stochastic Differential Equations (SDE) introduced in Dalalyan & Tsybakov (2012) and is mainly based on a variant of Grisanov's Theorem for change of measures (Lipster & Shiryaev, 2001 ) and Pinsker's Inequality that bounds total variation in terms of Kullback-Leibler divergence.

Theorem 1 (S-PLMC and S-ProxLMC Mixing Properties). Under Assumption C, we have:

1. S-PLMC Convergence. Let µ PLMC K be the distribution of the random variable X K obtained by iterating PLMC Eq. 5, and µ S-PLMC K be the distribution of the random variable Y K obtained by iteration S-PLMC given in Eq. 9. We have:

EQUATION

)

2. S-ProxLMC Convergence. Let µ ProxLMC K be the distribution of the random variable X K obtained by iterating ProxLMC Eq. 6, and µ S-ProxLMC K be the distribution of the random variable Y K obtained by iterating S-ProxLMC given in Eq. 10. We have:

EQUATION

)

From Theorem 1, we see that it suffices to approximate the potential gradient ∇ x U (X) (and not the potential U (X)) in order to guarantee convergence of surrogate-based Langevin sampling. Using the triangle inequality, and combining Theorem 1 and bounds in Eqs 7 and 8 we obtain: Theorem 2. (Convergence of Surrogate Constrained LMC to the Gibbs distribution.) Under assumptions A,B and C we have:

1. Assume in S-PLMC that there exists

δ > 0 such that E G(Y k ) -∇ x U (Y k ) 2 ≤ δ, ∀k ≥ 0. Set λ = 1, and η = Θ(min(R 2 /K, α/K 2 )) where α = 1/(δ + β 2 R 2 ) . Then for K = Ω(ε -12 d 12 ),

we have:

T V (µ S-PLMC K , π) ≤ ε. 2. Assume in S-ProxLMC that there exists δ > 0 such that E G(X k ) -∇ x U (X k ) 2 ≤ δ, ∀k ≥ 0.

Set λ = 1, and η = min(γ(1 + β 2 γ 2 ) -1 , 1 δK 2 ). Then for K = Ω(ε -6 d 5 ) we have:

T V (µ S-ProxLMC K , π) ≤ ε.

5. ZERO-ORDER CONSTRAINED LANGEVIN DYNAMICS:

In zero-order optimization (Nesterov & Spokoiny, 2017; Duchi et al., 2015; Ghadimi & Lan, 2013; Shen et al., 2019) , one considers the Gaussian smoothed potential U ν defined as U ν (x) = E g∼N (0,I d ) U (x + νg), and its gradient is given by ∇

EQUATION

where g 1 , . . . g n are i.i.d. standard normal vectors.

Zero-Order sampling from log-concave densities was recently studied in Shen et al. (2019) . We extend it here to the constrained sampling case of log-concave densities with compact support. We define Constrained Zero-Order Projected LMC (Z-PLMC) and Zero-Order Proximal LMC (Z-ProxLMC) by setting G(x) = Ĝn U(x) in Eq. 9 and Eq. 10 respectively. Lemma 2 (Zero-Order Gradient Approximation (Nesterov & Spokoiny, 2017; Shen et al., 2019) ).

Under Assumption B, we have for all x ∈ Ω:

EQUATION

Thanks to Lemma 2 that ensures uniform approximation of gradients in expectation, we can apply Theorem 2 and get the following corollary for Z-PLMC and Z-ProxLMC: Corollary 1 (Zero-order Constrained Langevin approximates the Gibbs distribution). Under As-

sumptions A and B, let δ ∈ [0, 1], for n ≥ βν(d + 2) 3/2 + (d + 1) 1 2 L 2

/δ, we have the following bounds in expectation:

1. Set λ = 1, and η = Θ(min(R 2 /K, α/K 2 )) where α = 1/(δ + β 2 R 2 ) . For K = Ω(ε -12 d 12 ),

we have:

EQUATION

)

Remark 3. For simplicity, we state the above bound in terms of expectations over the randomness in estimating the gradients. It is possible to get finite-sample bounds using the Vector Bernstein concentration inequality, coupled with covering number estimates of Ω but omit them due to space.

6. SURROGATE MODEL BASED CONSTRAINED LANGEVIN DYNAMICS:

Despite its theoretical guarantees, zero-order constrained Langevin (Z-PLMC and Z-ProxLMC) has a prohibitive computation cost as it needs O(nK) black-box queries (in our case, invocations of a nonlinear PDE solver). To alleviate this issue, we introduce in this Section a neural surrogate model as an alternative to the gradient of the true potential.

6.1. HERMITE LEARNING OF GRADIENTS: JACOBIAN MATCHING OF ZERO-ODER ESTIMATES:

From Theorem 2, we saw that in order to guarantee the convergence of constrained Langevin dynamics, we need a good estimate of the gradient of the potential of the Gibbs distribution. Recall that the potential given in Lemma 1 depends on ψ j and φ k , which are scalar outputs of computationally heavy PDE solvers in our material design problem. To avoid this, we propose to train surrogate neural network models approximating each PDE output and their gradients. Concretely, suppose we are given a training set S for a PDE solver for the property ψ (dropping the index j for simplicity):

S = {(x i , y i = ψ(x i ), ỹi = Ĝn ψ(x i )), x i ∼ ρ Ω i.i.d., i = 1, . . . , N },

where ρ Ω is the training distribution and Ĝn ψ(.) is the zero-order estimate of the gradient of ψ given in Eq. 13. We propose to learn a surrogate model belonging to a function class H θ , fθ ∈ H θ , that regresses the value of ψ and matches the zero-order gradient estimates as follows:

EQUATION

The problem in Eq. 17 was introduced and analyzed in Shi et al. (2010) where H θ is a ball in a Reproducing Kernel Hilbert Space (RKHS). Following Shi et al. (2010) , we refer to this type of learning as Hermite Learning. In the deep learning community, this type of learning is called Jacobian matching and was introduced in Srinivas & Fleuret (2018); Czarnecki et al. (2017) where H θ is a deep neural network parameterized with weights θ. When f θ is a deep network, we can optimize this objective efficiently using common deep learning frameworks (PyTorch, TensorFlow). (Shi et al., 2010) have shown that when H θ is an RKHS ball and when ỹi = ∇ x ψ(x i ) are exact gradients, for a sufficiently large training set with N = O(1/ 1/(2rζ) ) (where r, ζ are exponents in [0, 1] that depend on the regularity of the function ψ). Under the assumption that ψ ∈ H θ we have:

Ω ∇ x f θ (x) -∇ x ψ(x) 2 ρ Ω (x)dx ≤ .

Since we are using inexact zero-order gradients, we will incur an additional numerical error that is also bounded as shown in Lemma 2.

6.2. TAYLOR LEARNING OF GRADIENTS:

While Jacobian matching of zero-order gradients is a sound approach, it remains expensive to construct the dataset, as we need for each point to have 2n + 1 queries of the PDE solver. We exploit in this section the Taylor learning framework of gradients that was introduced in Mukherjee & Zhou (2006) ; Mukherjee & Wu (2006), and Wu et al. (2010) . In a nutshell, Mukherjee & Zhou (2006) suggests to learn a surrogate potential f θ and gradient G Λ that are consistent with the first-order taylor expansion. Given a training set 2010) suggest the following objective:

S = {(x i , y i = ψ(x i )), x ∼ ρ Ω , i = 1 . . . N }, Wu et al. (

EQUATION

where

w σ ij = exp -xi-xj 2 σ 2

, H θ is an RKHS ball of scalar valued functions, and H d Λ is an RKHS ball of vector valued functions.

Under mild assumptions, Mukherjee & Zhou (2006) shows that we have for

N = O(1/ d/2 ): Ω G Λ (x) -∇ x ψ(x)

2 ρ Ω (x)dx ≤ . We simplify the problem in Eq. 18 and propose the following two objective functions and leverage the deep learning toolkit to parameterize the surrogate f θ :

EQUATION

min f θ ∈H θ 1 N N i=1 {(y i -f θ (x i )) 2 + λ N 2 i,j w σ ij (y i -y j + ∇ x f θ (x i ), x j -x i ) 2 }, (Taylor-Reg) . (20)

The objective in Eq. 19 uses a single surrogate to parameterize the potential and its gradient. The objective in Eq. 20 is similar in spirit to the Jacobian matching formulation in the sense that it adds a regularizer on the gradient of the surrogate to be consistent with the first-order Taylor expansion in local neighborhoods. The advantage of the Taylor learning approach is that we do not need to perform zero-order estimation of gradients to construct the training set and we rely instead on first-order approximation in local neighborhood.

6.3. SURROGATE MODEL CONSTRAINED LMC:

Consider the surrogate model f θ obtained via Hermite Learning (Eq. 17) or via Taylor learning (Eqs 18, 19, 20) . We are now ready to define the surrogate model LMC by replacing

G(x) = ∇ x f θ (x)

in the constrained Langevin dynamics in Eqs 9 and 10.

Both Hermite and Taylor learning come with theoretical guarantees when the approximation function space is an RKHS under some mild assumptions on the training distribution and the regularity of the target function ψ. In Hermite learning (Theorem 2 in Shi et al. ( 2010)) we have: 2ζr) ) (where exponents ζ, r ∈ [0, 1] depend on regularity of ψ). In Taylor Learning with the objective function given in Eq. 18 (Proposition 7 in Wu et al. ( 2010) we have:

E x∼pΩ ∇ x f θ (x) -∇ x ψ(x) 2 ≤ for sufficiently large training set N = O(1/ 1/(

E x∼ρΩ G Λ (x) -∇ x ψ(x) 2 ≤ for N = O(1/ d/2

). In order to apply Theorem 2 we need this gradient approximation error to hold in expectation on all intermediate distributions in the Langevin sampling. Hence, we need the following extra-assumption on the training distribution p Ω :

Assumption D: Assume we have a learned surrogate G on training distribution ρ Ω such that E x∼ρΩ G(x) -∇ x U (x)

2 ≤ . Assume ρ Ω (x) > 0, ∀x ∈ Ω and that it is a dominating measure of Langevin (PLMC, S-PLMC, Prox-LMC, S-ProxLMC ) intermediate distributions µ k , i.e. there exists C > 0 such that:

µ k (x) ≤ Cρ Ω (x), ∀x ∈ Ω, ∀k = 0, . . . K -1. Under Assumption D, it follows immediately that E G(X k ) -∇U (X k ) 2 = Ω G(x) -∇U (x 2 µ k (x) ρ Ω (x) ρ Ω (x) ≤ C

and hence we can apply Theorem 2 for δ = C , and we obtain ε-approximation of the target Gibbs distribution in terms of total variation distance. In what follows we refer to surrogate constrained LMC, as x-PLMC or x-ProxLMC where x is one of four suffixes ({Z-Hermite, Taylor-2, Taylor-1, Taylor-Reg}).

7. RELATED WORK:

Zero-Order Methods. Zero-order optimization with Gaussian smoothing was studied in Nesterov & Spokoiny (2017) and Duchi et al. (2015) in the convex setting. Non-convex zero order optimization was also addressed in Ghadimi & Lan (2013) . The closest to our work is the zero-order Langevin Shen et al. (2019) introduced recently for black-box sampling from log concave density. The main difference in our setting, is that the density has a compact support and hence the need to appeal to projected LMC (Bubeck et al., 2015) and Proximal LMC (Brosse et al., 2017) . It is worth nothing that Hsieh et al. (2018) introduced recently mirror Langevin sampling that can also be leveraged in our framework.

Gradients and Score functions Estimators. We used the approach of gradient distillation (Srinivas & Fleuret, 2018) and learning gradients of (Wu et al., 2010) , since they are convenient for training on different constraints and they come with theoretical guarantees. However, other approaches can be also leveraged such as the score matching approach for learning the gradient of the log likelihood (Hyvärinen, 2005) and other variants appealing to dual embeddings (Dai et al., 2018) . Estimating gradients can be also performed using Stein's method as in (Li & Turner, 2017) , or via maintaining a surrogate of the gradient as in Stein descent without gradient (Han & Liu, 2018) .

Optimization approaches. Due to space limitation, we restrict the discussion to the optimization methods that are most commonly and recently used for optimal material (or molecule) design. A popular approach to deal with optimization of expensive black-box functions is Bayesian Optimization (BO) (Mockus, 1994; Jones et al., 1998; Frazier, 2018) The advantage of our framework is that we obtain a distribution of optimal configurations (as opposed to a single optimized sample) that does not rely on training on a specific pre-existing dataset and can be further screened and tested for their optimality for the task at hand.

8. EXPERIMENTS:

In this section, we demonstrate the usability of our black-blox Langevin sampling approach for the design of nano-porous configurations. We first show the performance of the surrogate models in learning the potential function, showcasing the results using four different variants: standard regression, Taylor regularization, Taylor-1 and Taylor-2. We then show how well the surrogate-based Langevin MC generates new samples under the thermal and mechanical constraints. We compare the sample quality on multiple criteria between the surrogate and zero-order approaches with either projection or proximal update step.

Data. We want to learn surrogate models to approximate the gradient of the potential from data.

To this end, we generate a dataset of 50K nano-porous structures, each of size 100nm × 100nm.

One such example is displayed in Fig. 1 . Number of pores is fixed to 10 in this study and each pore is a square with a side length of 17.32nm. We sample the pore centers uniformly over the unit square and construct the corresponding structure after re-scaling them appropriately. Then, using the solvers OpenBTE (Romano & Grossman, 2015) and Summit ( MIT Development Group, 2018), we obtain for each structure x a pair of values: thermal conductivity κ and von Mises stress σ. Finally, we collect two datasets: {(x i , κ i )} N i=1 and {(x i , σ i )} N i=1 with the same inputs x i 's and N = 50K samples. More details are given in Appendices B and C on the PDEs and their corresponding solvers.

Features. The pore locations are the natural input features to the surrogate models. Apart from the coordinates, we also derive some other features based on physical intuitions. For example, the distances between pores and the alignment along axes are informative of thermal conductivity (Romano & Grossman, 2016) . As such, we compute pore-pore distances along each coordinate axis and add them as additional features.

Surrogate gradient methods. We use feed-forward neural networks to model the surrogates since obtaining gradients for such networks is efficient thanks to automatic differentiation frameworks. We use networks comprised of 4 hidden layers with sizes 128, 72, 64, 32 and apply the same architecture to approximate the gradients for κ and σ separately. The hidden layers use ReLU activations whereas sigmoid was used at the output layer (after the target output is properly normalized). For the Taylor-2 variant (in Eq. 18), we have an additional output vector of the same size as the input for the gradient prediction. The networks are trained on the corresponding objective functions set up earlier by an Adam optimizer with learning rate 10 -4 and decay 1.0. We fine-tune the networks with simple grid-search and select the best models for comparison. Due to the space constraint, we present the results in Appendix A and emphasize that Z-Hermite is not included in the entire comparison but in a small experiment performed with a more lightweight OpenBTE version.

Incorporating constraints and comparison metrics. We demonstrate the usability of our proposed black-box Langevin sampling for the design of nano-configurations under thermal conductivity and mechanical stability constraints that are provided by the corresponding PDE solvers. To compare sampling outcomes, we use the following metrics. We report the minimum value of κ and Monte Carlo estimates for both κ and σ to compare the samples generated by different sampling methods and surrogate models. The Monte Carlo estimates are computed on 20 samples.

Single constraint. Our first task is to design nano-configurations under the thermal conductivity constraint where we want κ as low as possible in order to achieve high thermo-electric efficiency.

From the posterior regularization formulation Section 2, we pose the constraint satisfaction as sampling from the following Gibbs distribution:

EQUATION

where p 0 (x) is the uniform distribution over the unit square, which is equivalent to the Poisson process of 10 pores on the square, and κ(x) is the thermal conductivity we want to minimize. Starting from 20 samples initialized from p 0 (x), we run our proposed black-box Langevin MCs and obtain 20 new realizations from the target distribution π(x). We use four different surrogates (including simple regression, Taylor-Reg, Taylor-1 and zero-order) and each surrogate with either projection or proximal update. We show the summary statistics of these samples in Table 1 . The regression-PMLC in the first row and regression-ProxLMC in the fifth represent the sampling where the surrogate model are fitted on solely the mean square error objective. In all methods, we set λ = 100, the step size η = 1e-3 and the exponential decay rate 0.8. Since keeping track of the true κ value is expensive, we stop after K = 10 iterations. We first observe that the regression-based method (PLMC, ProxLMC) is less effective than the others simply because they do not have an implicit objective for approximating the gradients. Taylor-Reg and Taylor-1 demonstrate its effectiveness in approximating the gradient and are able to achieve lower thermal conductivity. In particular, Taylor-1-ProxLMC and Zero-order-PLMC perform in the similar range in terms of the minimum achieved, but the learned surrogate offers 17x speed up (per sample) over zero order methods. Due to the space limit, we do not report Taylor-2 results in Multiple constraints. Achieving the minimal thermal conductivity can be fulfilled without much difficulty (e.g. structures with all pores aligned along the vertical axis), but such structures are often mechanically unstable. In the next step, we study whether adding more (conflicting) constraints helps us design better nano-configurations. Hence, we consider both thermal conductivity κ and mechanical stability provided via von Mises stress σ. We want a sample x that minimizes κ(x) to achieve high thermo-electric efficiency while maintaining σ(x) less than some threshold (which we explain below). Like the single constraint case, we pose this as sampling from the following Gibbs distribution:

EQUATION

where p 0 (x) is the same as above, σ(x) is the von Mises stress and τ is a threshold on the maximum value of σ. With this framework, we relax the inequality constraint to the Hinge loss term on von Mises stress. The results are summarized in Table 2 . Note that all the surrogate Langevin MCs are initialized from the same set of 20 samples as above. In this experiment, we set τ = 0.5, λ 1 = 100, λ 2 = 10 the step size η = 1e-3 and the exponential decay rate 0.8. Comparing with Table 1 , one can see that not only better κ be achieved but also the σ can be reduced simultaneously. These results suggest that our approach can effectively sample new configurations under multiple competing constraints. Examples of new nano-configurations are show in Fig. 1 and Appendix A Fig. 

9. CONCLUSION:

In this paper we introduced Surrogate-Based Constrained Langevin Sampling for black-box sampling from a Gibbs distribution defined on a compact support. We studied two approaches for defining the surrogate: the first through zero-order methods and the second via learning gradient approximations using deep neural networks. We showed the proofs of convergence of the two approaches in the log-concave and smooth case. While zero-order Langevin had prohibitive computational cost, learned surrogate model Langevin enjoy a good tradeoff of lightweight computation and approximation power. We applied our black-box sampling scheme to the problem of nano-material configuration design, where the black box constraints are given by expensive PDE solvers, and showed the efficiency and the promise of our method in finding optimal configurations. Among different approaches for approximating the gradient, the zero-order ones (PLMC, ProxLMC) show overall superior performance, at a prohibitive computational cost. We established that the deep the surrogate (Taylor-1 ProxLMC) is a viable alternative to zero-order methods, achieving reasonable performance, and offering 15x speedup over zero-order methods. 

A SUPPLEMENTAL EXPERIMENTAL RESULTS:

Surrogate gradient methods We use feed-forward neural networks to model the surrogates since obtaining gradients for such networks is efficient thanks to automatic differentiation frameworks. We use networks comprised of 4 hidden layers with sizes 128, 72, 64, 32 and apply the same architecture to approximate the gradients for κ and σ separately. The hidden layers compute ReLU activation whereas sigmoid was used at the output layer (after the target output is properly normalized). For the Taylor-2 variant (in Eq. 18), we have an output vector for the gradient prediction. The networks are trained on the corresponding objective functions set up earlier by Adam optimizer with learning rate 10 -4 and decay 1.0. We fine-tune the networks with simple grid-search and select the best models for comparison.

As emphasized throughout, our focus is more on approximating the gradient rather than learning the true function. However, we need to somehow evaluate the surrogate models on how well they generalize on a hold-out test set. Like canonical regression problems, we compare the surrogate variants against each other using root mean square error (RMSE) on the test set. Figures 2 and 3 shows the results. The left figure shows RMSE for predicting κ and the right one shows RMSE for the von Mises stress σ. We can see that the Taylor-Reg generalizes better and also converges faster than Taylor-1 and Taylor-2 to target RMSE for κ, while all methods result similarly for σ prediction. This is reasonable because the objectives of Taylor-1 and Taylor-2 are not to optimize the mean square error, which we evaluate on here. Figure 3 shows the learning in terms of sample complexity. Again, Taylor-Reg outperforms Taylor-1 and Taylor-2 for κ prediction. In contrast, most models work similarly for σ regression, particularly when the training size is reduced to 50% (25K). Effectiveness of Z-Hermite learning Notice that Z-Hermite learning is not included in this comparison and as a surrogate model in the black-blox Langevin sammpling in Section 8. The reason is that apart from the usual sample pair (x i , y i ), we need the gradient ỹi (See Eq. 17). Since we can query the solvers, this gradient can only be estimated using finite difference. For both κ and σ in our experiment, obtaining such data is extremely expensive. As a consequence, we do not have the full results of the Z-Hermite model. Instead, we ran a separate study to show the effectiveness of Z-Hermite surrogate LMC on a smaller data with a lightweight OpenBTE version (0.9.55). The results in Table 3 shows the working of Z-Hermite learning in learning the gradient of κ(x). Here, the entropy is based nearest neighbor estimate to demonstrate the diversity of the pore centers in the unit square. With the (x p , y p )-coordinates of each pore p, the entropy estimate is given by: A hybrid algorithm between zero-order and Taylor-1 surrogate We can see in Tables 1, 2 and 3 the trade-off between computation and accuracy of our approach. While zero-order PLMC and ProxLMC can achieve the lowest thermal conductivity, their computational costs are prohibitive.

H = 1 n n i=1 log(n min j =i p i -p j ) + log 2 + C.

In contrast, deep surrogate models (including Taylor-Reg, Taylor-1) are far more time-efficient but slightly worse in terms of achieving the optimal κ. To mitigate the trade-off, we propose a simple hybrid method that combines the best of the zero-order and Taylor-1 surrogate models. The algorithm is shown in Figure A that alternates between using the gradient from the zero-order estimate and the gradient of the deep surrogate depending on whether taking this step would decrease the potential function (i.e. κ). We show and compare the achieved κ and running time in Table 3 . Examples of the samples generated by Zero-order PLMC, Taylor-1 PLMC and the hybrid method are also depicted in Figure 4 . The hybrid achieves the thermal conductivity that is lower than Taylor-1 PMLC while running almost 2x faster than zero-order PLMC. This suggests that the hybrid strategy offers a better trade-off in accuracy and computation. One way to further improve the hybrid is to collect the zero-order gradients while mixing and re-update the surrogate with Z-Hermite learning.

Algorithm 1 A hybrid PLMC algorithm alternating between zero-order and Taylor-1 surrogate gradients.

Train a network f θ (x) with Taylor-1 Randomly sample x 0 from the uniform p(x) Perform a Langevin dynamic step for t = 1, 2, . . . , K do if κ(x -η∇ x f θ (x)) < κ(x) then x ← P Ω (x -η∇ x f θ (x) + √ 2ηξ) else estimate ∇κ(x) using zero-order method update x ← P Ω (x -η ∇κ( B BACKGROUND ON MODELING NANOSCALE HEAT TRANSPORT At the nanoscale, heat transport may exhibit strong ballistic behaviour and a non-diffusive model must be used (Chen, 2005) . In this work we use the Boltzmann transport equation under the relaxation time approximation and in the mean-free-path (MFP) formulation (Romano & Grossman, 2015 )

EQUATION

where T (Λ) is the effective temperature associated to phonons with MFP Λ and direction ŝ; the notation . stands for an angular average. The coefficients α(Λ ) are given by

EQUATION

where K(Λ ) is the bulk MFP distribution. In general, such a quantity can span several orders of magnitude; however, for simplicity we assume the gray model, i.e. all phonons travel with the same MFP, Λ 0 . Within this approximation, we have K(Λ) = κ bulk δ(Λ -Λ 0 ). In this work we choose Λ 0 = 10 nm, namely as large as the unit cell, so that significant phonons size effects occur. With no loss of generality, we set κ bulk = 1 Wm -1 K -1 . Eq. 23 is an integro-differential PDE, which is solved iteratively for each phonon direction over an unstructured mesh (Romano & Di Carlo, 2011) . We apply periodic boundary conditions along the unit cell while imposing a difference of temperature of ∆T = 1 K along the x-axis. At the pores' walls we apply diffusive boundary conditions. Upon convergence, the effective thermal conductivity is computed using Fourier's law, i.e.

EQUATION

where J = (κ bulk /Λ 0 ) T (Λ 0 )ŝ n is the heat flux, L is the size of the unit cell, A is the area of the cold contact (with normal n). Throughout the text we use the quantity κ = κ eff /κ bulk as a measure of phonon size effects.

C BACKGROUND ON MODELING MECHANICAL STRESS:

We model mechanical stress by using the continuum linear elasticity equations

EQUATION

where f i is the body force (which is zero in this case), and σ ij is the stress tensor. Note that we used the Einstein notation, i.e. repeated indexes are summed over. The strain kl is related to the stress via the fourth-rank tensor elastic constant C ijkl

EQUATION

)

The strain is then related to the displacement u via kl = 1 2

EQUATION

We apply periodic boundary conditions along the unit-cell and applied solicitation is a small in-plane expansion. Once the stress tensor is calculated, we compute the von Mises stress as

EQUATION

)

where σ i are the principal stress axis. As a mechanical stability estimator we use σ = max x∈D (σ V M ) where D is the simulation domain. To avoid material's plasticity, σ needs to be smaller than the yield stress of a given material. For mechanical simulation we used the SUMIT code ( MIT Development Group, 2018). where X 0 = Y 0 is an F 0 measurable random variable. Suppose that the non-anticipative functionals α t (x) and β t (x) are such that a unique continuous strong solutions exits for both processes. If for any T > 0:

EQUATION

E BACKGROUND ON ZERO-ORDER OPTIMIZATION (GRADIENT-FREE)

Consider the smoothed potential U ν defined as follows:

U ν (x) = E g∼N (0,I d ) U (x + νg)

its gradient is given by:

∇ x U ν (x) = E g U (x + νg) -U (x) ν g, A monte carlo estimate of ∇ x U ν (x) is: Ĝn (x) = 1 n n j=1 U (x + νg j ) -U (x) ν g j ,

where g 1 , . . . g n are iid standard Gaussians vectors.

Using known results in zero order optimization under assumptions on smoothness and bounded gradients of the gradients we have for all x ( (Nesterov & Spokoiny, 2017; Shen et al., 2019) ):

E g Ĝ1 (x) -∇ x U (x) 2 ≤ βν(d + 2) 3/2 + (d + 1) 1 2 ∇ x U (x) 2 ≤ βν(d + 2) 3/2 + (d + 1) 1 2 L 2

Finally by independence of u 1 , . . . u n we have:

EQUATION

F PROOFS

Proof of Lemma 1. Define the Lagrangian:

L(q, η) = Ω log q(x) p 0 (x) q(x)dx + Ce j=1 λ j Ω (ψ j (x) -y j ) 2 q(x)dx + Ci k=1 λ k x∈Ω (φ k (x) -b k ) + q(x)dx + η 1 - x∈Ω q(x)

Setting first order optimality conditions on q, we have for x ∈ Ω:

log q(x) p 0 (x) + 1 + C j=1 λ j (ψ j (x) -y j ) 2 + Ci k=1 λ k (φ k (x) -b k ) + -η = 0

Hence we have:

q(x) = p 0 (x) exp - Ce j=1 λ j (ψ j (x) -y j ) 2 - Ci k=1 λ k (φ k (x) -b k ) + e exp -η

, x ∈ Ω and q(x) = 0, x / ∈ Ω, First order optimality on η give us: Ω q(x) = 1, we conclude by setting e exp(-η) = Z.

Proof of Theorem 1 1) Projected Langevin. Let us define following continuous processes by interpolation of X k and Y K (Piecewise constant):

d Xt = P Ω ( Ũt ( X)dt + √ 2λdW t )

where Ũt ( X) = - Let µ X T be the law of Xt∈[0,T ] . Same for µ Ỹ T . The proof here is similar to the proof of Lemma 8 in (Bubeck et al., 2015) . By the data processing inequality we have: 

G(Y kη ) -∇ x U (X kη ) 2 dt = η 4λ K-1 k=0 E G(Y kη ) -∇ x U (X kη ) 2 = η 4λ K-1 k=0 E G(Y kη ) -∇ x U (Y kη ) + ∇ x U (Y kη ) -∇ x U (X kη ) 2 ≤ η 2λ K-1 k=0 E G(Y kη ) -x U (Y kη ) 2 + E ∇ x U (Y kη ) -∇ x U (X kη ) 2

where in the last inequality we used the fact that ||a -b|| 2 ≤ 2(||a|| 2 + ||b|| 2 ). Note that we have by smoothness assumption on U :

∇ x U (Y kη ) -∇ x U (X kh ) 2 ≤ β 2 X kh -Y kh 2

Let R be the diameter of Ω, we can get a bound as follows:

EQUATION

Gradient approximation error

+β 2 K-1 k=0 E X kh -Y kh 2       ≤ η 2λ K-1 k=0 E G(Y kη ) -∇ x U (Y kη ) 2 + Kβ 2 R 2

Now using Pinsker inequality we have:

T V (µ X T , µ Ỹ T ) 2 ≤ 2KL(µ X T , µ Ỹ T ) ≤ η λ K-1 k=0 E G(Y kη ) -∇ x U (Y kη ) 2 + Kβ 2 R 2

Hence for T = Kη we have: where G t ( Ỹ ) = -∞ k=0 (G( Ỹkη ) + 1 γ ( Ỹkη -P Ω ( Ỹkη )))1 t∈[kη,(k+1)η] (t). Now applying Grisanov's Theorem for diffusions (Theorem 4) we have:

EQUATION

KL(µ X T , µ Ỹ T ) = 1 4λ E P X T T 0 U t ( X) -G t ( X) 2 dt = 1 4λ E K-1 k=0 (k+1)η kη G( Xkη ) -∇ x U ( Xkη ) 2 dt = η 4λ K-1 k=0 E G( Xkη ) -∇ x U ( Xkη ) 2 = η 4λ K-1 k=0 E G(X k ) -∇ x U (X k ) 2 .

Now using Pinsker inequality we have:

T V (µ T X , µ T Ỹ ) 2 ≤ 2KL(µ T X , µ T Ỹ ). Hence for T = Kη we have:

EQUATION

)

Proof of Theorem 2 . S-PLMC. If we set λ = 1, η ≤ α/K 2 , where α = 1/(δ + β 2 R 2 ), in this Corollary we obtain that : T V (µ S-P LM C K , µ P LM C K ) ≤ 1 √ K . Assuming A, B and C we consider η ≤ min(R 2 /K, α/K 2 ), and K = Ω(ε -12 d 12 ). Now using the triangle inequality together with the bounds in Eq.s 7 we have: T V (µ S-P LM C K , π) ≤ T V (µ S-P LM C K , µ P LM C K ) + T V (µ P LM C , π) ≤ ε + 1 √ K . S-ProxLMC. We conclude with a similar argument for T V (µ S-P roxLM C K , π) using Eq.s 8. Considering η = min(γ(1 + β 2 γ 2 ) -1 , 1 δK 2 ), and K = Ω(ε -6 d 5 ), we obtain (ε + 1 √ K ) approximation in TV of the target Gibbs distribution.

Proof of Corollary 1. Z-PLMC: We have:

T V (µ X T , µ Ỹ T ) ≤ η λ K-1 k=0 E G n U (Y kη ) -∇ x U (Y kη ) 2 + Kβ 2 R 2

Taking the expectation we have:

E g1...gn T V (µ X T , µ Ỹ T ) ≤ E g1...gn η λ K-1 k=0 E G n U (Y kη ) -∇ x U (Y kη ) 2 + Kβ 2 R 2 ≤ η λ K-1 k=0 E Y E g1...gn G n U (Y kη ) -∇ x U (Y kη ) 2 + Kβ 2 R 2 (Jensen inequality)

Note now that we have:

E g1...gn G n U (Y kη ) -∇ x U (Y kη ) 2 ≤ δ, ∀Y kη .

For n ≥ βν(d + 2) 3/2 + (d + 1)

1 2 L 2 /δ

The rest of the proof is an application of Theorem 2.

Z-ProxLMC. A similar argument holds.

:

Note that both properties κ and σ for a given configuration are obtained by numerically solving highly non-linear PDEs. The material configuration is defined by the pore locations, the material used, and the response of the material to heat (thermal) or stress (mechanical) flows.

Figures and Tables:

Remark 4. Assumption D on the -approximation of the gradient can be achieved for a large enough training set N , when we use Hermite learning in RKHS under mild assumptions and in Taylor learning. The assumption on the dominance of the training distribution is natural and means that we need a large training set that accounts to what we may encounter in Surrogate LMC iterations.

Figure 1: Example of nano-porous structures with corresponding heat flux shown using a color gradient. Yellow regions indicate high phonons flux. The thermal conductivity κ and von Mises stress σ are reported below each structure. The arrows show the moving directions of the pores. (Left) A random sample. (Middle) The sample obtained by Taylor-Reg PMLC starting from the left structure with κ constraint. (Right) The sample obtained by Taylor-Reg PMLC with both κ and σ constraints.

Figure 2: Comparison of the surrogate variants in testing RMSE. (Left) prediction accuracy for the thermal conductivity κ. (Right) prediction accuracy for mechanical stability σ. Note the difference in scale of κ and σ.

Figure 3: Comparison of the surrogate models in RMSE on the same test set when the training size is varied. Note the scale difference in the figures due to the different range of values.

Figure 4: Samples from by Zero-order PLMC (left), Taylor-1 PLMC (middle) and the hybrid algorithm of Zero-order and Taylor-1 PLMC (right). All are run with the κ constraint.

k=0 ∇ x U ( Xkη )1 t∈[kη,(k+1)η] (t). Similarly let us define :d Ỹt = P Ω (G t ( Ỹ )dt + √ 2λdW t )whereG t ( Ỹ ) = -∞ k=0 G( Ỹkη )1 t∈[kη,(k+1)η] (t).It is easy to see that we have :X k = Xkη and Y k = Ỹkη .Let π T X and π T Ỹ be the distributions of ( Xt ) t∈[0,T ] and ( Ỹ ) t∈[0,T ] . Note that :d Ỹt = P Ω Ũt ( Xt )dt + √ 2λ(dW t + 1 √ 2λ (G t ( Ỹt ) -Ũt ( Xt ))dt) Let d Wt = dW t + 1 √ 2λ (G t ( Ỹt ) -Ũt ( Xt ))dtHence we have :d Ỹt = P Ω Ũt ( X) + √ 2λd Wt , Assume that X 0 = Y 0 there exists Q such that , X T = Q({W t } t∈[0,T ]) and Y T = Q(( Wt ) t∈[0,T ] ).

KL(µ XT , µ Ỹ T ) ≤ KL(W t∈[0,T ] , Wt∈[0,T ] ), Now using Grisanov's Theorem for change of measure of Brownian Motion (Theorem 3) we have:KL(W t∈[0,T ] , Wt∈[0,T ] ) = 1 4λ E T 0 |G t ( Ỹt ) -Ũt ( Xt )| 2 dtConsider T = Kη, hence we have (with some abuse of notation we drop tilde as Y k = Ỹkη ):

Proof of Theorem 1 2) Proximal LMC. Let us define the following continuous processes by interpolation of X k and Y K (Piecewise constant):d Xt = Ũt ( X)dt + √ 2λdW t where Ũt ( X) = -∞ k=0 (∇ x U ( Xkη ) + 1 γ ( Xkη -P Ω ( Xkη )))1 t∈[kη,(k+1)η] (t). Similarly let us define : d Ỹt = G t ( Ỹ )dt + √ 2λdW t

. The standard BO protocol is comprised of estimating the black-box function from data through a probabilistic surrogate model, usually a Gaussian process, and maximizing an acquisition function to decide where to sample next. BO is often performed over a latent space, as in(Gómez-Bombarelli et al., 2018).Hernández-Lobato et al. (2016) proposed an information-theoretic framework for extending BO to address optimization under black-box constraints, which is close to current problem scenario. Genetic Algorithms (GA), a class of meta-heuristic based evolutionary optimization techniques, is another widely used approach for generating (material) samples with desired property(Jennings et al., 2019) and has been also used for handling optimization under constraints(Chehouri et al., 2016). However, GA typically requires a large number of function evaluations, can get stuck in local optima, and does not scale well with complexity. Finally, Zhou et al. (2019) has used deep reinforcement learning technique of Deep Q-networks to optimize molecules under a specific constraint using desired properties as rewards.

Table content: <table/>

Statistics of 20 new samples obtained by running different surrogate-based Langevin MCs on π with the thermal conductivity constraint (Eq. 21). We show the min and mean over the generated samples and the per-sample time. Initialized samples have min κ = 0.0619 and mean σ = 0.1268.

Table content: <table><tr><td>Model</td><td>Min κ</td><td>Mean κ</td><td>Per-sam. time (s)</td></tr><tr><td>Regression-PLMC</td><td>0.0757</td><td>0.1206 ± 0.0480</td><td>1055</td></tr><tr><td>Taylor-Reg-PLMC</td><td>0.0638</td><td>0.1196 ± 0.0495</td><td>899</td></tr><tr><td>Taylor-1-PLMC</td><td>0.0637</td><td>0.1278 ± 0.0610</td><td>852</td></tr><tr><td>Zero-order-PLMC</td><td colspan="2">0.0510 0.1093 ± 0.0271</td><td>14967</td></tr><tr><td colspan="2">Regression-ProxLMC 0.0646</td><td>0.1282 ± 0.0531</td><td>1107</td></tr><tr><td colspan="2">Taylor-Reg-ProxLMC 0.0712</td><td>0.1205 ± 0.0455</td><td>899</td></tr><tr><td>Taylor-1-ProxLMC</td><td colspan="2">0.0575 0.1297 ± 0.0543</td><td>874</td></tr><tr><td colspan="2">Zero-order-ProxLMC 0.0719</td><td>0.1112 ± 0.0363</td><td>14938</td></tr></table>

Summary statistics of 20 new samples obtained by our sampling method on π(x) with κ and σ constraints Eq. 22. The starting samples are reused from the single constraint case (min κ = 0.0759, mean κ = 0.1268, and mean σ = 0.8181; note that σ can be as high as 16.)

Table content: <table><tr><td>5, 6 and 7.</td></tr></table>

Qiang Wu, Justin Guinney, Mauro Maggioni, and Sayan Mukherjee. Learning gradients: Predictive models that infer geometry and statistical dependence. J. Mach. Learn. Res., 2010. Zhenpeng Zhou, Steven Kearnes, Li Li, Richard N Zare, and Patrick Riley. Optimization of molecules via deep reinforcement learning. Scientific reports, 9(1):10752, 2019.

Table content: <table/>

Z-Hermite learning is sample efficient at training time as well as for the Langevin sampling, but collecting the training set is prohibitive. Zero-order and deep surrogate can work in a hybird way and offer a better trade-off between accuracy and computation.

Table content: <table><tr><td>Model</td><td colspan="3">Mean κ Mean entropy Per-sam. time (s)</td></tr><tr><td>Zero-order PLMC</td><td>0.0676</td><td>1.960</td><td>3658</td></tr><tr><td>Taylor-1 PLMC</td><td>0.0988</td><td>1.745</td><td>253</td></tr><tr><td>Z-Hermite PLMC</td><td>0.0946</td><td>1.739</td><td>227</td></tr><tr><td colspan="2">Hybrid (Zero-order + Taylor-1) 0.0786</td><td>1.867</td><td>2136</td></tr></table>

Then the random process : d Wt = dW t -β t dt or written equivalently:Wt = W t -t 0 β s ds, is a Wiener process with respect to F t , t ∈ [0, T ]. Let P W T = L (W [0,T ]), and P W T = L ( W[0,T ] ) the densities are given by: Let (X t ) t≥0 and (Y t ) t≥0 dX t = α t (X)dt + dW t dY t = β t (Y )dt + dW t

Table content: <table><tr><td>dP W T dP W T</td><td>= exp</td><td colspan="3">T 0 β s , dW s -1 2</td><td colspan="2">T 0</td><td>β s</td><td>2 ds . It follows that:</td></tr><tr><td/><td colspan="2">KL(P W T , P W T ) =</td><td>1 2</td><td>E P W T</td><td>0</td><td>T</td><td>β s</td><td>2 ds</td><td>(30)</td></tr><tr><td colspan="9">Theorem 4 (Grisanov Theorem, Change of Measure for Diffusion Processes, (Lipster &amp; Shiryaev,</td></tr><tr><td>2001), ()).</td><td/><td/><td/><td/><td/><td/><td/></tr><tr><td colspan="9">D BACKGROUND ON STOCHASTIC DIFFERENTIAL EQUATIONS (SDE):</td></tr><tr><td colspan="9">CHANGE OF MEASURE AND GRISANOV'S FORMULA</td></tr></table>

Peer review (iclr2020):
The paper presents novel Langevin sampling implementations with (i) zero order and (ii) deep neural network gradient approximations for constraint satisfaction problems. These approximations improve the computational efficiency by a factor of 17 compared to related black box sampling approaches (Shen et al.). As a non-domain expert I cannot evaluate the importance of the contributions or the complexity of the tasks. However, given the well written and clearly structured presentation, the theoretical proofs and the generality of the constraint solver (with learned approximated gradients), I would vote for accepting the paper. Open questions: - How does the approach relate to *multi-objective Bayesian optimization*. The result of this multi objective optimization problem can be also used to further investigate different configurations in different tasks (generalization via pareto front sampling)? - The authors note that their approach is closely related to Shen et al. (2019). However, there is no algorithmic comparison in the experiments. Can the work of Shen et al. (2019) also be applied to the constraint setting? - How crucial or restricting is the assumption of log-concave and smooth target distributions? Which constraint satisfaction problems fall into this class? Can you provide some examples? - The authors use for approximating the gradients a neural network with 4 hidden layers (with 128, 72, 64 and 32 ReLU neurons) and fine tune the networks in a grid search manny. Which parameters were fine tuned? What were the input and output dimensions. I assume 10 kappa and 10 sigma values were the predicted outputs in individual networks. Minor Issues: - Some paragraph titles to no end with a dot, e.g., *Optimization approaches* in page 8 - The last sentence in page 8 is incomplete. - The computational speed up is a factor of 17 (page 9) or 15 (page 10)? - Typo in page 13: after the target output is properly normalization