Exploration of new chemical materials using black-box optimization with the D-wave quantum annealer

In materials informatics, searching for chemical materials with desired properties is challenging due to the vastness of the chemical space. Moreover, the high cost of evaluating properties necessitates a search with a few clues. In practice, there is also a demand for proposing compositions that are easily synthesizable. In the real world, such as in the exploration of chemical materials, it is common to encounter problems targeting black-box objective functions where formalizing the objective function in explicit form is challenging, and the evaluation cost is high. In recent research, a Bayesian optimization method has been proposed to formulate the quadratic unconstrained binary optimization (QUBO) problem as a surrogate model for black-box objective functions with discrete variables. Regarding this method, studies have been conducted using the D-Wave quantum annealer to optimize the acquisition function, which is based on the surrogate model and determines the next exploration point for the black-box objective function. In this paper, we address optimizing a black-box objective function containing discrete variables in the context of actual chemical material exploration. In this optimization problem, we demonstrate results obtaining parameters of the acquisition function by sampling from a probability distribution with variance can explore the solution space more extensively than in the case of no variance. As a result, we found combinations of substituents in compositions with the desired properties, which could only be discovered when we set an appropriate variance.


INTRODUCTION
Black-box optimization is a method to optimize a function that does not have an explicit objective function in the mathematical form.In the real world, this optimization problem appears in various fields, including material informatics, robotics (Deisenroth, 2011), machine learning (Snoek et al., 2012), and recommendation systems (Vanchinathan et al., 2014).Bayesian optimization is one of the solutions for black-box optimization problems (Jones et al., 1998).Taking the exploration of chemical materials as an example, a surrogate model is constructed using an existing dataset to predict the relationship between the combinations of substituents in the chemical materials and the corresponding property values.Based on this surrogate model, an acquisition function is defined.The combination of substituents obtained through optimizing this acquisition function is then used as the next input point for the black-box objective function, enabling the evaluation of the actual property values.The relationship between the inputted combination of substituents and the actual property value is then added to the existing dataset, then updating the surrogate model.Repeating this process is to explore the combinations of substituents that yield the desired property values.Especially for black-box optimization problems involving discrete variables, discrete variables are included in both the surrogate model and the acquisition function.Therefore, even optimizing the acquisition function often proves to be NP-hard, and the solutions obtained through optimization are generally approximate.In a previous study, Bayesian optimization of combinatorial structures (BOCS) was proposed as the promising algorithm for such problems (Baptista and Poloczek, 2018).In this algorithm, the acquisition function was assumed as quadratic unconstrained binary optimization (QUBO) problem.
Quantum annealing (Kadowaki and Nishimori, 1998) is a heuristic algorithm to solve QUBO problems by driving binary variables through quantum fluctuations.Many well-known combinatorial optimization problems can be encoded into QUBO problems (Lucas, 2014).Practical applications of quantum annealing can be found in various fields, including traffic flow optimization (Neukart et al., 2017;Inoue et al., 2021;Shikanai et al., 2023), manufacturing (Ohzeki et al., 2019;Haba et al., 2022), finance (Rosenberg et al., 2015;Venturelli and Kondratyev, 2019), steel manufacturing (Yonaga et al., 2022), decoding problems (Ide et al., 2020;Arai et al., 2021), and algorithms in machine learning (Amin et al., 2018;O'Malley et al., 2018;Urushibata et al., 2022;Hasegawa et al., 2023;Goto and Ohzeki, 2023).Furthermore, quantum annealing, which utilizes the quantum tunneling effect, is expected to find the optimal solution for several combinatorial optimization problems more rapidly than algorithms such as simulated annealing (Kirkpatrick et al., 1983).This advantage is investigated from the perspective of energy landscape characteristics (Das and Chakrabarti, 2008) and through numerical computation (Denchev et al., 2016).In addition, there are discussions about the characteristics of solutions obtained in cases where multiple optimal solutions exist (Yamamoto et al., 2020;Maruyama et al., 2021).With these backgrounds, quantum annealing has recently attracted attention, both for its potential applications and for validating the fundamental aspects of quantum effects.
Studies that employ quantum annealing in some algorithms for black-box optimization problems involving discrete variables exist.These include benchmark tests (Koshikawa et al., 2021) that have examined the presence or absence of quantum superiority in optimizing acquisition functions.In terms of practical applications, there are case studies that have achieved significant screening in the exploration of chemical materials within the search chemical space (Hatakeyama-Sato et al., 2021;Takuro Tanaka et al., 2023), as well as instances of designing complex metamaterials (Kitai et al., 2020).
In the exploration of chemical materials, it is necessary not only to discover molecules with the desired property values but also to be concerned about scenarios in actual synthesis where molecules with specific substructures may become entirely unfeasible to synthesize.Drawing inspiration from previous studies and practical needs, we demonstrate a method for proposing diverse compositions of chemical materials with desired properties, targeting a black-box optimization problem that includes discrete variables in actual chemical material exploration.In more detail, we show results that by obtaining parameters of the surrogate model and acquisition function from sampling a probability distribution with an appropriate variance and optimizing the acquisition function, we explored the solution space more extensively while optimizing the black-box objective function.The method used in this paper is generally referred to as Thompson sampling (THOMPSON, 1933;Chapelle and Li, 2011).In this sense, it can be said that our research results evaluate the impact of the magnitude of the variance of the posterior probability distribution in Thompson sampling.
The remaining sections of this paper are organized as follows: In the next section, Section 2, we explain the problem setting in this paper and the method we propose.In Section 3, we demonstrate the results of the experiments related to the actual exploration of chemical materials.Finally, Section 4 summarizes our research and discusses this paper and future research directions.

MATERIALS AND METHODS
In this section, we introduce the problem settings based on the search for chemical materials, which is the focus of this paper.Subsequently, in Bayesian optimization, we explain the construction of the surrogate model in the QUBO form, which is well-known in prior research, along with the construction of the acquisition function.We provide this explanation in conjunction with our method aim.

Problem settings
In this paper, we define the binding of substituents to specific sites of the molecular frame as the composition of chemical materials.We aim to propose various combinations of substituents through Bayesian optimization while maximizing a target material property value.To align our description with other literature focusing on black-box optimization problems, we define our goal as a minimization problem, utilizing the fact that maximization and minimization problems can be transformed into each other by reversing the sign of the objective function.

Methods
We express the assignment of substituents using a binary vector.In particular, for substituents that can bind to each site, we encode them by converting the 0-indexed substituent number to binary.Thus, we set a binary vector ⃗ x (µ) ∈ {0, 1} N as input, and the corresponding target material property value y (µ) as output.We aim to find ⃗ x that minimizes a black-box objective function.Since we cannot know an explicit form of the black-box objective function, we construct a surrogate model as QUBO form following the previous studies.We utilize an existing dataset D = {⃗ x (µ) , y (µ) } D µ=1 to sample the parameters of the surrogate model from a probability distribution we discuss later and construct it.Based on the surrogate model, we construct an acquisition function and propose a combination of substituents that optimize the acquisition function using the D-Wave quantum annealer.Subsequently, we input the proposed combination of substituents as the next exploration point ⃗ x (new) and obtain output y (new) from the black-box objective function.Then we append {⃗ x (new) , y (new) } to the existing dataset as new data and reconstruct the surrogate model.By repeating this process, we aim to obtain diverse combinations of substituents with desired target material property values.

Construction of surrogate model function
We construct the surrogate model f surrogate (⃗ x) in the QUBO form in this paper.
For simplicity, we set the surrogate model parameters {α i , α ij } = ⃗ α ∈ R p .Note that p = 1 + N + N (N − 1)/2.Defining X ∈ {0, 1} D×p as the design matrix and denoting the µ-th row in the design matrix X as X (µ) , we have the following expression N .Furthermore, we set the output vector ⃗ y ∈ R D and I as the identity matrix.Then, we assume a prior distribution of surrogate model parameters P (⃗ α) with a variance σ 2 α I and a likelihood function over the surrogate model parameters ⃗ α with a variance σ 2 y I.We give the prior distribution and likelihood function as following multivariate Gaussian distributions.
At this time, the posterior distribution of the surrogate model parameters ⃗ α is computed and given by a multivariate Gaussian distribution, similar to the prior distribution and the likelihood function.
We sample the surrogate model parameters ⃗ α ∈ R p from the multivariate Gaussian distribution described in (4).σ 2 is a hyperparameter indicating the magnitude of fluctuations from the mean vector ⃗ µ when sampling the surrogate model parameters.λ is also a hyperparameter.Note that λ corresponds to the coefficient of the regularization term during ridge regression.

Construction of acquisition function
The acquisition function f acquisition (⃗ x) is constructed in the same QUBO form as the surrogate model, and the next exploration point ⃗ x (new) is proposed by optimizing the acquisition function.
) is a function with modified specific parameters from the surrogate model f surrogate (⃗ x) described in 2.2.1.This modification is like a penalty method, designed to ensure that binary vectors with substituent numbers that do not exist at each site do not become the optimal points of the acquisition function.Parameters that are not modified are identical to those in the surrogate model f surrogate (⃗ x).For example, when six potential substituents can bind at a specific site, representing the 0-indexed substituent numbers in binary requires three bits (x 1 , x 2 , x 3 ).In this context, x 1 x 2 x 3 = (000, 001, 010, 011, 100, 101) 2 corresponds to valid substituent numbers from 0 to 5 in decimal.However, each combination x 1 x 2 x 3 = (110, 111) 2 is equivalent to substituent numbers 6-7 in decimal, rendering them inappropriate as optimal point candidates.To prevent the substituent combinations with substituent numbers 6-7 at this site from being proposed as the optimal points of the acquisition function, we adjust the surrogate model parameters.
In this example, we modify the coefficient of x 1 x 2 in the surrogate model function to a positive constant C, and the other coefficients are kept the same as in the surrogate model.The next exploration point of the black-box objective function is determined by the optimization of the acquisition function f acquisition (⃗ x).
The search space explored varies greatly depending on how the acquisition function is constructed and how the acquisition function is optimized.As described, our method samples the parameters of the surrogate model and the acquisition function from a probability distribution with variance.The hyperparameter σ 2 indicates the magnitude of the variance.The larger this hyperparameter σ 2 is, the more significant the variance of the acquisition function, potentially allowing for exploration across a broader solution space and avoiding resampling the previously explored points.

