Surrogate-Based Constrained Langevin Sampling With Applications to Optimal Material Configuration Design
Authors: Thanh V Nguyen
Youssef Mroueh
Samuel C. Hoffman
Payel Das
Pierre Dognin
Giuseppe Romano
Chinmay Hegde
Venue:  (2020)
Keywords: Black-box Constrained Langevin sampling
surrogate methods
projected and proximal methods
approximation theory of gradients
nano-porous 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 andshow that this constraint satisfaction problem can be formulated as sampling froma Gibbs distribution.  The main challenges come from the black-box nature ofthose physical constraints, since they are obtained via solving highly non-linearPDEs. 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.

Under review as a conference paper at ICLR 2020
SURROGATE-BASED CONSTRAINED LANGEVIN SAMPLING WITH APPLICATIONS TO OPTIMAL MATERIAL CONFIGURATION DESIGN
Anonymous authors Paper under double-blind review
ABSTRACT
We consider the problem of generating configurations that satisfy physical con- straints for optimal material nano-pattern design, where multiple (and often con- flicting) 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 thermoclectric efficiency. To that end, we leverage the posterior regularization framework and show that this constraint satisfaction problem can be formulated as sampling from Gibbs distribution. The main challenges those physical constraints, since they are obtained via solving highly non-linear PDEs. To overcome those difficulties, we introduce Surrogate- -base 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-Ordee 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 distri- bution is log-concave and smooth. We show the effectiveness both approaches in designing optimal nano moro-porous material configurations, where the goal is produce nano-pattern templates with low thermal conductivity and reasonable mechanical stability.
from the black-box nature of
INTRODUCTION
In many real-world design problems, the optimal design needs to simultaneously satisfy multiple constraints, which can be expensive 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 their inner variables their gradients.
We pose this design-uunder-constrainss 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 I 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 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

---

