Bayesian Optimization for Materials Design

Abstract

We introduce Bayesian optimization , a technique developed for optimizing time-consuming engineering simulations and for fitting machine learning models on large datasets. Bayesian optimization guides the choice of experiments during materials design and discovery to find good material designs in as few experiments as possible. We focus on the case when materials designs are parameterized by a low-dimensional vector. Bayesian optimization is built on a statistical technique called Gaussian process regression , which allows predicting the performance of a new design based on previously tested designs. After providing a detailed introduction to Gaussian process regression, we describe two Bayesian optimization methods: expected improvement, for design problems with noise-free evaluations; and the knowledge-gradient method, which generalizes expected improvement and may be used in design problems with noisy evaluations. Both methods are derived using a value-of-information analysis, and enjoy one-step Bayes-optimality .

Derivations and Proofs

This section contains derivations and proofs of equations and theoretical results found in the main text.

1.1 Proof of Proposition 1

Proof

Using Bayes’ rule, the conditional probability density of $\theta _{[2]}$ at a point $u_{[2]}$ given that $\theta _{[1]} = u_{[1]}$ is

$$\begin{aligned} p(\theta _{[2]} = u_{[2]}&\mid \theta _{[1]} = u_{[1]}) = \frac{p(\theta _{[1]} = u_{[1]}, \theta _{[2]}=u_{[2]})}{p(\theta _{[1]} = u_{[1]})} \propto p(\theta _{[1]} = u_{[1]}, \theta _{[2]}=u_{[2]}) \nonumber \\&\propto \exp \left( -\frac{1}{2} \begin{bmatrix} u_{[1]} - \mu _{[1]} \\ u_{[2]} - \mu _{[2]}\end{bmatrix}^T \begin{bmatrix}\varSigma _{[1,1]}&\varSigma _{[1,2]} \\ \varSigma _{[2,1]}&\varSigma _{[2,2]}\end{bmatrix}^{-1} \begin{bmatrix} u_{[1]} - \mu _{[1]} \\ u_{[2]} - \mu _{[2]}\end{bmatrix} \right) . \end{aligned}$$

(3.19)

To deal with the inverse matrix in this expression, we use the following identity for inverting a block matrix: the inverse of the block matrix $\begin{bmatrix} A&B \\ C&D\end{bmatrix}$, where both A and D are invertible square matrices, is

$$\begin{aligned} \begin{bmatrix} A&B \\ C&D\end{bmatrix}^{-1} = \begin{bmatrix} (A - BD^{-1}C)^{-1}&-(A-BD^{-1}C)^{-1}BD^{-1} \\ -(D-CA^{-1}B)^{-1}CA^{-1}&(D-CA^{-1}B)^{-1}\end{bmatrix}. \end{aligned}$$

(3.20)

Applying (3.20) to (3.19), and using a bit of algebraic manipulation to get rid of constants, we have

$$\begin{aligned} p(\theta _{[2]}=u_{[2]} \mid \theta _{[1]} = u_{[1]}) \propto \exp \left( -\frac{1}{2}(u_{[2]} - \mu ^{\text {new}})^T (\varSigma ^{\text {new}})^{-1} (u_{[2]} - \mu ^{\text {new}}) \right) , \end{aligned}$$

(3.21)

where $\mu ^{\text {new}} = \mu _{[2]} - \varSigma _{[2,1]}\varSigma _{[1,1]}^{-1}(u_{[1]} - \mu _{[1]})$ and $\varSigma ^{\text {new}} = \varSigma _{[2,2]}-\varSigma _{[2,1]}\varSigma _{[1,1]}^{-1}\varSigma _{[1,2]}$.

We see that (3.21) is simply the unnormalized probability density function of a normal distribution. Thus the conditional distribution of $\theta _{[2]}$ given $\theta _{[1]} = u_{[1]}$ is multivariate normal, with mean $\mu ^{\text {new}}$ and covariance matrix $\varSigma ^{\text {new}}$.