RESULTS
In this section, we describe detailed problem settings and experimental conditions and then show the experimental results obtained by applying our method.In particular, we compare and discuss based on the magnitude of the hyperparameter σ 2 .Our discussion centers on two main points of interest in this paper.The first point is whether our method has brought diversity to the proposed substituent combinations.The second point is whether our method has optimized the black-box objective function.

Detailed problem settings and experimental conditions
We set the number of substituent binding sites as four, and for convenience in the description, we call each binding site R1, R2, R3, and R4, respectively.The number of possible substituents that can bind at each site is R1: 6, R2: 29, R3: 64, and R4: 64, respectively.Therefore, the size of the chemical space is calculated as 6 × 29 × 64 × 64 = 712704.Moreover, the number of bits necessary to represent the number of each substituent is R1: 3, R2: 5, R3: 6, and R4: 6.Consequently, the binary vector ⃗ x dimension is calculated as N = 3 + 5 + 6 + 6 = 20.The substituent number at R1 is represented in 0-indexed form using x 1 to x 3 , similarly, x 4 to x 8 represent the substituent number at R2, x 9 to x 14 represent the substituent number at R3, and x 15 to x 20 represent the substituent number at R4.To illustrate with a concrete encoding example, suppose the substituent numbers at each site are R1: 0, R2: 2, R3: 10, and R4: 63.In this case, the binary vector ⃗ x would be represented as ⃗ x = (0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1).We set the hyperparameter λ at 10 −2 and the hyperparameter σ 2 , which indicates the magnitude of fluctuation from the mean vector ⃗ µ when sampling surrogate model parameters, to {0, 4 × 10 −3 , 8 × 10 −3 , 12 × 10 −3 }.We set the surrogate model's parameter correction for R1 in the acquisition function as C = α 12 = 2 × max(⃗ α) at each after sampling ⃗ α.We used D-Wave Advantage 4.1 as the quantum annealer, setting the annealing time to 2000µs, and the number of samples is 300.The quantum adiabatic theorem ensures that it is possible to find the nontrivial ground state at the end of the quantum annealing if the transverse field changes sufficiently slowly (Suzuki and Okada, 2005;Morita and Nishimori, 2008;Ohzeki and Nishimori, 2011).On the other hand, when quantum annealing is carried out on a physical device D-wave quantum annealer, it operates at a finite temperature and is subject to external noise.Due to these factors, the annealing time is often short in many studies.Considering these theoretical and experimental backgrounds, we set the annealing time to be longer in our setting because we observed a tendency for the results to stabilize, possibly due to the effects of ambient temperature.The number of samples in the initial dataset is 992.For comparison as a baseline, we also conducted an experiment where the optimization part of the acquisition function was replaced with random sampling.Due to the nature of this study, which is conducted in the context of actual chemical material exploration, the computational cost of the black-box objective function is exceptionally high, resulting in an experiment of only one instance.We defined one loop as carrying out the following steps (i) through (v), and we performed 20 loops.
(i) By sampling the surrogate model parameters ⃗ α from a multivariate Gaussian distribution N (⃗ µ, Σ) described in (4), construct the surrogate model.(ii) Construct the acquisition function by partially correcting the surrogate model parameters as explained in 2.2.2.(iii) Optimize the acquisition function by quantum annealing and select the top 10 points of the acquisition function as the next exploration points for the black-box objective function.In the random sampling used as a baseline, 10 sampling points are randomly selected.Note that at this time, the top 10 points exclude combinations of substituents that are already present in the existing dataset and combinations of substituents that include non-existent substituent numbers, such as substituent numbers 6-7 in R1 and substituent numbers 29-31 in R2, through screening.(iv) Take the next exploration points obtained in (iii) as inputs and get outputs, carrying out the evaluation of target material property values, which is the computation of the black-box objective function, through DFT (Density Functional Theory) calculations.The detailed calculation method is described in the Additional Requirements.(v) Append the new samples {⃗ x (new) , y (new) } obtained in (iv) to the existing dataset and return to (i).

