Front. Comput. Sci. Frontiers in Computer Science Front. Comput. Sci. 2624-9898 Frontiers Media S.A. 10.3389/fcomp.2023.1286226 Computer Science Original Research Exploration of new chemical materials using black-box optimization with the D-wave quantum annealer Doi Mikiya 1 * Nakao Yoshihiro 2 Tanaka Takuro 2 Sako Masami 2 Ohzeki Masayuki 1 3 4 5 1Graduate School of Information Sciences, Tohoku University, Sendai, Japan 2LG Japan Lab Ltd., Yokohama, Japan 3International Research Frontier Initiative, Tokyo Institute of Technology, Tokyo, Japan 4Department of Physics, Tokyo Institute of Technology, Tokyo, Japan 5Sigma-i Co. Ltd., Tokyo, Japan

Edited by: David Esteban Bernal Neira, Purdue University, United States

Reviewed by: Ma Hongyang, Qingdao University of Technology, China; Farshud Sorourifar, The Ohio State University, United States

 *Correspondence: Mikiya Doi mikiya.doi.t2@dc.tohoku.ac.jp 12 12 2023 2023 5 1286226 31 08 2023 06 11 2023 Copyright © 2023 Doi, Nakao, Tanaka, Sako and Ohzeki. 2023 Doi, Nakao, Tanaka, Sako and Ohzeki

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

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.

 quantum annealing quantum computing black-box optimization combinatorial optimization problem materials informatics section-at-acceptance Theoretical Computer Science 1 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., 2016; 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; Goto and Ohzeki, 2023; Hasegawa et al., 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; 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.

 2 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.

2.1 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.

 2.2 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(μ)}μ=1D 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.

2.2.1 Construction of surrogate model function

We construct the surrogate model fsurrogate(x→) in the QUBO form in this paper.

fsurrogate(x→)=α0+∑iαixi+∑i<jαijxixj

For simplicity, we set the surrogate model parameters {αi,αij}=α→∈ℝ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 X(μ)=(1,x1(μ),...,xN(μ),x1(μ)x2(μ),x1(μ)x3(μ),...,xN-1(μ)xN(μ)). Furthermore, we set the output vector y→∈ℝD and I as the identity matrix. Then, we assume a prior distribution of surrogate model parameters P(α→) with a variance σα2I and a likelihood function over the surrogate model parameters α→ with a variance σy2I. We give the prior distribution and likelihood function as following multivariate Gaussian distributions.

P(α→)=N(0→,σα2I) P(y→|α→,X)=N(Xα→,σy2I)

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.

α→|y→,X~N(μ→,Σ)μ→=(XTX+λI)-1XTy→Σ=σ2(XTX+λI)-1s.t.σ2=σy2,λ=σy2σα2

We sample the surrogate model parameters α→∈ℝ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.

 2.2.2 Construction of acquisition function