1.2 Derivation of Equation (3.16)

Since $f(x) \sim \text {Normal}(\mu _n(x), \sigma _n^2(x))$, the probability density of f(x) is $p(f(x)=$ $z) = \frac{1}{\sqrt{2\pi }} \exp \left( (z - \mu _n(x))^2/2 \sigma _n(x)^2\right)$. We use this to calculate $\mathrm {EI}(x)$:

$$\begin{aligned} \begin{aligned} \mathrm {EI}(x)&= E_n[ (f(x)-f^*_n)^+] \\&= \int _{f^*_n}^{\infty } (z-f^*_n) \frac{1}{\sqrt{2\pi }\sigma _n(x)} e^{\frac{-(z - \mu _n(x))^2}{2 \sigma _n^2(x)}} dz \\&= \int _{f^*_n}^{\infty } z \frac{1}{\sqrt{2\pi }\sigma _n(x)} e^{\frac{-(z - \mu _n(x))^2}{2 \sigma _n^2(x)}} dz - f_n^*\left( 1-\varPhi \left( \frac{f_n^*-\mu _n(x)}{\sigma _n(x)}\right) \right) \\&= \int _{f^*_n}^{\infty } \left( \mu _n(x) + (z-\mu _n(x)) \right) \frac{1}{\sqrt{2\pi }\sigma _n(x)} e^{\frac{-(z - \mu _n(x))^2}{2 \sigma _n^2(x)}} dz - f_n^*\left( 1-\varPhi \left( \frac{f_n^*-\mu _n(x)}{\sigma _n(x)}\right) \right) \\&= \int _{f^*_n}^{\infty } \left( z-\mu _n(x) \right) \frac{1}{\sqrt{2\pi }\sigma _n(x)} e^{\frac{-(z - \mu _n(x))^2}{2 \sigma _n^2(x)}} dz + (\mu _n(x)-f_n^*)\left( 1-\varPhi \left( \frac{f_n^*-\mu _n(x)}{\sigma _n(x)}\right) \right) \\&= \sigma _n(x)\frac{1}{\sqrt{2\pi }} e^{\frac{-(f^*_n - \mu _n(x))^2}{2 \sigma _n(x)^2}} + (\mu _n(x)-f_n^*)\left( 1-\varPhi \left( \frac{f_n^*-\mu _n(x)}{\sigma _n(x)}\right) \right) \\&= (\mu _n(x)-f_n^*)\left( 1-\varPhi \left( \frac{f_n^*-\mu _n(x)}{\sigma _n(x)}\right) \right) + \sigma _n(x) \varphi \left( \frac{f_n^*-\mu _n(x)}{\sigma _n(x)}\right) \\&= (\mu _n(x)-f_n^*)\varPhi \left( \frac{\mu _n(x) - f_n^*}{\sigma _n(x)}\right) + \sigma _n(x) \varphi \left( \frac{\mu _n(x) - f_n^*}{\sigma _n(x)}\right) . \end{aligned} \end{aligned}$$

1.3 Calculation of the KG factor

The KG factor (3.18) is calculated by first considering how the quantity $\mu ^*_{n+1} - \mu ^*_n$ depends on the information that we have at time n, and the additional datapoint that we will obtain, $y_{n+1}$.

First observe that $\mu ^*_{n+1} - \mu ^*_n$ is a deterministic function of the vector $[\mu _{n+1}(x) : x\in A_{n+1}]$ and other quantities that are known at time n. Then, by applying the analysis in Sect. 3.3.5, but letting the posterior given $x_{1:n},y_{1:n}$ play the role of the prior, we obtain the following version of (3.10), which applies to any given x,

$$\begin{aligned} \mu _{n+1}(x) = \mu _n(x) + \frac{\varSigma _n(x,x_{n+1})}{\varSigma _n(x_{n+1},x_{n+1}) + \lambda ^2}\left( y_{n+1} - \mu _n(x_{n+1})\right) . \end{aligned}$$