Histogram of substituent numbers in combinations added by end of the experiment
We show the histogram of substituent numbers at the binding sites R1, R2, R3, and R4 for the combinations of substituents added to the dataset by the end of the experiment in Figure 1.In the case of σ 2 = 0, we observed a tendency in R3 and R4 where specific substituent numbers were frequently proposed.However, as σ 2 increases, it can be observed that diversity is brought into the combinations of substituents proposed for R3 and R4.This difference is particularly pronounced when comparing σ 2 = 0 and σ 2 = 12 × 10 −3 .From these results, we can infer that we realized the proposal of various combinations of substituents by sampling parameters of the surrogate model and the acquisition function from probability distributions with variance.By sampling parameters from probability distributions with larger variances, the optimal points and the shape of the acquisition function change significantly in each loop.We believe that this approach allowed us to explore the solution space without getting trapped by some specific approximate solutions and without resampling the previously explored points.

Relationship between the number of loops and the R 2 of the surrogate model
We show the transition of the coefficient of determination R 2 in the surrogate model at each loop in Figure 2. The coefficient of determination R 2 is calculated from the initial dataset sample points, 992 points, and the sample points appended up to each loop.Note that R 2 , plotted in Figure 2, represents the results of mean-based regression.This result is equivalent to the regression of the maximum a posteriori (MAP) estimation.As σ 2 becomes larger, a tendency for R 2 at each loop to become smaller was observed.We speculate that we can attribute this result to the tendency shown in Figure 1, where the larger σ 2 is, the more diverse the combinations of substituents that the optimization of the acquisition function proposes become.When σ 2 is small, R 2 improves by fitting to similar input vectors and outputs.However, to improve R 2 when σ 2 is large, it is necessary to fit diverse input vectors and outputs.We speculate that this difficulty is why there was the tendency for the coefficient of determination, R, to be smaller when σ 2 is larger.

