Bayesian Optimization for Materials Design with Mixed Quantitative and Qualitative Variables

In plain words

Designing a material often means choosing discrete options, such as which ingredient or which surface treatment to use, and also setting continuous knobs, such as a thickness or a ratio. Each simulation or experiment is expensive, so you want to find a good design in as few runs as possible. Standard optimizers that model the problem with a smooth curve handle the continuous knobs well but struggle with discrete options, because there is no natural distance between "option A" and "option B".

The usual workaround is to turn each option into a yes/no switch and let the model treat every switch as a separate, unrelated problem. This paper instead gives each option a learned position on a flat 2D map. Options that behave alike end up close together, and options that behave differently end up far apart. The model uses these positions to make its predictions, and an optimizer then decides which design to simulate next by balancing "looks promising" against "still uncertain". The authors show that the learned map can also be read to see which options are similar.

TL;DR

• Standard Gaussian-process Bayesian optimization needs numeric inputs; the common fix encodes each qualitative level as dummy variables, which the paper argues is theoretically restrictive p. 6.
• LVGP maps every qualitative level to a 2D latent point learned by maximum likelihood, and uses those points inside a standard GP correlation p. 6, eq. 3.
• Expected improvement is the acquisition function, with the loop given in tab. 3.
• Author-reported: LV-EI beats the MATLAB bayesopt baseline (dummy variables) on two test functions and two materials problems p. 7, p. 13, p. 15–16.
• No code is released ("upon request") p. 29.

Problem

Materials and device design mix quantitative variables (thickness, ratio, microstructure statistics) with qualitative ones (constituent type, surface treatment, solvent) p. 2. A GP with 0/1 dummy inputs for each qualitative level amounts to a separate response surface per combination of levels. The authors argue that this is restrictive when the number of levels is large and cannot represent correlations between levels p. 5–6. Each evaluation is an expensive simulation, so the method must work with few samples p. 3–4.

The method assumes that the effect of any qualitative level on the response is produced by some underlying quantitative physical variables. A low-dimensional latent representation (here 2D) is used as a stand-in for them p. 6, p. 19–20.

Method

Standard GP. The response is modelled as a constant mean plus a zero-mean GP with a stationary covariance, using a Gaussian correlation over the quantitative inputs eq. 1, eq. 2. The correlation parameters φ are estimated by maximum likelihood p. 18.

LVGP correlation. Each qualitative factor t_j with m_j levels is assigned a 2D latent point z^j(t_j) for each of its levels p. 18–19. The correlation between two points adds a squared Euclidean distance between their latent points to the usual quantitative term eq. 3. The latent points are unknown and are estimated together with μ, σ² and φ p. 20. The paper uses a Gaussian correlation here; it states that LVGP also works with power-exponential, Matérn and lifted Brownian correlations for the quantitative part p. 6, p. 18.

# LVGP correlation between two design points  [eq. 3]
# x: quantitative inputs (p,); t: qualitative levels (q,)
# phi: weights for quantitative inputs (p,)
# Z[j][level]: 2D latent point for level of factor j  (2,)

def lvgp_correlation(x, t, x_prime, t_prime, phi, Z):
    # squared distance over the quantitative inputs  [eq. 3]
    quant_term = sum(phi[i] * (x[i] - x_prime[i]) ** 2 for i in range(p))
    # squared distance between latent points, per qualitative factor  [eq. 3]
    qual_term = 0.0
    for j in range(q):
        qual_term += squared_norm(Z[j][t[j]] - Z[j][t_prime[j]])
    return exp(-quant_term - qual_term)

Likelihood and fitting. Given n observations, the correlation matrix R is built from eq. 3 for all pairs. μ and σ² are profiled out in closed form eq. 5, which leaves a function of φ and the latent points Z only eq. 6. If the simulation is noisy, a nugget λ is added to the diagonal of R and estimated together with φ and Z p. 21.

# Negative profile log-likelihood for LVGP  [eq. 5], [eq. 6]
# X: (n, p); T: (n, q) level indices; y: (n,)

def neg_profile_loglik(phi, Z, nugget):
    # correlation matrix over all training points  [eq. 3]
    R = correlation_matrix(X, T, phi, Z) + nugget * identity(n)
    # closed-form mean and variance  [eq. 5]
    mu = (ones @ inv(R) @ y) / (ones @ inv(R) @ ones)
    sigma2 = (y - mu) @ inv(R) @ (y - mu) / n
    # profile log-likelihood, negated for minimisation  [eq. 6]
    return n * log(sigma2) + log(det(R))

The optimiser, its starting values, bounds on φ, and the treatment of the latent points are not stated p. 21.

Prediction. At a new point x*, the mean is the constant mean plus the correlation vector times R⁻¹ applied to the centred data eq. 7. The predictive variance is given by eq. 8 (see Pitfalls for a typo in this equation).

Acquisition. Expected improvement (EI) is used for the deterministic examples and for the noisy solar-cell example p. 7, p. 10. With y_min the lowest value observed so far, Δ(x) = y_min − ŷ(x), and φ and Φ the standard normal density and distribution eq. 10:

EI(\mathbf{x}) = \hat s(\mathbf{x})\,\phi\!\left(\frac{\Delta(\mathbf{x})}{\hat s(\mathbf{x})}\right) + \Delta(\mathbf{x})\,\Phi\!\left(\frac{\Delta(\mathbf{x})}{\hat s(\mathbf{x})}\right) eq. 10

# Expected improvement for minimisation  [eq. 10]
# mean, std: predictive mean and standard deviation at x
# y_min: lowest observed value so far