(3.22)

In this expression, $\mu _n(\cdot )$ and $\varSigma _n(\cdot ,\cdot )$ are given by (3.13) and (3.14).

We see from this expression that $\mu _{n+1}(x)$ is a linear function of $y_{n+1}$, with an intercept and a slope that can be computed based on what we know after the nth measurement.

We will calculate the distribution of $y_{n+1}$, given what we have observed at time n. First, $f(x_{n+1}) | x_{1:n},y_{1:n} \sim \text {Normal}\left( \mu _n(x_{n+1}), \varSigma _n(x_{n+1}, x_{n+1})\right)$. Since $y_{n+1} = f(x_{n+1}) + \varepsilon _{n+1}$, where $\varepsilon _{n+1}$ is independent with distribution $\varepsilon _{n+1} \sim \text {Normal}(0, \lambda ^2)$, we have

$$\begin{aligned} y_{n+1} | x_{1:n}, y_{1:n} \sim \text {Normal}\left( \mu _n(x_{n+1}), \varSigma _n(x_{n+1}, x_{n+1}) + \lambda ^2 \right) . \end{aligned}$$

Plugging the distribution of $y_{n+1}$ into (3.22) and doing some algebra, we have

$$\begin{aligned} \mu _{n+1}(x) | x_{1:n}, y_{1:n} \sim \text {Normal} \left( \mu _n(x), \widetilde{\sigma }^2(x, x_{n+1}) \right) , \end{aligned}$$

where $\widetilde{\sigma }(x, x_{n+1}) = \frac{\varSigma _n(x, x_{n+1})}{\sqrt{\varSigma _n(x_{n+1}, x_{n+1}) + \lambda ^2}}$. Moreover, we can write $\mu _{n+1}(x)$ as

$$\begin{aligned} \mu _{n+1}(x) = \mu _n(x) + \widetilde{\sigma }(x, x_{n+1}) Z, \end{aligned}$$

where $Z=(y_{n+1}-\mu _n(x_{n+1}))/ \sqrt{\varSigma _n(x_{n+1}, x_{n+1}) + \lambda ^2}$ is a standard normal random variable, given $x_{1:n}$ and $y_{1:n}$.

Now (3.18) becomes

$$\begin{aligned} \begin{aligned} \mathrm {KG}_n(x)&= E_n \left[ \max _{x'\in A_{n+1}} \mu _n(x') + \widetilde{\sigma }(x', x_{n+1}) Z \mid x_{n+1}=x \right] - \mu ^*_n. \end{aligned} \end{aligned}$$

Thus, computing the KG factor comes down to being able to compute the expectation of the maximum of a collection of linear functions of a scalar normal random variable. Algorithm 2 of [24], with software provided as part of the matlabKG library [53], computes the quantity

$$\begin{aligned} h(a,b) = \mathbb {E}\left[ \max _{i = 1, \ldots , |a|} (a_i + b_i Z)\right] - \max _{i = 1, \ldots , |a|} a_i \end{aligned}$$

for arbitrary equal-length vectors a and b. Using this ability, and letting $\mu _n(A_{n+1})$ be the vector $[\mu _n(x'): x'\in A_{n+1}]$ and $\widetilde{\sigma }(A_{n+1},x)$ be the vector $[\widetilde{\sigma }(x', x): x'\in A_{n+1}]$, we can write the KG factor as

$$\begin{aligned} \mathrm {KG}_n(x) = h(\mu _n(A_{n+1}),\widetilde{\sigma }(A_{n+1},x)) + \left[ \max (\mu _n(A_{n+1})) - \mu ^*_n\right] . \end{aligned}$$

If $A_{n+1}=A_n$, as it is in the versions of the knowledge-gradient method proposed in [24, 25], then the last term $\max (\mu _n(A_{n+1})) - \mu ^*_n$ is equal to 0 and vanishes.