Analysis of target material property values
Finally, we show the target material property values evaluated by DFT calculations, corresponding to the combinations of substituents proposed through the optimization of the acquisition function as the next exploration point of the black-box objective function in Figure 3 and Figure 4.In Figure 3, we plot the transition of the best target material property values in the existing dataset up to each loop.Although we could only experiment once because of the extremely high computational cost of the black-box objective function, in the case of optimizing the acquisition functions, we confirm that it is possible to search for combinations of substituents with higher target material property values than the best value in the initial dataset.Under the conditions set in this study, using random sampling in the optimization part of the acquisition function, we could not find any combination of substituents that exhibited a property value To reiterate, the objective of black-box optimization in this study was to maximize the target material property value while bringing diversity to the combinations of substituents.Therefore, we listed the combinations of substituents whose target material property values exceeded our criteria of 0.880 or higher in Table 1,Table 2, Table 3 and Table 4. From the perspective of the number of combinations of substituents with property values that exceed our criteria, the number of proposed combinations was the highest at 25 combinations when σ 2 = 0.However, considering the diversity of proposed combinations of substituents, which is one of the aims of this paper, the advantage can be found when σ 2 ̸ = 0. Especially in the case of σ 2 = 4 × 10 −3 , it was possible to discover combinations of substituents with the property values that exceed our criteria, which have the substituent number of R4:0, a combination not discovered in case of σ 2 = 0.