The acquisition function facquisition(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.

x→(new)=argminx→{facquisition(x→)}

facquisition(x→) is a function with modified specific parameters from the surrogate model fsurrogate(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 fsurrogate(x→). For example, when six potential substituents can bind at a specific site, representing the 0-indexed substituent numbers in binary requires three bits (x1, x2, x3). In this context, x1x2x3 = (000, 001, 010, 011, 100, 101)2 corresponds to valid substituent numbers from 0 to 5 in decimal. However, each combination x1x2x3 = (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 x1x2 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 facquisition(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.

 3 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.

3.1 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 x1 to x3, similarly, x4 to x8 represent the substituent number at R2, x9 to x14 represent the substituent number at R3, and x15 to x20 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 2,000 μ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.

By sampling the surrogate model parameters α→ from a multivariate Gaussian distribution N(μ→,Σ) described in (4), construct the surrogate model.

Construct the acquisition function by partially correcting the surrogate model parameters as explained in 2.2.2

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.

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.

Append the new samples {x→(new),y(new)} obtained in (iv) to the existing dataset and return to (i).

 3.2 Experimental results 3.2.1 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.

Histogram of substituent numbers for combinations of substituents added to the dataset by the end of the experiment each σ2. Top left is σ2 = 0, top right is σ2 = 4 × 10−3, bottom left is σ2 = 8 × 10−3, bottom right is σ2 = 12 × 10−3.

 3.2.2 Relationship between the number of loops and the <italic>R</italic><sup>2</sup> of the surrogate model

We show the transition of the coefficient of determination R2 in the surrogate model at each loop in Figure 2. The coefficient of determination R2 is calculated from the initial dataset sample points, 992 points, and the sample points appended up to each loop. Note that R2, 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 R2 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, R2 improves by fitting to similar input vectors and outputs. However, to improve R2 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.

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

 3.2.3 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 Figures 3, 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 exceeding the best target material property value in the initial dataset. In Figure 4, we show the histograms of the target material property values for all combinations of substituents in the initial dataset and those added to the dataset by the end of each experiment.

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

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.

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 Tables 1–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.

Desired target material property value in σ2 = 0.

R1 R2 R3 R4 Target material property value
2 5 0 15 0.881
4 7 8 15 0.889
4 5 8 15 0.889
4 5 0 14 0.902
0 5 0 14 0.897
4 6 8 15 0.894
0 5 8 15 0.887
2 5 8 15 0.902
1 5 0 14 0.905
0 7 8 14 0.882
0 6 0 14 0.898
0 7 0 14 0.905
0 5 14 14 0.911
2 5 0 14 0.88
4 6 0 14 0.915
2 7 0 14 0.888
2 7 8 15 0.892
1 6 0 14 0.91
2 6 0 14 0.89
3 6 0 14 0.884
3 5 8 15 0.894
1 7 8 14 0.887
1 6 8 14 0.891
3 7 8 15 0.891
1 7 0 14 0.914

Desired target material property value in σ2 = 4 × 10−3.

R1 R2 R3 R4 Target material property value
0 5 0 15 0.894
3 5 0 14 0.882
3 5 0 15 0.88
2 5 0 15 0.881
4 7 0 0 0.895
3 5 0 0 0.882
1 5 0 14 0.909
0 5 0 14 0.901
3 6 8 15 0.899
3 7 0 0 0.881
2 7 0 0 0.883

Desired target material property value in σ2 = 8 × 10−3.

R1 R2 R3 R4 Target material property value
2 5 0 15 0.881
2 5 0 14 0.88
4 5 0 15 0.88
4 5 8 15 0.886
0 5 0 14 0.897
1 5 0 14 0.906

Desired target material property value in σ2 = 12 × 10−3.

R1 R2 R3 R4 Target material property value
2 7 14 14 0.901
4 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.

 5 Additional requirements 5.1 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.19 parameters from the DFT calculation were used to reproduce the experimental values. Here, a prediction model was created using random forest regression.

 Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

 Author contributions

MD: Data curation, Formal analysis, Investigation, Methodology, Validation, Visualization, Writing—original draft, Writing—review & editing. YN: Data curation, Funding acquisition, Investigation, Project administration, Resources, Software, Visualization, Writing—review & editing, Methodology. TT: Data curation, Funding acquisition, Investigation, Project administration, Resources, Software, Visualization, Writing—review & editing, Methodology. MS: Investigation, Methodology, Resources, Software, Writing—review & editing. MO: Conceptualization, Funding acquisition, Methodology, Project administration, Resources, Supervision, Writing—review & editing.

 Funding

Our study receives financial support from the MEXT-Quantum Leap Flagship Program Grant No. JPMXS0120352009, as well as Public\Private R&D Investment Strategic Expansion PrograM (PRISM) and programs for Bridging the gap between R&D and the IDeal society (society 5.0) and Generating Economic and social value (BRIDGE) from Cabinet Office.

This study used the TSUBAME3.0 supercomputer at Tokyo Institute of Technology. The authors thank the fruitful discussions with Y. Nishikawa.

 Conflict of interest

YN, TT, and MS are employed by LG Japan Lab Ltd. MO is employed by Sigma-i. The remaining author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

 Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

 References Amin M. H. Andriyash E. Rolfe J. Kulchytskyy B. Melko R. (2018). Quantum Boltzmann machine. Phys. Rev. X 8, 021050. 10.1103/PhysRevX.8.021050 Arai S. Ohzeki M. Tanaka K. (2021). Mean field analysis of reverse annealing for code-division multiple-access multiuser detection. Phys. Rev. Res. 3, 033006. 10.1103/PhysRevResearch.3.033006 Baptista R. Poloczek M. (2018). “Bayesian optimization of combinatorial structures,” in Proceedings of the 35th International Conference on Machine Learning (PMLR) (Stockholm: ICML), 462–471. Chapelle O. Li L. (2011). “An empirical evaluation of Thompson sampling,” in Proceedings of the 24th International Conference on Neural Information Processing Systems, eds J. Shawe-Taylor, R. Zemel, P. Bartlett, F. Pereira, and K. Weinberger (Red Hook, NY: Curran Associates, Inc.), 2249–2257. Das A. Chakrabarti B. K. (2008). Colloquium: Quantum annealing and analog quantum computation. Rev. Modern Phys. 80, 1061–1081. 10.1103/RevModPhys.80.1061 Deisenroth M. P. (2011). A survey on policy search for robotics. Found. Trends Robot. 2, 1–142. 10.1561/2300000021 Denchev V. S. Boixo S. Isakov S. V. Ding N. Babbush R. Smelyanskiy V. . (2016). What is the computational value of finite-range tunneling? Phys. Rev. X 6, 031015. 10.1103/PhysRevX.6.031015 Frisch M. J. Trucks G. W. Schlegel H. B. Scuseria G. E. Robb M. A. Cheeseman J. R. . (2019). Gaussian 16 Revision C.01. Wallingford, CT: Gaussian Inc. Goto T. Ohzeki M. (2023). Online calibration scheme for training restricted Boltzmann machines with quantum annealing. arXiv preprint arXiv:2307.09785. 10.48550/ARXIV.2307.09785 Haba R. Ohzeki M. Tanaka K. (2022). Travel time optimization on multi-AGV routing by reverse annealing. Sci. Rep. 12, 17753. 10.1038/s41598-022-22704-036273242 Hasegawa Y. Oshiyama H. Ohzeki M. (2023). Kernel Learning by quantum annealer. arXiv preprint arXiv:2304.10144. 10.48550/ARXIV.2304.10144 Hatakeyama-Sato K. Kashikawa T. Kimura K. Oyaizu K. (2021). Tackling the challenge of a huge materials science search space with quantum-inspired annealing. Adv. Intell. Syst. 3, 2000209. 10.1002/aisy.202000209 Ide N. Asayama T. Ueno H. Ohzeki M. (2020). “Maximum likelihood channel decoding with quantum annealing machine,” in 2020 International Symposium on Information Theory and Its Applications (ISITA) (IEEE). Inoue D. Okada A. Matsumori T. Aihara K. Yoshida H. (2021). Traffic signal optimization on a square lattice with quantum annealing. Sci. Rep. 11, 1–12. 10.1038/s41598-021-82740-033568714 Jones D. R. Schonlau M. Welch W. J. (1998). Efficient global optimization of expensive black-box functions. J. Glob. Opt. 13, 455–492. 10.1023/A:1008306431147 Kadowaki T. Nishimori H. (1998). Quantum annealing in the transverse Ising model. Phys. Rev E 58, 5355–5363. Kirkpatrick S. Gelatt C. D. Vecchi M. P. (1983). Optimization by simulated annealing. Science 220, 671–680. Kitai K. Guo J. Ju S. Tanaka S. Tsuda K. Shiomi J. . (2020). Designing metamaterials with quantum annealing and factorization machines. Phys. Rev. Res. 2, 013319. 10.1103/PhysRevResearch.2.013319 Koshikawa A. S. Ohzeki M. Kadowaki T. Tanaka K. (2021). Benchmark test of black-box optimization using D-wave quantum annealer. J. Phys. Soc. Japan 90, 064001. 10.7566/JPSJ.90.064001 Lucas A. (2014). Ising formulations of many NP problems. Front. Phys. 2, 5. 10.3389/fphy.2014.00005 Maruyama N. Ohzeki M. Tanaka K. (2021). Graph minor embedding of degenerate systems in quantum annealing. arXiv preprint arXiv:2110.10930. Morita S. Nishimori H. (2008). Mathematical foundation of quantum annealing. J. Math. Phys. 49, 125210. 10.1063/1.2995837 Neukart F. Compostella G. Seidel C. Von Dollen D. Yarkoni S. Parney B. (2017). Traffic flow optimization using a quantum annealer. Front. ICT 4, 29. 10.3389/fict.2017.00029 Ohzeki M. Miki A. Miyama M. J. Terabe M. (2019). Control of automated guided vehicles without collision by quantum annealer and digital devices. Front. Comput. Sci. 1, 9. 10.3389/fcomp.2019.00009 Ohzeki M. Nishimori H. (2011). Quantum annealing: an introduction and new developments. J. Comput. Theoret. Nanosci. 8, 963–971. 10.1166/jctn.2011.1776963 O'Malley D. Vesselinov V. V. Alexandrov B. S. Alexandrov L. B. (2018). Nonnegative/binary matrix factorization with a D-wave quantum annealer. PLoS ONE 13, e0206653. 10.1371/journal.pone.020665330532243 Rosenberg G. Haghnegahdar P. Goddard P. Carr P. Wu K. de Prado M. L. (2016). Solving the optimal trading trajectory problem using a quantum annealer. IEEE J. Select. Top. Signal Processs. 10, 1053–1060. 10.1109/JSTSP.2016.2574703 Shikanai R. Ohzeki M. Tanaka K. (2023). Traffic signal optimization using quantum annealing on real map. arXiv preprint arXiv:2308.14462. 10.48550/ARXIV.2308.14462 Snoek J. Larochelle H. Adams R. P. (2012). “Practical Bayesian optimization of machine learning algorithms,” in Proceedings of the 25th International Conference on Neural Information Processing Systems - Volume 2 (Red Hook, NY: Curran Associates, Inc.), 2951–2959. Suzuki S. Okada M. (2005). Residual energies after slow quantum annealing. J. Phys. Soc. Japan 74, 1649–1652. 10.1143/JPSJ.74.1649 Tanaka T. Sako M. Chiba M. Lee C. Cha H. Ohzeki M. (2023). Virtual screening of chemical space based on quantum annealing. J. Phys. Soc. Japan 92, 023001. 10.7566/JPSJ.92.023001 Thompson W. R. (1933). On the likelihood that one unknown probability exceeds another in view of the evidence of two samples. Biometrika 25, 285–294. Urushibata M. Ohzeki M. Tanaka K. (2022). Comparing the effects of Boltzmann machines as associative memory in generative adversarial networks between classical and quantum samplings. J. Phys. Soc. Japan 91, 074008. 10.7566/JPSJ.91.074008 Vanchinathan H. P. Nikolic I. Bona F. D. Krause A. (2014). “Explore-exploit in top-n recommender systems via Gaussian processes,” in Proceedings of the 8th ACM Conference on Recommender Systems (New York, NY: ACM). Venturelli D. Kondratyev A. (2019). Reverse quantum annealing approach to portfolio optimization problems. Quant. Mach. Intell. 1, 17–30. 10.1007/s42484-019-00001-w Yamamoto M. Ohzeki M. Tanaka K. (2020). Fair sampling by simulated annealing on quantum annealer. J. Phys. Soc. Japan 89, 025002. 10.7566/JPSJ.89.025002 Yonaga K. Miyama M. Ohzeki M. Hirano K. Kobayashi H. Kurokawa T. (2022). Quantum optimization with Lagrangian decomposition for multiple-process scheduling in steel manufacturing. ISIJ Int. 62, 1874–1880. 10.2355/isijinternational.ISIJINT-2022-019