Under review as a conference paper at ICLR 2020
mixing of the Langevin dynamics towards the target Gibbs distribution. We devise several objective functions, as well as deep neural-netword architectures for parameterizing the approximating function class, for learning the gradient of the potential function. (ii) We then the surrogate gradient model in the constrained Langevin dynamics in lieu of the black-box potential. Using the surrogate enables more efficient sampling, since 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 the target Gibbs distribution.
3. We introduce Surrogate Model Based Projected Langevin Monte Carlo via learning the gradient 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 nano- porous 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 x) and desired mechanical stability (von Mises Stress a S T, where is some preset threshold).
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 of samples satisfying list of equality and inequality constraints: (x)) = j 1... Ce. and Φx (x) S by: k: 1... .. C, where E 52 and f2 IRd is bounded domain. We assume prior distribution pa (whose analytical form is known). The main challenge in black-box optimization is that the functions 1/ and øR can be only evaluated point-wise, and neither do we have functional forms nor access to their gradients. For example, 1h 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 Ay > 0 and obtain the following relaxed objective:
Ce
C (o
min KL.q.Po) od d()==1
+ Σ ( - )22) +
bx)+
(I)
y=1
A=1
The formulation in Eq. is similar in spirit 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 Q is bounded), and (ii) we assume black-box setting. We first obtain:
Lemma 1 (Constraint Satisfaction Sampling from Gibbs Distribution). The solution to the distribution learning problem given in Eq. is given by:
T(x) exp( exp(-U(xx) U(x)) 1 1Een
Z
where U(x) = -logp((z) + EE (x) yx)² + E Ae(o(()) bx)+ and Z Jeen exp (-U(z)) dr.
Lemma shows that the constraint satisfaction problem formulated in Eq. amounts to sampling from a Gibbs distribution defined on a compact support given in Eq. 2. Sampling from a Gibbs
Note that both properties k and a for given configuration are obtained by numerically solving highly non-linear PDFs. The material configuration is defined by the pore locations, the material used, and the response of the material to heat (thermal) stress (mechanical) flows.
2

---

Under review as a conference paper at ICLR 2020
distribution (also known as Boltzmann distribution) has 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 compact domain Ω and their mixing properties were actively studied Bubeck et al. (2015); Brosse et al. (2017). In the next Section, provide more detailed review.
Remark (Relation Bayesian Optimization). While 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 al., 2019) for more details.
For the rest of the paper, we will assume po to the uniform distribution on S2 which means that its gradients are zero on the support of the domain Q. Otherwise, if po is known and belongs to, for instance, an exponential family or generative model prior (such as normalizing flws), can sample from x using mixture of black-box sampling on the constraints (175, Vx)) and white-box sampling on log(Po)-
Remark 2.
3 WHITE-BOX SAMPLING: CONSTRAINED LANGEVIN DYNAMICS
We review in this section Langevin dynamics in the unconstrained case (S2 = Rd) and the constrained setting (22 C (R.). Below, ||-Il denotes the Euclidean norm unless otherwise specified. We are interested in sampling from
(T)) 1 exp(-U(<))
(3)
Preliminaries. We give here assumptions, definitions and few preliminary known facts that will be useful later. Those assumptions 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 f2 is a convex such that 0 E Q, Z contains a Euclidean ball of radius r, and Ω is contained in a Euclidean ball of radius R. (For example, Q might encode box constraints.) The projection onto Ω2 Pa(xr) is defined follows: for all Ω Pa(x) arg mingen || =1|²
Let R = supz,r'n en || - x'1l < cc.
2. Assumption B: We assume that U is convex, B-smooth, and with bounded gradients:
[V)U(x) v,U(W) s 8 | yll
Vx, y E Ω (3-smoothness).
[V(())) S L, Vo E Ω (Boundedness).
The Total Variation (TV) distance between two measures # is defined as follows: TV (p, v) sup [x(A) v(A)). Pinsker Inequality relates KL divergence to TV: TV (p, S (2KL(H,,) v).
Unconstrained Langevin Dynamics. In the unconstrained case, the goal sample from Gibbs distribution πT (2) exp(-U Uxx))/ that has unbounded support. This sampling can be done via the Langevin Monte Carlo (LMC) algorithm, which is given by the following iteration:
X Xx vXLn) + 22 = 1-....K11 (LMC),
where & ~ (O, Ta), 1y is the learning rate, and À a variance term.
Constrained Langevin Dynamics. In the constrained case, the goal is to sample from T(x) = exp(-U()))/21eer)) 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:
Xk+1* Ps (x
22XT(x)
k (PLMC),
(5)
In essence, PLMC consists of a single iteration of LMC, followed by a projection on the set using the operator Pa-
3

---

Under review as a conference paper at ICLR 2020
Proximal Langevin Dynamics. Similar proximal methods in constrained optimization, Brosse et al. (2017) introduced Proximal LMC (ProxLMC) that uses the iteration:
XR+1 (1
Xx- XV- " Px(x))
+ 22r = 1.....-1.. (ProxLMC) (6)
where 1j is the step size and y is regularization parameter. In essence, ProxLMC (Brosse al.. 2017) performs ordinary LMC on (x) = (x) + (z),, where (ar) is the proximal operator:
% (1) = inf (n + (2>) -1 || y/l 2 (2y) -1 | Pa(x))
where (x) 0 for x S2 and (z)
for Q. Therefore, the update in Eq. i regular Langevin update (as in Eq. 4) with potential gradient (x) (=)(z)++7) (r
(z))).
We denote by PK PLMC and fk PM the distributions of Xς obtained by iterating Eq. and Eq. respectively. Under Assumptions A and B. both these distributions converge to the target Gibbs distribution in the total variation distance. In particular, Bubeck al. (2015) showed that for 77 (R² [)), we obtain:
PLMC . S e for K = 12) TV (uK Likewise, Brosse et al. (2017) showed that for < S 7(1 + /3 -1, we obtain:
(7)
TV TV
(uK (u TT) S for K (0(=-dd"),
(8)
where the notation On = (3,),) means that there exists c E R. C > 0 such that on 2 CBBn (B,))
CONSTRAINED LANGEVIN DYNAMICS IN THE BLACK-BOX SETTING
We now introduce 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, 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
Rd be a vector valued function that approximates the gradient of the potential, U.U. We make:
Assumption C. The surrogate gradient G satisfies E [[G(Ya)] < <o, Vk.
Surrogate Projected Langevin Dynamics. Given Yo. the Surrogate Projected LMC (S-PLMC) replaces the potential gradient V in Eq. 5 with the surrogate gradient G:
Ys+1 Pn Yv GY2) + 22Xn(x) - K - (S-PLMC)
(9)
Surrogate Proximal Langevin Dynamics. Similarly, the Surrogate Proximal LMC (S-ProXLMC) replaces the unknown potential gradient V,U in Eq. 6 with the gradient surrogate G:
1 (1-2) _ Y - 7)G,) " Pa(e)) + 22 & 0... K - (S-ProxLMC) (10)
Ye+1
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 variant of Grisanov's 's Theorem for change of measures (Lipster & Shiryaev, 2001) and Pinsker's Inequality that bounds total variation in terms of Kullback-Leiblee divergence.
Theorem (S-PLMC and S-ProXLMC Mixing Properties). Under Assumption C, we have:
1. S-PLMC Convergence. Let PLMC be the distribution of the random variable X obtained by iterating PLMC Eq. 5, and p*x
SPPMC be the distribution of the random variable Yg obtained by
iteration S-PLMC given in Eq. 9. We have:
1
K-1
TV (1O |KK
S-PLMC PLMC <
V
Σ
[|((Y)) V.U(Y.) K3RR2
(11)
A==
4

---

Under review as a conference paper at ICLR 2020
2. S-ProxLMC Convergence. Let PeK PK M.R.d.MC. be the distribution of the random variable Xx obtained by iterating Mr Eg. 6, and jKK by iterating S-ProxLMC given in Eq. 10.
SS-PPoLLCC be the distribution of the random variable Yg obtained
We have:
Σ BII((x) V. :
TV TV(ux M-PoLLLCC
Pro .AK
? 2X
(12)
From Theorem 1, we see that it suffices to approximate the potential gradient (x)(()) (and not the potential (X)) in order to guarantee convergence of surrogate- based Langevin sampling. Using the triangle inequality, and combining Theorem and bounds in Eqs and we obtain:
Theorem 2. (Convergence of Surrogate Constrained LMC to the Gibbs distribution.) Under assump- tions A,B and C we have:
1. Assume in S-PLMC that there exists 8 > O such that B. B)))) V (VL)I < a, vk 2 0. Set À = 1, and 1) = Ommnn( R2 K, oz K²)) where 1/(6 8 R²)). Then for K = 122)
we have;
TV
S-PLMC
.T)) E.
2. Assume in S-ProxLMC that there exists > O such that E ||G(Xx)) - VaU(Xa)I S 6, vk 2 0.
Set X = L, and = min(7(1 + 3212)-1 oK?
Thenfor K = 52(E -dd) ) we have:
TV (uK
S. ProXLMC
TT) E.
5 ZERO-ORDER CONSTRAINED LANGEVIN DYNAMICS
In zero-order optimization (Nesterov & Spokoiny, 2017: Duchi et al., 2015; Ghadimi & Lan, 2013; Shen et al., 2019), considers the Gaussian smoothed potential Un defined as (x) (a-or(0.1a) (x + vg), and its gradient given by (zU,(z) Eg U(xtug)-U(x) g. The following a Monte Carlo estimate of V (z); N
U(x vgz) U(x))
G,U(s)
9j.
(13)
72
V
where g1.... In are i.i.d. standard normal vectors.
Zero-Order sampling from log- conccce densities was recently studied in Shen et al. (2019). We extend it here to the constrained sampling case of log- eggccoccve densities with compact support. We define Constrained Zero-Orda Projected LMC (Z-PLMC) and Zero-Order Proximal LMC (Z-ProXL.MC) by setting G(x) G, U(x) in Eq. and Eq. 10 respectively.
Lemma (Zero-Order Gradient Approximation & Spokoiny, 2017; Shen et al., 2019)).
Under Assumption B, we have for all Ω2:
Egi.....h [ < (Bvld+ 223/22222222 (d (d+1)+1)) /n
-
(14)
Thanks to Lemma that ensures uniform approximation of gradients expectation, we can apply Theorem 2 and get the following corollary for Z-PLMC and Z-ProxLMC:
Corollary (Zero-order Constrained Langevin approximates the Gibbs distribution). Under As- sumptions A and B, let 8 E [0, for n 2 (Bv(d 22/32 + (d+ 1)+ 15, we have the following bounds in expectation:
1. Set X = 1, and (min R2 /K,//²2) where a= = 1/(6 (6 + R²) For K = 12 we have:
Eg ....... TV (W (uK
ZPLMC
TT) SE
(15)
2. Set À = 1, and r = min(7(1 (1 + 22 2- 522). For = 52(E 6dD) we have:
Eg....g.. TV M-ProxLMC ,T) <
(16)
Remark 3. For simplicity, state the above bound in terms of expectations over the randomness in estimating the gradients. is possible to get Enite-sample bounds using the Vector Bernstein concentration inequality, coupled with covering number estimates of but omit them due space.
5

---

Under review as a conference paper at ICLR 2020
6 SURROGATE MODEL BASED CONSTRAINED LANGEVIN DYNAMICS
Despite its theoretical guarantees, zero-order constrained Langevin (Z-PLMC and L-ProxLMMC) has a prohibitive computation cost it needs (nnK) black-box queries (in our case, invocations of 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 depends f/y and Φx. 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 given a training set S for a PDE solver for the property &b (dropping the index j for simplicity):
S ( V4 (xx)), 8. Gn (xx))) paiid., i 1, N...
where Pn is the training distribution and G₂ (.) is the zero-order estimate of the gradient of a given in Eq. 13. We propose to learn a surrogate model belonging to function class Hs, fe E Me. that regresses the value of di and matches the zero-order gradient estimates as follows:
1
min E.E
N (v: fo(m:))2 (o. - ūe 1² (Z-Hermite Learning)
(17)
The problem in Eq. 17 was introduced and analyzed in Shi et al. (2010) where My is ball in 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 (2088;; Czarnecki et al. (2017) where My is a deep neural network parameterized with weights O. When So 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 Hi is an RKHS ball and when Ø; (2,((z);) are exact gradients, for sufficiently large training with N = ((1/2/(2r)) (where r; Ç are exponents in [0, 1] that depend on the regularity of the function v³). Under the assumption that v/s E He we have: fa IIV fe(()) V (=) (z()d < e. 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 nutshell, Mukherjee & Zhou (2006) suggests to learn a surrogate potential fo and gradient GA that are consistent with the first-order taylor expansion. Given a training set S = ((ri, Vi = (xxx)),, ~ io 11. N....N. Wu et al. (2010)
suggest the following objective:
min ToE.H.GGAEE GAE
Σ wgf(u so(2z) (GA(I.), -xa))((Taylor-.))
(18)
No
where w. = exp ( -112;=001 Mo is an RKHS ball of scalar valued functions, and Hd is an RKHS ball of vector valued functions.
Under mild assumptions, Mukherjee & Zhou (2006) shows that we have for N = 4/1) Jo [z) V=b(2))l pa(x)drr S c. We simplify the problem in Eq. 18 and propose the following two objective functions and leverage the deep learning toolkit to parameterize the surrogate fo:
min Jo
Σ w. fo(2z) + (V..olz;), (Taylor)),
(19)
N2

---

Under review as a conference paper at ICLR 2020
min 1 JoE N
A N 1, 2.)) (Taylor-Reg) (20)
Σ 1 f
The objective in Eq. uses a single surrogate to parameterize the potential and its gradient. The objective in Eq. 20 is similar in spirit the Jacobian matching formulation in the sense that it adds regularizer on the gradient of the surrogate to 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 fo 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) Vxfo(x) in the constrained Langevin dynamics in Eqs 9 and 10.
Both Hermite and Taylor learning come with theoretical guarantees when the approximation func- tion space is an RKHS under some mild assumptions the training distribution and the regu- larity of the target function db. In Hermite learning (Theorem in Shi et al. (2010)) we have: [x)](= S e for sufficiently large training set = O ((1/e1/((r)) (where
Ez-pa V fo(r) exponents S, E [0, depend regularity of v/). In Taylor Learning with the objective function given in Eq. 18 (Proposition
(x))] s e for N = 412)) 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 Pri
Wu al. (2010) we have: Ex~pa [[G
Assumption D: Assume have a learned surrogate G' on training distribution pn such that Eg-pp [G(x) (=U(2))] < c. Assume pa(x)) >> 0, Va E Q and that it is dominating mea- sure of Langevin (PLMC, S-PLMC, Prox-LMC, S-ProxLMC) ) intermediate distributions Pkx i.c. there exists C such that:
(x(x)= S Cpn(1), Vx E Vk = K. K- 1.
Under Assumption D, follows immediately that
(Xx+)) - VU(X |IG(x) (xll ma(x) pa(x)
pa(a) S Ce
and hence we can apply Theorem for == Cc, and we obtain e-approximation of the target Gibbs distribution in terms of total variation distance.
Remark 4. Assumption D on the c-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 Taylor learning. The assumption on the dominance of the training distribution is natural and means that we need large training set that accounts to what we may encounter in Surrogate LMC iterations.
In what follows we refer to surrogate constrained LMC, as X-PLMC or x-ProxI where x is one of four suffixes ([Z-Hermite, Taylor-2, Taylor-1, Taylor-Reg)),
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 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, that the density has a compact support and hence the need to appeal to projected LMC (Bubeck al., 2015) and Proximal LMC (Brosse et al., 2017). It is worth nothing that Hsich et al. (2018) introduced recently mirror Langevin sampling that can also be leveraged in framework.
7

---

Under review as a conference paper at ICLR 2020
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 dual embeddings (Dai et al., 201 8). Estimating gradients can be also performed using Stein's '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 al., 1998; Frazier, 2018). 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. BÔ is often performed over a latent space, in (Gómez-Bonid Bombarelli et al., 20.8). Hernández-Lobato et al. (2016) proposed an information-theoree black-box constraints, which is close to current problem scenario. Genetic Algorithms (GA), class of meta- hett-heuristic based evolutionary optimization techniques, is another widely used approach for generating (material) samples with desired property (Jennings 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. The advantage of our framework is that we obtain a distribution of optimal configurations (as opposed to 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.
framework for extending BO to address optimization under
8 EXPERIMENTS
In this section, we demonstrate the usability of our black-blox Langevin sampling approach for the design of mano-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 T 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 proximal update step.
Data. We want to learn surrogate models to approximate the gradient of the potential from data. To this end, we generate dataset of 50K nano-porous structures, each of size 100nm x 100nm. One such example is displayed in Fig. 1. Number of pores is fixed to 10 in this study and each pore is square with 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 ( Development Group, 2018), we obtain for each structure x pair of values: thermal conductivity x and von Mises stress a. Finally, we collect two datasets: {(T,, Kx) ) and {2,, oy { with the same inputs x; 's and N 50K samples. More details are given in Appendices B and C on the PDEs and their corresponding solvers.
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 aligmment 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.
Features.
Surrogate gradient methods. We use feed-fon hor-foorard 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 A and a 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
8

---

Under review as a conference paper at ICLR 2020
6 = 0.0871, = 0.4826
K = 0.0802, a = 0.6681
R. = 0.0732, = 0.4401
Figure 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 a 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 x constraint. (Right) The sample obtained by Taylor-Reg PMLC with both x and o constraints.
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 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 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 x 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:
क x④x) po(r) exp( )x)
(21)
7
where Po(x) is the uniform distribution over the unit square, which is equivalent to the Poisson process of 10 pores on the square, and k(x) is the thermal conductivity we want to minimize. Starting from 20 samples initialized from (r(z), we run our proposed black-box Langevin MCs and obtain 20 new realizations from the target distribution (xx). We use four different surrogates (including simple regression, Taylor-Reg, Taylor-I 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-PMICC in the first row and regression-Pront ProxLLCC in the fifth represent the sampling where the surrogate model are fitted on solely the mean square error objective. step size 7 = le-3 and the exponential decay rate 0.8. Since keeping track of the true value expensive, we stop after K 10 iterations. We first observe that the regression- (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- demonstrate its effectiveness in approximating the gradient and are able to achieve lower thermal conductivity. In particular, Taylor- I-CroXLCC and ceee--rrerr--MMM perform in the similar range in terms the minimum achieved, but the learned surrogate offers 17x speed up (per sample) to the space limit, we do not report Taylor-2 results in Table 1, and note that Taylor-2 works in the similar vein T
all methods, we set 100, the
method
zero order methods. Due

---

Under review as a conference paper at ICLR 2020
Model
Min x
Mean K
Per-sam. time (s)
Regression-PL.Ne
0.0757
0.1206 + 0.0480
1055
Cioor-Ree-PLMG
0.0638
0.1196 + 0.0495
899
Taylor-IPPLMC IPLMC
0.0637
0.1278 + 0.0610
852
7CMC 0.0510
0.1093 + 0.0271
14967
Regression-ProxxLMee
0.0646
0.1282 + 0.0531
1107
Taylor-Reg-ProxLLMe
0.0712
0.1205 + 0.0455
899
TTylor- 1PMoxLC
0.0575
0.1297 = 0.0543
874
Pemoo-rrdr--PrxxLM.. 0.0719
0.1112 + 0.0363
14938
Table I: Statistics of 20 new samples obtained by running different surrogate-basedd 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 x 0.0619 and mean o 0.1268.
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 cano-configurations. Hence, we consider both thermal conductivity x and mechanical stability provided via Mises stress We want sample that minimizes k(z) to achieve high thermo-electric efficiency while maintaining o (2) less than some threshold (which we explain below). Like the single constraint case, we pose this sampling from the following Gibbs distribution:
π(x) = po(r) exp(- exp(=>in(x)² {2(a(x) - -]++) A
(22)
Z
where po(x)) is the same above, (I) is the von Mises stress and is threshold on the maximum value of a. 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 T = 0.5, A = 100, A2 10 the step size 1 1e-3 and the exponential decay rate 0.8. Comparing with Table 1, one can see that not only better x. be achieved but also the a can be reduced simultaneously. These effectively sample new configurations under multiple competing constraints. Examples of new nano-configurations are show in Fig. and Appendix A Fig. 5, 6 and 7.
results suggest that our approach
Model
Min K
Mean K
Mean 0
Per-sam. time (s)
Taylor-Reg PLMC
0.0613
0.1256 + 0.0538
0.6590 + 0.2261
952
Taylor- I-PLMC
0.0611
0.1278 + 0.0610
0.6380 = 0.9898
852
Zero-order- PLMC
0.0471
0.1148 + 0.0475
0.6511 = 0.19
15677
Taylor-Reg- Pro
0.0666
0.1195
+ 0.0534
0.6402 + 0.1464
856
Taylor- Taylor-I-Prox 1-orox--PPxMM
0.0548
0.1298
0.0610
0.6156
0.1463
972
Pero-order-ProxLMMG 0.0354
0.1080 = 0.0384
0.6029 - 0.1376
15080
Table 2: Summary statistics of 20 new samples obtained by our sampling method on π(x) with and constraints Eq. 22. The starting samples reused from the single constraint case (min K = 0.0759, mean ń = 0.1268, and mean a = 0.8181; note that a can be as high as 16.)
9 CONCLUSION
In this paper we introduced Surrogate- Based Constrained Langevin Sampling for black-box sampling from a Gibbs distribution defined on 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 mano-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 prohibitive computational cost. We established that the deep the surrogate (Taylor-1
10

---

Under review as a conference paper at ICLR 2020
ProxLMC) is viable alternative zero-order methods, achieving reasonable performance, and offering 15x speedup over zero-order methods.
REFERENCES
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 Learning Theory, pp. 319-342, 2017.
Sebastien Bubeck, Ronen Eldan, and Joseph Lehec. Finite-time analysis of projected langevin monte carlo. In Advances in Neural Information Processing Systems, pp.
2015.
Adam Chehouri, Rafic Younes, Jean Perron, and Adrian Ilinca. A coonstraint-handiie technique for genetic algorithms using a violation factor. arXiv preprint arXiv:1610000997, 2016.
Gang Chen. Nanoscale energy transport and conversion: a parallel treatment of electrons, molecules, phonons, and photons. Oxford University Press, USA, 2005. URL https: /ttp:p::://ww///// /www. amazon. com/ Trr org
dd/011955942X.
Ming-Hui Chen and Bruce Schmeiser. Toward black-box sampling: A random-direction interior-point markov chain approach. Journal of Computational and Graphical Statistics, 7(1):1-22, 1998.
Wojciech M. Czarnecki, Simon Osindero, Max Jaderberg, Grzegorz Swirszez, and Razvan Pascanu. Sobolev training for neural networks. In Advances in Neural Information Processing Systems 30, pp. 02778-287... 2017.
Bo Dai, Hanjun Dai, Arthur Gretton, Le Song, Dale Schuurmans, and Niao He. Kernel exponential family estimation via doubly dual embedding, arXiv preprint 1rXiv:1811.02228. 2018
A. Dalalyan and A. B. Tsybakov. Sparse regression learning by aggregation and langevin monte- carlo. J. Comput. Syst. Sci., 2012.
Arnak S Dalalyan. Theoretical guarantees for approximate sampling from smooth and log-concave densities. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 79(3): 651-676, 2017.
John C Duchi, Michael Jordan, Martin J Wainwright, and Andre Wibisono. Optimal rates for zero-order convex optimization: The power of two function evaluations. IEEE Transactions on Information Theory, 61 (1) 2015.
Alain Durmus, Eric Moulines, et al. High- digh--dimessional bayesian inference via the unadjusted langevin algorithm. Bernoulli, 25(44)):854-2882 2019.
Peter Frazier. tutorial on bayesian optimization. arXiv preprint /rXiv: 2018.
Kuzman Ganchev, Jennifer Gillenwater, Ben Taskar, et Posterior regularization for structured latent variable models. Journal of Machine Learning Research, (1(uu))2001-2049, 2010.
Saeed Ghadimi and Guanghui Lan. Stochastic first-and zeroth-order methods for nonconvex stochastic programming. SIAM Journal on Optimization, (3(4))2341-2368 2013.
Rafael Bómezz-Bm Jennifer N Wei, David Duvenaud, José Miguel Hernández-L Lengeling, Dennis Sheberla, Jorge Ryan P Adams, and Alán Aapurr-Guzzk... Automatic chemical design using data-driven continuous representation of molecules. ACS central science, 4(2):268 276, 2018.
Benjamín Sánchez-L Sánchez-Lengeling,
Timothy D Hirzel,
Jun Han and Qiang Liu. Stein variational gradient descent without gradient. In Proceedings of the 35th International Conference on Machine Learning, pp. 1900-1908, 2018
José Miguel Hernández- Lobato. Michael A. Gelbart, Ryan P. Adams, Matthew W. Hoffiman, and Zoubin Ghahramani. A general framework for constrained bayesian optimization using information- based search. Journal of Machine Learning Research, 17(160):1-55 2016. URL http: :/ /mlr.. orrg http:/pprss/ww1//11 61 html.
11

---

Under review as a conference paper at ICLR 2020
Ya- Ping Hsieh, Ali Kavis, Paul Rolland, and Volkan Cevher. Mirrored langevin dynamics. In Advances in Neural Information Processing Systems. pp. 2878 277882887222 2018.
Zhiting Hu, Zichao Yang, Ruslan R Salakhutdinov, LIANHUI Qin, Xiaodan Liang, Haoye Dong, and Eric Xing. Deep generative models with learnable knowledge constraints. In Advances in Neural Information Processing Systems, pp. 10501-10512, 2018.
Aapo Hyvärinen. Estimation of non-normalized statistical models by score matching. J. Mach. Learn. Res., 2005.
Paul C Jennings, Steen Lysgaard, Jens Strabo Hummelshoj, Tejs Vegge, and Thomas Bligaard. Genetic algorithms for computational materials discovery accelerated by machine learning, npj Computational Materials, 5(1); 2019.
Donald R Jones, Matthias Schonlau, and William J Welch. Efficient global optimization of expensive black-box functions. Journal of Global optimization, 13(4):455-49 1998.
Yingzhen Li and Richard E Turner. Gradient estimators for implicit models. arXiv preprint 17X.v:1705.07107 2017.
Robert Lipster and Albert Shiryaev. Statisties of random processes. Springer, 2001.
Jonas Mockus. Application of bayesian approach to numerical methods of global and stochastic optimization. Journal of Global Optimization, 4(4):347-365, 1994.
Sayan Mukherjee and Qiang Wu. Estimation of gradients and coordinate covariation classification.
J. Mach. Learn. Res., 2006.
Sayan Mukherjee and Ding-Xuan Zhou. Learning coordinate covariances via gradients. J. Mach. Learn. Res., 2006.
Radford M Neal. Slice sampling. The Annals of Statistics, 31, 2003.
Yurii Nesterov and Vladimir Spokoiny. Random gradient Found. Comput. Math., 2017.
tradient-free free minimization of convex functions.
Giuseppe Romano and Aldo Di Carlo. Multiscale electrothermal modeling of nanostructured devices. IEEE Trans. Nanotechmol, (006)) 2011. URL htt.:ttteeeepepl.ooco ieee org/ http://////000////ammbeeeeeeee
Giuseppe Romano and Jeffrey C Grossman. Heat conduction in nanostructured materi- als predicted by phonon bulk mean free path distribution. J. Heat Transf.. 137(7): 071302, 2015. URL P a
asme.
Giuseppe Romano and Jeffrey C Grossman. Tuning thermal transport in disordered porous materials via phonon bottleneck identication. Phys. Rev. 1:3434, 2016.
Lingqing Shen, Krishnakumar Balasubramanian, and Saeed Ghadimi. Non-asymptotic results for langevin monte carlo: Coordinate-wee and black-box sampling. arXiv preprint 900.v:190..01333, 2019.
Lei Shi, Xin Guo, and Ding-Xuan Zhou. Hermite learning with gradient data. J. Comput. Appl. Math., 2010.
Suraj Srinivas and François Fleuret. Knowledge transfer with jacobian matching. arXiv preprint arXiv: (803.00443. 2018.
The MIT Development Group. mii, a scalable computational framework for large-scale simula- tion of complex mechanical response of materials, 2018. URL http:// //smm//smmmmm i t. mi t. edu.
Shinya Suzuki. Shion Takeno, Tomoyuki Tamura, Kazuki Shitara, and Masayuki Karasuyama. Multi- objective bayesian optimization using pareto-frontic entropy. arXiv preprint 17X1v:1966.00127, 2019.
12

---

Under review as a conference paper at ICLR 2020
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 Zare, and Patrick Riley. Optimization of molecules via deep reinforcement learning. Scientific reports, 9():: 2019.
A SUPPLEMENTAL EXPERIMENTAL RESULTS
Surrogate gradient methods We use forward neural networks 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 A and a separately. The hidden layers compute ReLU activation whereas sigmoid was used at the output layer (after the target output is properly normalized). For the T 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-1 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 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 and shows the results. The left figure shows RMSE for predicting and the right one shows RMSE for the von Mises stress a. We can see that the Taylor- generalizes better and also converges faster than Taylor-1 and Taylor-2 to target RMSE for K, while all methods result similarly for prediction. This is reasonable because the objectives of Taylor-l 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 I and Taylor-2 for prediction. In contrast, most models work similarly for regression, particularly when the training size is reduced to 50% (25K).
Performance in RMSE of the surroçate model for
Parformance in RMSE of the surrogate model for a
Nogresion
Rog ussin Taylr- Taybr)
Tayler Rog
Tayler
Taylor-2
Taylo12
DISSE (100 kmm
0.0
-12
-13
-0.2
-14
-0.4
-15
-0.6
10
15
20
25
30
Epects
Spocna
Figure 2: Comparison of the surrogate variants in testing RMSE. (Left) prediction accuracy for the thermal conductivity k. (Right) prediction accuracy for mechanical stability Note the difference in scale of and
Effectiveness of Z-Hermite learning Notice that Z-Hermite learning is not included in this com- parison and as a surrogate model in the black-blox Langevin sammpling in Section 8. The reason that apart from the usual sample pair (x,, v:). we need the gradient y; (See Eq. 17). Since we can query the solvers, this gradient can only be estimated using finite difference. For both k and o 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 separate study to show the effectiveness of Zeremmitt 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 k(x). Here, the entropy is based nearest neighbor estimate to demonstrate the diversity of the pore centers in the unit square. With the (Ip, Up) cooordiatss of each pore p. the entropy estimate is given by:
n log(n n min lps -Pyll) + log 2 C. a
12
13

---

Under review as a conference paper at ICLR 2020
Sample complexity of the surrogate models for
Sample complexity of the surrogate models for
Regular
-1500
Teylo-foz
-5.20
Tayrer Rog
Tople13
Toplor2
-0.25
Teyrer-
Tayler22
N902 Jogantomn
2.35
-4.35
-3.40
-2.40
-4.50
-1.50
-2.15
12
25
20
25
20
33
20
15
20
30
25
Numberaa tamples D00)
Number samples (X3800)
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.
Model
Mean K
Mean entropy
Per-sam. time (8)
Zero-order PLMC
0.0676
1.960
3658
Taylor-1 1 PLMC
0.0988
1.745
253
7-Hermite PLMC
0.0946
1.739
227
Hybrid (Zero-order Taylor-1)
0.0786
1.867
2136
Table 3: Z-Hermite learning is sample efficient 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 better trade-off between accuracy and computation.
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 PLMM and ProxLMC can achieve the lowest thermal conductivity, their computational costs are prohibitive. In contrast, deep surrogate models (including Taylor-Reg, Taylor-I) are far more time-efficient but slightly worse in terms of achieving the optimal To mitigate the trade-off, propose simple hybrid method that combines the best of the zero-order and Taylor-1 surrogate models. The algorithm 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. k). 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 The hybrid achieves the thermal conductivity that is lower than Taylor- PMLC while running almost 2x faster than zero-order PLMC. This suggests that the hybrid strategy offers better trade-off 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 A hybrid PLMC algorithm alternating between zero-order and Taylor- surrogate gradients.
Train a network fo fo(x) with Taylor-1 Randomly sample Zig from the uniform p(x) Perform a Langevin dynamic step
for t 2, 2,.,..... K do if x( - V (z)) < k(x) then x Pa(x (z) 22n))
else
estimate V*(x)) using zero-order method update <- Pa VK(x) + 22n
end if
end for
Return new sample
14

---

Under review as a conference paper at ICLR 2020
K = 0.063
K 0.090
R 0.076
Figure 4: Samples from by Zero-order PLMC (left), Taylor-1 PLMC (middle) and the hybrid algorithm of Zero-order and Taylor- PLMC (right). All are run with the constraint.
Additional generated samples We show additional configurations generated by our sampling ap- proach (Taylor-Reg ProxLMC, Taylor-1 ProxLMC and Zero-order ProxLMC) in Fig. 5, 6 and
0.0871, a 0.4826
0.0834, 0.4638
0.0942, 0.4566
Figure 5: Example of nano-porous structures with corresponding heat flux shown using color gradient. Yellow regions indicate high phonons flux. The thermal conductivity x and von Mises stress reported below each structure. The arrows show the moving directions of the pores. (Left) random sample. (Middle) The sample obtained by Taylor-Reg ProxLMC starting from the left structure with constraint. (Right) The sample obtained by Taylor-Reg ProxLMC with both and constraints.
15

---

Under review as a conference paper at ICLR 2020
6 = 0.0871, a = 0.4826
K = 0.0834, 0 = 0.4638
A. = 0.0942, 0 = 0.4566
Figure 6: 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 a are reported below each structure. The arrows show the moving directions of the pores. (Left) A random sample. (Middle) The sample obtained by Taylor- ProxLMC starting from the left structure with x constraint. (Right) The sample obtained by Taylor- ProxLMC with both x and constraints.
K = 0.0871, o 0.4826
K = 0.0719, o = 0.4424
K = 0. = 0.5110
Figure 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 a are reported below each structure. The arrows show the moving directions the pores. (Left) A random sample. (Middle) The sample obtained by Zero-order ProxLMC starting from the left structure with x constraint. (Right) The sample obtained by Zero-order ProxLMC with both k and o constraints.
B BACKGROUND ON MODELING NANOSCALE HEAT TRANSPORT
At the nanoscale, heat transport may exhibit strong ballistic behaviour and 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-patee -atch (MFP) formulation (Romano & Grossman, 2015)
As VT(A) +T()) = ((A))((A)))N))
(23)
where Т(A) the effective temperature associated to phonons with MFP and direction š; the notation (-) stands for angular average. The coefficients a<A') are given by
16

---

Under review as a conference paper at ICLR 2020
a(A) K(N) A / K(A")
dA"
(24)
A"
where (A') 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, Ag. Within this approximation, we have K (A) Kbulk Ao). In this work we choose Ao 10 nm, namely as large as the unit cell, so that sigmificant phonons size effects occur. With no loss of generality, we set Kbulk Wm -1 -K K-1 Eq. 23 is an integro-differente PDE, which is solved iteratively for each phonon direction over unstructured mesh (Romano & Di Carlo, 201 1). We apply periodic boundary conditions along the unit cell while imposing a difference of temperature of AT = 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.
ATA J. J-ndS, Add,
L
Koff
(25)
where J = (No (T(n)))))) the heat flux, is the size of the unit cell, A is the area of the cold contact (with normal fi). Throughout the text we use the quantity k Keff /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
dz, 0ij fa:
(26)
where fa the body force (which is zero this case), and Biy the stress tensor. Note that we used the Einstein notation, i.c. repeated indexes are summed over. The strain egt is related to the stress via the fourth-rank tensor elastic constant Ciga!
Øif = CiJREERI
(27)
The strain is then related to the displacement u via
1 Oux
dui
(28)
2
ax
Ouk
We apply periodic boundary conditions along the unit-cell and applied solicitation is small in in-plane expansion. Once the stress tensor calculated, we compute the von Mises stress as
Ov v / (os o2)² (as a1)² ( (02- o1)²
(29)
where o; are the principal stress axis. As mechanical stability estimator we use = m where D is the simulation domain. To avoid material's plasticity, o needs to be smaller than the yield stress of given material. For mechanical simulation we used the SUMIT code MIT Development Group, 2018).
D BACKGROUND ON STOCHASTIC DIFFERENTIAL EQUATIONS (SDE): CHANGE OF MEASURE AND GRISANOV'S FORMULA
Theorem 3 (Grisanov Theorem, Change of Measure for Brownian Motion (Lipster & Shiryaev, 2001), Theorem 6.3 page 257). Let (W,, Fi) be Wiener process (Brownian motion) and (Br, F,) random process such that for any T
1'" 1|&,112 dt < co a.s
17

---

Under review as a conference paper at ICLR 2020
Then the random process dW, dW, Bidt written equivalently: W, W, - So 3yds, is Wiener process with respect to F. t E [0, Let PT = (..1)) and PK - 2( (0.11) the
densities are given by: aply = exp (Jo' (Bs, 3 113.112 ds) It follows that:
aPw
KL(PT Pr E pr [ 11B&11 ds
(30)
Theorem 4 (Grisanov Theorem, Change of Measure for Diffusion Processes, (Lipster & Shiryaev, 2001), O). Let (Xt)120 and (Y))e20
dx, = dx(X) + dW,
dY, (Y dW.
where Xo = Yo is an Fo measurable random variable. Suppose that the non-anticipative functionals (z(2) and (x( are such that a unique continuous strong solutions exits for both processes. If for any T > 0:
ds < co(a.s) and S'
< co(a.s).
Let PX (X(oo.)) and P. (YYo.))
dPY (X) exp (-/-- (aa(x) Bx(),a()), dx.) + [' [[]a.(x)1-1188xxx
dP
L/ (xx)) ds
KL(PX,PH - Epp
(31)
E BACKGROUND ON ZERO-ORDER OPTIMIZATION (GRADIENT-FREE)
Consider the smoothed potential U, defined as follows:
U²(2) (o
its gradient is given by:
U(x vg) U(+)
V (x) R.
9.
monte carlo estimate VZUn(x) is:
C,(())
1 n
U(x vgg) U(x))
95>
where 91.... In. 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 ((Nesterov & Spokoiny, 2017; Shen et al., 2019)):
Eg - V.U(r)) < (Buld+ + (d 1 V=U(()))) < (Su(d+ 2)3/2 (d + 1) 1)+2)
(((z))
Finally by independence of uI, Un have:
BBlld (old++2)3/2 + (d (d+1) 1)²
Ey......n GG((z) -(.(z)--v.U()))
v.U('))
(32)
n
18

---

Under review as a conference paper at ICLR 2020
F PROOFS
Proof of Lemma 1. Define the Lagrangian:
Ce
L(q, 1) = S log
((x)
((r)drr + Σ
( (15(z)) n13 (x())))))))))))))))))))) (())dd
po(x))
5=1
C;
(du(x) bb)) + (1-/seh JuEN g(0)
Ak
Setting first order optimality conditions g, we have for x E s2:
((x)
log
- + +1+ Σ -1) + (x( bk)+ 17
(x
k=1
Hence we have:
exp (- EE (z)))(z)) 2; - {c(dx(0)- - -b2)+)
((x) PT(
, $2
ecxp -1
and
g(x) 0, $ Ω
First order optimality on give us: Jn g(x) = 1, we conclude by setting e eex((nn) = Z.
Proof of Theorem Projected Langevin. Let us define the following continuous processes by interpolation of X and Yx (Piecewise constant):
dX. Pa(Ui + V2XdW,)
where = V kn (1+1)) (t). Similarly let define
d². = Pn.G. + v2XWW,)
where - - G(Y,n) ( (1)-
It is easy to see that have X2 Xxn and Yv Ykn-
Let nT and oे be the distributions of (Xi)ee(o.T) and (Y)e
Note that:
d', Pn (ūi(Xx)dtt v2X(WW,+ + 1 (Go(Y)) (.(()))
V2A
Let
dW, dW, (G.(Y,) (,(X)))
22A
Hence we have
dÝ Po (Ú,(X)+ + v2AdW.).
Assume that Xo- Yo there exists Q such that XT= Q((W:) and Yr Q((W)
Let yux be the law of X Same for p.y. The proof here is similar to the proof of Lemma 8 in (Bubeck et al., 2015). By the data processing inequality we have:
KL x < KL(WiE Wie(o,r))
Now using Grisanov's Theorem for change of measure of Brownian Motion (Theorem 3) we have:
4A = [G.(Y) (x))))
KL(WW. Wie(o.r)
19

---

Under review as a conference paper at ICLR 2020
Consider T Kn,, hence we have (with some abuse of notation drop tilde Y₁ Ykn)):
Kn
KL.(H) x .1x)
-
E
[G.(V)
(X.))
4X
K-1
(k+1)1
E
[[G(Yun) VaU((X dt
4X
k=0
kn
E [((()))) VXU((l)))
41
E (|G(())))))
V.U((on) (/U()))) -V.UXx))
4X
Σ
(||G((
V + (xk(X)))
2A
km
where in the last inequality we used the fact that ||a - b11 $ 2([|a|1 + [[b[1)). Note that we have by smoothness assumption on U:
[V (UU(N))) s Yeall
Let R be the diameter of Q we can get a bound as follows:
K-1 E B|((())
K --1
(L (1 -
(Y 4.5
E
Mun - YYall?
2A
A==
km
crror Gradient approximation error
K
17
E E|(((Y))) - vVU(()))) K32
Now using Pinsker incquality we have:
K-1
TV < ( S
Σ
B|G((Y)) V.U - KBRR
#=0
Hence for T Kn we have:
TV (EK JKK
S-PLCC PLMC
S V
Bl|G(Y)) v.U(K)1 -ICCX)--TUUX]]]+3331))
(33)
Proof of Theorem 1 2) Proximal LMC. Let us define the following continuous processes by interpo- lation of X₃, and YK (Piccewise constant):
dx: Üi(X)dt + V2XdW,
where Ü,(8) - Ee V. + (Xkn define
Similarly let us
dř, G.(P))) ddl 22XWW,
20