DISCUSSION
In this study, we achieved the exploration of diverse approximate solutions in black-box optimization, which has the background of new chemical material discovery, by considering appropriate fluctuations in the parameters of the surrogate model and the acquisition function.Although the validity of the result is debatable because of the one-instance experiment, our result indicates that quantum annealing can accelerate the discovery of diverse chemical materials with desired material property values in materials informatics.More generally, our results demonstrate the advantages and disadvantages of varying the magnitude of the variance when sampling the parameters of the surrogate model from a probability distribution in optimizing a black-box objective function.In this paper, we explored a broader solution space by devising the construction of the surrogate model and the acquisition function.As an alternative approach, we are considering optimizing the acquisition function using a different method from quantum annealing, such as simulated annealing.Our method in this paper, which encodes combinations of substituents as a binary vector, can be applied even in a more vast chemical space.Future challenges include verifying the performance in such cases and investigating the computational time advantage of using quantum annealing.

DFT (Density Functional Theory) calculations
For the proposed substituents by the D-Wave quantum annealer, the energy value of ground and excited states were calculated by optimizing the geometry based on DFT calculation.DFT calculations were performed using the supercomputer TSUBAME 3.0 with Gaussian16, Revision C.01 software (Frisch et al., 2019), with the functional B3LYP and basis functions 6-31G.19parameters from the DFT calculation were used to reproduce the experimental values.Here, a prediction model was created using random forest regression.

Figure 2 .
Figure 2. Relationship between the number of loops and the coefficient of determination R 2 in the surrogate model each σ 2 and random sampling.

Figure 3 .
Figure 3.The transition of the best target material property values in the existing dataset up to each loop.

Figure 4 .
Figure 4.The histogram of the target material property values for the combinations of substituents in the initial dataset and those added to the dataset by the end of each experiment.The red dotted line shows the cutoff value (0.880), which we defined as a desired target material property value.