def expected_improvement(mean, std, y_min):
    # improvement relative to the best value so far  [eq. 10]
    delta = y_min - mean
    z = delta / std
    return std * pdf(z) + delta * cdf(z)

For noisy responses, the paper states that EI with a "plug-in" value μ_min(x) is used, and says this is introduced in the Methods section p. 10. The Methods section does not define it (see Pitfalls).

Optimisation loop. Table 3 gives the loop: build an initial set, refit the LVGP, pick the design maximising EI, evaluate it, and add it to the data tab. 3, p. 24. The paper describes the same four steps in words (create data, fit, choose the next point, evaluate and augment) p. 22.

# Bayesian optimisation loop with LVGP  [tab. 3]
# data: list of (x, t, y) observations; y is minimised

def lvgp_bo(initial_data, budget):
    # (0) start from the initial design set  [tab. 3]
    data = initial_data
    # (1) repeat for each iteration  [tab. 3]
    for n in range(budget):
        # (2) refit the latent-variable GP on all data so far  [tab. 3]
        model = fit_lvgp(data)                    # [eq. 4]-[eq. 6]
        # (3) choose the design with the largest expected improvement  [tab. 3]
        next_x = argmax_over_design(model, expected_improvement)
        # (4) run the simulation or experiment  [tab. 3]
        next_y = evaluate(next_x)
        # (5) augment the data set  [tab. 3]
        data = data + [(next_x, next_y)]
    return data

The inner maximisation (argmax_over_design) is not described for either continuous or discrete design spaces p. 24.

Hyperparameters

Table 1: Refractive index for the candidate materials p. 8

Table 2: Range of design variables in the solar-cell problem p. 9

Results

Author-reported. Each row names the setup and the baseline. The paper gives no numeric values for the boxplot insets, so they are cited to the figures.

Other claims in the paper: LVGP gives better RMSE than existing GP models with qualitative inputs across several examples. This is cited to the authors' earlier study and is not re-tested here p. 6.

Paper vs code

No official code. The paper says that source code will be made available "upon request to the corresponding author" p. 29. The only repository it links is clancyLab/NCM2018, which belongs to the authors of ref. 35 (the HOIP data source), not to this paper's authors, and is not used here.

GitHub searches for the method name "LVGP" and for "latent variable Gaussian process Bayesian optimization" returned no repository from the authors. A search for lvgpytorch also returned nothing. Other repositories named LVGP belong to unrelated owners and were not used.

Pitfalls

• No code. Any implementation must be built from eq. 3–eq. 10 and tab. 3. Optimiser, starting values, bounds on φ and the latent points, and the design-space search are all not stated p. 21, p. 24.
• Equation 8 as printed. The first term reads r(x*), which is a vector, so the expression is dimensionally inconsistent. The standard kriging variance has a scalar constant in that place. The paper does not say which form it intends; check against the GP reference it cites (ref. 37) before implementing p. 21.
• Missing equation 9. The numbering jumps from eq. 8 to eq. 10 p. 21, p. 24. Do not assume a missing formula; the EI formula is only given as eq. 10.
• Undefined noisy EI. The solar-cell section says the "plug-in" μ_min(x) is introduced in the Methods for noisy responses p. 10, but the Methods section does not define it. Only a reference to a benchmark study is given p. 22.
• Latent positions are identifiable only up to rigid motions. Eq. (3) depends only on squared distances between latent points, so the latent plots are defined only up to rotation, reflection and translation. The paper does not state this; it follows from eq. 3. Compare latent maps by distance, not coordinates.
• Latent-point weights are unstated. Eq. (3) gives no separate scale for the latent distance relative to the quantitative terms, so the relative weighting is fixed by the form of the model and not by a stated parameter eq. 3.
• Noise handling. The solar-cell response is averaged over three statistically equivalent microstructures per design point p. 9–10. The paper does not quantify the residual noise, and the nugget is only mentioned for non-deterministic responses p. 21.
• Sign convention in the HOIP example. The text uses "binding energy" and "−binding energy" interchangeably p. 12–14, and the figure axes show the negated quantity fig. 6a, fig. 6c. Check which sign an implementation maximises or minimises.
• HOIP initial set is filtered. The starting samples are restricted to designs with binding energy above −30, which the authors describe as making the search harder p. 13. Results will depend on this filter.
• Level counts. The HOIP problem has five qualitative variables with 3, 3, 3, 3 and 8 levels, giving 3 × 3 × 3 × 3 × 8 = 648 combinations. The paper says 240 are stable p. 12. Only stable combinations are used.
• Cross-reference errors in the paper. The refractive-index table is Table 1, but the text cites "Table 2" p. 10. The framework figure is Figure 2, but the text cites "Figure 1" p. 22. The text gives a 2000 nm unit cell p. 8, while the Figure 4 scale bar reads 1000 nm fig. 4.
• Duplicate reference. Refs. 27 and 28 are the same paper p. 27–28.
• Baseline settings. The MATLAB bayesopt configuration (kernel, fitting, initial design) is not stated beyond the use of EI and the dummy-variable model p. 7.
• Cost. The paper reports no run times. The loop refits the full model at every iteration tab. 3.
• Generality. The authors state the method is generic for mixed-variable problems with expensive simulations p. 1, but all tests are small: two test functions and two materials problems. Results on larger qualitative spaces are not shown.

Lineage

• Builds on: the LVGP model of Zhang et al. (Technometrics, 2019), cited as refs. 27–28 of this paper. Not in the wiki; no arXiv id given p. 6.
• Compared against: MC-EI, the MATLAB bayesopt baseline with dummy variables, equivalent to a multiplicative covariance model (refs. 24, 26) p. 7.
• Superseded by: not stated.