---

Under review as a conference paper at ICLR 2020
where GY(Y) -
+
Now applying Grisanov's
Theorem for diffusions (Theorem 4) we have:
(01(8)
KL(MT, m4) = 4A
Epp
K-1
(((xn)) - V-U(ien dt (k+1)m
E
4X
k=0
Σ
E
- VaU((())))
K-1
Σ
[G(xx)
VzU(Xxa]1*
41
k=
Now using Pinsker inequality we have:
1z)² S (
Hence for T Kn we have:
K-1
(V
MC
,KK
MC)
S
V
2X
π
Bl|G(X+) - -VaUx)))
(34)
km0
Proof of Theorem 2. S-PLMC. If we set = 1, S a/K², where a 1/(6 + B²R)). in this Corollary we obtain that TV SPSS-PMMC
PLMC .1K (.IMC) S VK Assuming A, and we consider
?) S min(R / (²),,, and K (22 Now using the triangle inequality together with the S-PLMC
bounds in Eq.s we have: (1x ( n)) ( .u PLMC) TV (LMMC m) S
S-PLMC
E + VK*
L-ProxLMC. We conclude with similar argument for T ( .) using Eq.s Consid- ering 1 = min(y(1 + - अ and K 52(= 9d²), (((), we obtain (a + vR approximation in TV of the target Gibbs distribution.
Proof of Corollary 1. Z-PLMC: We have:
K-1
TV (PT &
"
EIGG,UYYm) V.U KB
Á
k= 0
Taking the expectation we have:
Egg...... TV ( S F Ry....g
R - V R2
K- Ey .n ||,,U(YYn,)) - v.U(),)))) + K3
(Jensen inequality)
k=
Note now that we have:
Egi...g. [[,,U(Yknn)) v.U(Y,,))1 < VY
For 2 (Bu(d= - 2)3/22 (d + 1)² 18 The rest of the proof is an application of Theorem 2.
Z-ProxLMC. similar argument holds.
21


Reviews:
[{"review_id": "ryxGJnWP9S", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "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": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "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.). \nAs 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. \n\nOpen questions:\n- 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)?\n\n- 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? \n\n- 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?  \n\n- 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. \n\nMinor Issues:\n- Some paragraph titles to no end with a dot, e.g., \"Optimization approaches\" in page 8\n- The last sentence in page 8 is incomplete. \n- The computational speed up is a factor of 17 (page 9) or 15 (page 10)?\n- Typo in page 13: after the target output is properly normalization", "extra_scores": {"experience_assessment": "I do not know much about this area.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper at least twice and used my best judgement in assessing the paper.", "review_assessment:_checking_correctness_of_experiments": "I assessed the sensibility of the experiments.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I assessed the sensibility of the derivations and theory."}, "extra_text": {}}, {"review_id": "S1l9L6uQqS", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "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": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The paper considers the problem of sampling points from a constrained set in R^d where the constraints can only be accessed in a zero order fashion. They consider the specific situation where the constraints are a solution of a complicated PDE solver and hence the derivatives or specific functional forms of the constraint cannot be obtained. They repose the problem as sampling from a Gibbs distribution whose potential contains constraints as penalties in a Lagrangian fashion. They now wish to sample from the Gibbs distribution using Langevin diffusion.  The Langevin process requires a derivative of the gradient. The setting does not allow for that and therefore the authors propose two approaches - \n\n1. Constructing the gradient from zeroth order entries of a gaussian smoothed potential (much like works of Nesterov et al on zero order optimization). \n2. Using a parameteric function class (like an RKHS or a neural network) to learn a function which well approximates the gradient of the constraints as well given zeroth order constraint evaluations. \n\nFor the latter, the authors propose two approaches. Hermite learning which directly approximates the error on the gradient as well as the function evaluation. However this involves a separate estimate of the gradient (via a zero order approach such as Gaussian smoothing). An alternative which seems more sample efficient is to do direct learning from zero order samples by penalizing first-order Taylor approximation between given points. \n\nThe authors provide a theoretical analysis of the all the approaches and sub-approaches given above and further provide experiments to evaluate this on a problem from material design. I am not at all familiar with the domain of the experiments to comment on competing approaches. As far as I can see the authors also compare only their own proposed algorithms which are many. I will focus on the theoretical analysis which seems to form the bulk of the paper anyway. \n\nThe theoretical analysis seems quite rigorous as it begins by first providing a basic guarantee for constrained langevin sampling when gradients are computed with error. The non error gradient part of this analysis has been established before and the authors mention the references appropriately. I have a couple of questions regarding the precise statements of the theorem that i will ask towards the end of the review. Overall its hard to comment on the tightness of the analysis as the non-error versions are also unclear of the tightness of the bounds. Nevertheless the bounds achieved do not look much worse than the non-error counterparts and are easy to implement. The rest of the bounds focus on achieving low error in approximation of the gradients in various settings. Overall the theory in this part seems very loose in terms of bounds as exponential factors in dimensions start to appear and in that regard seems quite preliminary but its hard to comment on whether its natural or can be improved. \n\nOverall the paper is a rigorous treatment of multiple components that would go into the problem of sampling zero order constrained sets and merges many ideas which can all be useful. Nevertheless the paper is a little lacking of novelty in the sense that it brings together many existing ideas and provides an analysis of the effect of bringing them together but none of the theory significantly improves over the existing theory.\n\nSome questions I have regarding the theorem statements \n\nRegarding theorem 1 SPLMC convergence  (and corollaries of the theorem) -  I find it surprising that there is no lower bound assumption on eta in terms of K - only an upper bound. This seems wrong particularly as the theorem as stated then allows eta to be set extremely small while K is finite, in which case there should be no convergence theorems at all. The condition on eta should be a theta(f (K)) for some f type of condition like in the second part of the theorem. I would suggest the authors to look into the theorem - or claify why this is the case. \n\nI am confused by the presentation of the Shi et al results as there is no penalty appearing for the approximation error due to an RKHS, only a finite sample penalty. Does the result assume that phi belongs to the function class of the RKHS in question? Probably yes and in that case that should be specified. \n", "extra_scores": {"experience_assessment": "I have read many papers in this area.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper thoroughly.", "review_assessment:_checking_correctness_of_experiments": "I did not assess the experiments.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I assessed the sensibility of the derivations and theory."}, "extra_text": {}}, {"review_id": "H1xFSl8rFr", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "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": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This paper proposes a solution to overcome the challenges due to the black-box nature of physical constraints that are involved in the design of nano-porous templates with optimal thermoelectric efficiency. \n\nUnfortunately, I cannot comment on the overall scientific contribution of the paper, as I do not possess the expertise to judge it accurately. My expertise is so outside of this field that I will rely on the judgement of the other reviewers, whom I hope will have more experience and will better know the literature. \n\nI can only report that the proposed method does not seem to be a particularly good approximation to the zeroth-order method since the mean values for kappa and sigma in table 2 are quite a bit worse than those obtained with the baseline. Of course, the proposed approach is quite a bit faster. However, the paper does not provide a sense of whether these values are actually useful. In practice, would one want to wait longer to get a better quality result, or are the numbers obtained with the proposed approach usable? \n\nAlso, could the proposed approach be applied to other problems? It would be great to see at least one or two other areas where this could be applied to, since I doubt the general ICLR audience is well-versed in nano-porous templates. ", "extra_scores": {"experience_assessment": "I do not know much about this area.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I did not assess the derivations or theory.", "review_assessment:_checking_correctness_of_experiments": "I did not assess the experiments.", "review_assessment:_thoroughness_in_paper_reading": "I made a quick assessment of this paper."}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper is not overly well written and motivated. A guiding thread through the paper is often missing. Comparisons with constrained BO methods would have improved the paper as well as a more explicit link to multi-objective BO. It could have been interesting to evaluate the sensitivity w.r.t. the number of samples in the Monte Carlo estimate. What happens if the observations of the function are noisy? Is there a natural way to deal with this?
Given that the paper is 10+ pages long, we expect a higher quality than an 8-pages paper (reviewing and submission guidelines). 
Decision: Reject