Black-box optimization using factorization and Ising machines

Black-box optimization using factorization and Ising machines

Ryo Tamura,1, a) Yuya Seki,2 Yuki Minamoto,3 Koki Kitai,4 Yoshiki Matsuda,3, b) Shu Tanaka,2, c) and Koji
Tsuda5, d)
1) Center for Basic Research on Materials, National Institute for Materials Science, 1-1 Namiki, Tsukuba, Ibaraki 305-0044,
Japan
2) Graduate School of Science and Technology, Keio University, 3-14-1 Hiyoshi, Kouhoku-ku, Yokohama, Kanagawa 223-8522,
Japan
3) Fixstars Amplify Corporation, 3-1-1 Shibaura, Minato-ku, Tokyo 108-0023, Japan
4) Unprecedented-scale Data Analytics Center, Tohoku University, 6-3 Aoba, Aramakiaza, Aoba-ku, Sendai, Miyagi 980-8578,
Japan
5) Graduate School of Frontier Sciences, The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa 277-8561,
arXiv:2507.18003v1 [cond-mat.stat-mech] 24 Jul 2025

Japan
(*Electronic mail: tsuda@k.u-tokyo.ac.jp)
(*Electronic mail: shu.tanaka@keio.jp)
(*Electronic mail: tamura.ryo@nims.go.jp)
(Dated: 25 July 2025)
Black-box optimization (BBO) is used in materials design, drug discovery, and hyperparameter tuning in machine
learning. The world is experiencing several of these problems. In this review, a factorization machine with quantum
annealing or with quadratic-optimization annealing (FMQA) algorithm to realize fast computations of BBO using
Ising machines (IMs) is discussed. The FMQA algorithm uses a factorization machine (FM) as a surrogate model
for BBO. The FM model can be directly transformed into a quadratic unconstrained binary optimization model that
can be solved using IMs. This makes it possible to optimize the acquisition function in BBO, which is a difficult
task using conventional methods without IMs. Consequently, it has the advantage of handling large BBO problems.
To be able to perform BBO with the FMQA algorithm immediately, we introduce the FMQA algorithm along with
Python packages to run it. In addition, we review examples of applications of the FMQA algorithm in various fields,
including physics, chemistry, materials science, and social sciences. These successful examples include binary and
integer optimization problems, as well as more general optimization problems involving graphs, networks, and strings,
using a binary variational autoencoder. We believe that BBO using the FMQA algorithm will become a key technology
in IMs including quantum annealers.

CONTENTS 1. Encoding integer-variables 5

2. Binary variational autoencoder 6
I. Introduction 2 E. Ising machine 6
1. D-Wave quantum annealer 6
II. FMQA Algorithm 3 2. Fixstars Amplify Annealing Engine 6
A. Black-box optimization 3 3. Digital annealer 7
B. Factorization machine (FM) 4
C. FMQA algorithm for binary optimization III. Implementation and Tools 7
problem 5 1. Original FMQA package 7
D. Encoding method to binary variables 5 2. FMQA implementation using the Amplify
SDK and PyTorch 7
3. Python library: Amplify-BBOpt 8
a) Also at Graduate School of Frontier Sciences, The University of Tokyo, IV. Application Examples 8
5-1-5 Kashiwanoha, Kashiwa 2778561, Japan; Also at RIKEN Center for
Advanced Intelligence Project, 1-4-1 Nihonbashi, Chuo-ku, Tokyo 103-0027, A. Binary optimization problems without constraint 8
Japan 1. Metamaterial design 8
b) Also at Fixstars Corporation, 3-1-1 Shibaura, Minato-ku, Tokyo 108-0023,
2. Layered material design 9
Japan 3. Polymer design 9
c) Also at Department of Applied Physics and Physico-Informatics, Keio Uni-

versity, 3-14-1 Hiyoshi, Kohoku-ku, Yokohama-shi, Kanagawa 223-8522, 4. Analog circuit design 10
Japan; Also at Keio University Sustainable Quantum Artificial Intelligence 5. Feature selection 11
Center (KSQAIC), Keio University, Tokyo 108-8345, Japan; Also at Human B. Binary optimization problems with constraint 11
Biology-Microbiome-Quantum Research Center (WPI-Bio2Q), Keio Univer- 1. Magnetic tunnel junction structure design 11
sity, Tokyo 108-8345, Japan
d) Also at Center for Basic Research on Materials, National Institute for Mate- 2. Atomistic configuration optimization in
rials Science, 1-1 Namiki, Tsukuba, Ibaraki 305-0044, Japan; Also at RIKEN molecules 13
Center for Advanced Intelligence Project, 1-4-1 Nihonbashi, Chuo-ku, Tokyo 3. Structural design for resonance avoidance 13
103-0027, Japan C. Integer variable problems 13
Black-box optimization using factorization and Ising machines 2

1. Structural optimization of photonic crystal obtained, when the number of candidate inputs is small. How-
lasers 13 ever, combinatorial optimization problems are NP-hard23 , and
2. Wing shape optimization 15 as the number of candidate inputs increases, the optimization
3. Signal control optimization 15 of a surrogate model using an exhaustive search is impossi-
4. Design optimization of vehicle body structures ble. Consequently, effectively searching for better surrogate
for multiple vehicles 16 model solutions is difficult when huge number of candidate
D. Others 17 inputs are present.
1. Thermal emitter design 17 From another perspective, the world is full of combinato-
2. Molecular structure design 17 rial optimization problems such as scheduling, matching, and
3. Peptide design 18 route search problems24 . To address these problems, Ising
machines (IMs)25 have been developed that can effectively
V. Closing Statements and Outlook 19 explore the ground state of the Ising models26–28 . The Ising
model was introduced to address critical phenomena in statis-
Acknowledgments 19 tical physics29,30 . Nearly all combinatorial optimization prob-
lems can be expressed using the Ising model. The ground
AUTHOR DECLARATIONS 20 state of the Ising model corresponds to the optimal solution of
Conflict of Interest 20 the combinatorial optimization problem. The Hamiltonian of
Author Contributions 20 Ising model is defined as follows:
N
Data Availability Statement 20 HIsing = ∑ hi si + Ji j si s j , (1)
∑
i=1 1≤i< j≤N
REFERENCES 20
where si = ±1 denotes the classical spin variables and N de-
notes the number of variables. Parameters hi and Ji j can have
I. INTRODUCTION real values known as the magnetic field and magnetic inter-
action, respectively. Because the number of states is 2N , de-
Black-box optimization (BBO) methods1–4 can provide termining the optimal state with the minimum energy of the
better solutions to black-box functions5,6 , where the output Ising model is difficult when N increases. This model can be
value is determined only from the input value, without know- individually translated into
ing the functional form of the function. Examples of such op-
HQUBO = ∑ Qi j xi x j , (2)
timization problems include the optimization of experimental 1≤i≤ j≤N
conditions, simulation parameters, and hyperparameters in the
machine learning (ML) models, which appear everywhere in where xi is a 0/1 binary variable, and this type of model
the basic and social sciences. As a BBO method, the most fa- is known as quadratic unconstrained binary optimization
mous one is the Bayesian optimization7–10 , which has solved (QUBO)31 . To effectively obtain the minimum energy states
several optimization problems in physics11–14 , chemistry15–18 , in Eqs. (1) and (2), that is, the optimal solution to the
and materials science19–22 . combinatorial optimization problem, IMs implement clas-
In BBO methods, ML surrogate models are prepared from sical annealing32–35 , quantum annealing36,37 , or dynamic
training data, where both the input and output are already evolution38,39 as hardware. In addition, an algorithm for solv-
known by black-box functions, such that the output value can ing QUBO by gated quantum computing, known as quantum
be correctly predicted from any inputs using the surrogate approximate optimization algorithm (QAOA), has also been
models. By optimizing the surrogate model, the best solu- developed.40
tion to the surrogate model is obtained. If the surrogate model The use of IMs can solve the difficulties of BBO using
has a better prediction accuracy, then the best solution to the discrete variables41 . In other words, the IM was used to
surrogate model will also be a better solution to the black-box solve the combinatorial optimization problem based on sur-
function. This is the idea behind BBO. rogate model optimization with discrete variables, and better
Let us consider the optimization of the surrogate models. inputs for black-box functions can be identified at high speeds.
When the input is a continuous variable, continuous optimiza- The attempt to use IMs for discrete BBO is roughly catego-
tion methods, such as gradient descent and Newton’s method, rized according to the surrogate model used: factorization ma-
are used to optimize the surrogate models. However, obtain- chine (FM)42 , Bayesian optimization of combinatorial struc-
ing globally optimal solutions is difficult, because these algo- tures (BOCS)43–46 , and the Ising model Hamiltonian itself47 .
rithms tend to be trapped in local minima. In addition, when Because these surrogate models can be directly converted to
the input variables are discrete, the optimization of surrogate QUBO defined by Eq. (2), better solutions to the surrogate
models is a combinatorial optimization problem, and deter- models were instantly obtained by using IMs. This review fo-
mining the best solution is significantly more difficult. Of cuses on the case using an FM that has been used in several
course, because the computation time required to obtain the applications. For the first time, the use of an FM for BBO was
output value from the input using the surrogate models is quite proposed by Kitai et al.41 , and a D-Wave quantum annealer
short, exhaustive search is useful and an optimal solution is was used to solve the FM model. The proposed algorithm
Black-box optimization using factorization and Ising machines 3

IMs. The continuous input can be treated as an integer

optimization problem, if the variables are discretized.

Others: For this category, the input for the black-box

function is not binary and integer variables such as
graphs, networks, and strings. The binary variational
autoencoder (bVAE)56 is useful to convert these inputs
into binary variables. Using the converted binary vari-
ables, FM is trained, and the optimal solution to FM
is selected by IMs. The solution with binary variables
is converted back to inputs for the black-box function
using the decoder.

Section V discusses the future prospects of BBO methods us-

ing IMs.

II. FMQA ALGORITHM

FIG. 1. Flow of the FMQA algorithm that performs BBO using FM

and IMs. In this section, the strategy for combining FM and IMs to
effectively solve BBO problems is presented. First, the gen-
eral problem of BBO is defined. Thereafter, the FM, which
is known as FMQA (factorization machine with quantum an- is used as a surrogate model, is introduced as an elemen-
nealing). In this review, the BBO algorithm using an FM is tary technique, and the FMQA algorithm is described. Sub-
collectively referred to as FMQA (factorization machine with sequently, encoding methods from integer variables and oth-
quadratic-optimization annealing), if the IMs do not use quan- ers to binary variables are introduced to solve the general
tum technologies. This process is summarized in Fig. 1. optimization problem using the FMQA algorithm. In addi-
In Section II, the details of the FMQA algorithm are de- tion, some IMs that are solvers of the surrogate model in the
scribed, along with the elemental techniques used, that is, FM FMQA algorithm are introduced. Finally, the implementation
and IM. In Section III, we present several open-source imple- of FMQA is explained.
mentations of the FMQA algorithm In Section IV, applica-
tion studies of FMQA for the basic and social sciences and
engineering are introduced. For instance, materials designs A. Black-box optimization
using computational materials science methods have been ac-
tively performed41,48,49 . Furthermore, the crystal structure BBO is a method for optimizing black-box functions. The
prediction was considered using FMQA50–52 . In addition, input to the black-box function can be any vector of integer
peptide design in biotechnology53 , optimization of resonance and real values, graphs, networks, or strings, and the input is
avoidance in engineering54 , and traffic signal control in social expressed by x. A black-box function returns only the func-
science55 have also been reported. These application exam- tion value y(x) when x is inputted. In other words, it is a func-
ples can be classified into three types depending on the type tion for which no information useful for optimization com-
of input to the black-box function, as follows (see Fig. 2). putations, such as function forms and gradients, can be ob-
Binary optimization problem: For this category, the tained. Simulations and experimental results in basic science
input is given by binary variables. The FM can be used are generally considered to be black-box functions. The pre-
as a surrogate model for the black-box function, and the diction accuracy of neural networks is also a black-box func-
obtained solution can be used as an input to the black- tion, where the input is hyperparameters. Owing to the lack
box function as it is. of useful information for optimization, BBO uses a surrogate
model that approximates the black-box function. As the sur-
Integer optimization problem: For this category, the rogate model only needs to predict the output from the inputs,
input is given by integer variables. By converting a sin- various ML models can be used. The BBO procedure is sum-
gle integer variable into multiple binary variables using marized below.
an encoder, the FM can be trained using binary vari-
ables as inputs. The optimal solution to an FM using Step 1 For some inputs, the outputs of the black-box functions
binary variables should be decoded into integer vari- are prepared as initial training data. If the number of
ables that are input to the black-box function. To en- initial data is M, the training data are expressed as D =
sure that the solution with binary variables of an FM {(xm , ym )}m=1,...,M .
can be correctly decoded into integer variables for the
black-box function, it is necessary to search for the bet- Step 2 Using D, the surrogate model is trained to correctly pre-
ter solutions to the FM under appropriate constraints by dict the output value from the inputs.
 Black-box optimization using factorization and Ising machines 4

tion, which uses Gaussian process regression as a surrogate

model. Because Gaussian process regression is a stochastic
model, the predictions and their variances, which correspond
to the uncertainty of the predictions, can be evaluated. By
selecting candidates with better predictions, better candidates
for the black-box function are obtained; this is known as “ex-
ploitation.” In contrast, by selecting candidates with a large
variance, candidates far from the training data are selected.
This is known as “exploration.” In Bayesian optimization, an
acquisition function is introduced to simultaneously perform
exploitation and exploration with a single function using both
the predictions and variances obtained from the Gaussian pro-
cess. Rather than optimizing the surrogate model in Step 3,
Bayesian optimization selects the input that maximizes the ac-
quisition function as the next input for the black-box function.
Python packages can be used to perform Bayesian optimiza-
tion, such as Optuna57 , GpyOpt58 , Summit59 , BoTorch57,60 ,
and PHYSBO61,62 . However, as explained in the introduc-
tion, optimizing the acquisition function is a time-consuming
task.

B. Factorization machine (FM)

In the FMQA algorithm, the FM is used as a surrogate

model. The ML prediction model up to the quadratic term
is the FM proposed by Rendle63 . The FM is defined as

N N N
ȳ(x) = w0 + ∑ wi xi + ∑ ∑ ⟨vi , v j ⟩xi x j := HFM , (3)
i=1 i=1 j=i+1

where ⟨vi , v j ⟩ is defined as

K
⟨vi , v j ⟩ = ∑ vik v jk . (4)
k=1

FIG. 2. Three types of optimization problems depending on the type

of input to the black-box function. When the input is binary vari- Here, w0 , {wi }i=1,...,N and {vik }i=1,...,N,k=1,...,K are the model
ables, preparing the encoder and decoder is not necessary. For inte- parameters, and are determined such that ȳ(x) can correctly
ger optimization problems, integer variables should be converted into predict y by x. In this model, only a single hyperparameter
binary variables using an encoder. To perform optimization problems K determines the dimensionality of the model parameters. In
with variables other than binary or integer values, a bVAE is useful other words, the number of model parameters is NK + N + 1.
for the encoder and decoder. An important aspect of the FM is that the computation time
of the third term in Eq. (3) exhibits linear complexity. In
addition, the training of the model parameters can be easily
Step 3 The optimal input of the surrogate model is explored. performed using gradient descent methods with linear com-
In this step, the acquisition function defined by the sur- plexity. For instance, the root mean square loss is considered,
rogate model is often used instead of surrogate model. which is defined as follows:
Step 4 For the selected input, the output value is obtained via a
black-box function. Then, the number of training data L= ∑ (ȳ(xm ) − ym )2 , (5)
is increased as D = {(xm , ym )}m=1,...,M+1 . m=1,...,M

Step 5 Repeats from Steps 2 to 4.

when the training data are D = {(xi , yi )}i=1,...,M . To minimize
The most well-known BBO method is Bayesian optimiza- L, the gradients of y(xi ), which depend on the model parame-
 Black-box optimization using factorization and Ising machines 5

ters, are necessary, and can be easily evaluated as follows: D. Encoding method to binary variables

∂L
= 1, (6) The FMQA algorithm is a specialized optimization method
∂ w0
for binary variables. However, optimization problems ex-
∂L pressed in terms of integers and other variables can also be
= xi , (7)
∂ wi treated by converting them into binary variables using an ap-
∂L N propriate encoding method. This section describes encoding
= xi ∑ v jk − vik xi2 , (8) methods for integer variables and the use of a bVAE.
∂ vik j=1

The computational order for these gradient calculations is

1. Encoding integer-variables
O(1). Thus, all the parameter updates can be performed using
O(NK), and fast training of the FM is guaranteed. The com-
putation time increases linearly with the amount of training Here, we introduce three types of representative encoding
data M. Some packages are available for training FMs, such methods, from integers to binary variables, that is, binary,
as libFM63 , fastFM64 , PyTorch65 , and FMQA66 . one-hot, and domain-wall. In Ref.67 , the optimization per-
Considering the binary variables for xi = ±1, the model formances using three encoding methods were compared for
in Eq. (3) can be easily converted into QUBO. Because the benchmark problems. In addition, other encoding methods
constant term has no meaning for QUBO, except for w0 , the such as the Gray code68 are also used. Note that if continuous
following relationships are established between the FM and values are discretized, they can be treated as integer variables.
QUBO: Thus, the encoding methods introduced here can also handle
( continuous values.
wi (i = j) Binary encoding: A positive integer variable z is expressed
Qi j = , (9) by the binary variables {xi }i=1,...,N where xi = 0 or 1, and N
⟨vi , v j ⟩ (i ̸= j)
denotes the length of the binary variables, as follows
where xi2 = xi . Subsequently, the FM can be solved effectively
using IMs. z = 1 + x1 + 2x2 + ... + 2(N−1) xN . (10)

The integer variable z ranges from 1 to 2N . Depending on

C. FMQA algorithm for binary optimization problem the range of the integer, the dimension of N should be given.
If L integer variables are required to express the optimization
problem, L × N-dimensional binary variables are used as in-
Using the FM as a surrogate model for BBO, an IM can be puts for BBF. Therefore, the FM can be expressed as
used to quickly search for the optimal input for the surrogate
model. This is known as the FMQA algorithm, and the proce- LN LN LN
dure when the input is binary is as follows (see also Fig. 1). HFM = w0 + ∑ wi xi + ∑ ∑ ⟨vi , v j ⟩xi x j . (11)
i=1 i=1 j=i+1
Step 1 For some inputs that are expressed by binary variables
as x = (x1 , ..., xN ), the outputs of black-box functions Thus, by solving HFM using IMs as it is, the better solution to
are prepared as initial training data. When the number the integer optimization problem can be obtained.
of initial data is M, the training data are expressed as One-hot encoding: In addition to the integer encoding, a
D = {(xm , ym )}m=1,...,M . positive integer variable z is expressed by the binary variables
{xi }i=1,...,N as follows:
Step 2 Using D, the FM is trained to correctly predict the out-
put value from the inputs. N
z = ∑ ixi . (12)
Step 3 The optimal input of FM defined by Eq. (3) is explored i=1
by the IM.
under the following constraints:
Step 4 For the selected input, the output value is obtained via a
N
black-box function. Then, the number of training data
is increased as D = {(xi , yi )}i=1,...,M+1 . ∑ xi = 1. (13)
i=1
Step 5 Repeats from Steps 2 to 4. The integer variable z ranges from 1 to N. For instance, when
In Step 3, the input already included in D may become op- N = 4, {1, 0, 0, 0} = 1, {0, 1, 0, 0} = 2, {0, 0, 1, 0} = 3, and
timal. In this case, it is necessary to select a solution that {0, 0, 0, 1} = 4 (see also Fig. 2). If L integer variables are con-
is not included in D by randomly generating, slightly chang- sidered, then L×N-dimensional binary variables are prepared.
ing the optimal input, or adjusting the parameters of the IM. Thus, the FM can be prepared as follows:
In Ref.41 , which proposes the FMQA algorithm, the case of LN LN LN
random generation is suggested. This will work roughly as HFM = w0 + ∑ wi xi + ∑ ∑ ⟨vi , v j ⟩xi x j . (14)
“exploitation” in Bayesian optimization. i=1 i=1 j=i+1
Black-box optimization using factorization and Ising machines 6

In addition, we should consider the constraints given by Using the IMs, the lower energy state of HFM is obtained. Us-
Eq. (13). To prevent this constraint, the QUBO model solved ing the decoder of the trained bVAE, the input for BBF can
by IMs is given by be generated from the binary variables obtained by IMs in the
!2 latent space. Thus, by applying the FMQA algorithm in the
L N
latent space, the BBO can be performed.
HQUBO = HFM + α ∑ ∑ xi+N(l−1) − 1 , (15)
l=1 i=1

where α denotes a positive hyperparameter. The second term E. Ising machine

becomes zero when all constraints are satisfied. The first term
decreases when the surrogate model defined by the FM has a
smaller objective function value. Thus, by solving HQUBO us- The IMs connected to the FMQA algorithm include the
ing IMs, a better solution to the integer optimization problem D-Wave quantum annealer, the Fixstars Amplify Annealing
is obtained. For the one-hot encoding, Endo et al.69 reported Engine (Amplify AE), and the Fujitsu digital annealer (DA).
that when continuous values are considered, introducing func- The following sections describe the characteristics of the ma-
tion smoothing regularization during the training of the FM chines.
model leads to improved optimization performance.
Domain-wall encoding: Domain-wall encoding70,71 resem-
bles one-hot encoding, and the integer value of z can be ex- 1. D-Wave quantum annealer
pressed as follows:
N−1 The D-Wave One was discovered in 2011 as the first com-
z= ∑ i(xi − 2xi xi+1 + xi+1 ) + NxN + 1. (16) mercial quantum computer36 . D-Wave One has 128 qubits,
i=1 which can directly realize quantum annealing72,73 computa-
The integer variable z ranges from 0 to N. For instance, tions in the hardware and solve QUBO. A superconducting
when N = 3, {0, 0, 0} = 1, {1, 0, 0} = 2, {1, 1, 0} = 3, and circuit was used in the D-Wave system in which the right and
{1, 1, 1} = 4 (see also Fig. 2). The integer variable can be rep- left rotations of the superconducting current acted as binary
resented with one bit less than in the case of binary encoding. variables, representing 0 or 1, respectively. The superconduct-
If L integer variables are considered, then L × N-dimensional ing current is indeterminate until it is observed owing to quan-
binary variables are prepared. In addition, we should consider tum superposition, and the quantum annealing computation is
the constraint, and the QUBO model is defined as performed by gradually reducing the quantum effects. The di-
lution refrigerator produces a temperature of 12 mK, which is
HQUBO = HFM close to absolute 0 K, and the shielding layer blocks external
signals such as external magnetic fields and vibrations. The
!
L N N−1
+α ∑ ∑ xi+N(l−1) − ∑ xi+N(l−1) xi+1+N(l−1) D-Wave system operates at less than 25 kW of power, mostly
l=1 i=2 i=1 for the cooling system and the front-end server.
(17) The D-Wave 2000Q, a machine with 2048 qubits arranged
on a coarsely connected graph known as a chimera graph,
where α denotes a positive hyperparameter. The second term was launched in 2017. Although this is a coarsely connected
becomes zero when the all constraints are satisfied. graph, a fully connected graph can be created by making sev-
eral bits act as a single bit in a process known as embed-
2. Binary variational autoencoder
ding. The latest machine is the D-Wave Advantage74 , which
has more than 5,000 qubits and can solve large problems.
In addition, the advantage system and hybrid solver service
A bVAE with binary variables as the latent space is useful can also be combined to run up to 1,000,000 binary variable
as an encoding method for variables other than binary and in- problems for coarsely connected graphs, and up to 20,000
teger variables, such as graphs, networks, and strings to binary binary variable problems for fully connected graphs. In re-
variables. The structure of bVAE is shown in Fig. 2, and was cent years, state-of-the-art techniques such as fast annealing
constructed using two neural networks: an encoder and a de- computations75,76 have been implemented.
coder. The bVAE learns to make the decoded variables equal
to the input using numerous inputs without objective func-
tions. The key point is that the values of the objective func-
2. Fixstars Amplify Annealing Engine
tions are not required to train the bVAE. Using the encoder
part, transforming inputs into binary variables is possible. If
the dimension of the latent space is N, the input variable can The Amplify AE is a GPU-based IM developed and oper-
be expressed as {xi }i=1,...,N where xi = 0 or 1. Using a dataset ated by Fixstars Amplify Corporation in Japan77 . Amplify
in which the objective functions are already known, the FM is AE have been available to the public since 2021, and as of to-
trained as day (October 2024), they have processed more than 45 million
N N N user requests. The IM can solve QUBO using up to 131,072
HFM = w0 + ∑ wi xi + ∑ ∑ ⟨vi , v j ⟩xi x j . (18) (fully connected graphs) and 262,144 (sparse graphs) binary
i=1 i=1 j=i+1 variables. Thus, several practical combinatorial optimization
Black-box optimization using factorization and Ising machines 7

problems, including FMQA, are solvable in terms of the deci- III. IMPLEMENTATION AND TOOLS
sion variable dimensions.
The power of the Amplify AE can be exploited most ef- Several implementations of the FMQA algorithm are avail-
fectively when the users formulate their combinatorial opti- able as open source under an MIT license. These points are
mization problems using the Amplify SDK development en- discussed in this section.
vironment. The Amplify SDK realizes fast and easy formula-
tion of optimization problems. It supports more than ten dif-
ferent solvers, including gate-based quantum computers with 1. Original FMQA package
QAOA, quantum annealing machines, and IMs of other types,
including the Amplify AE. The SDK also supports software-
based optimization solvers, such as Gurobi. The SDK bridges The original implementation of the FMQA algorithm is
the user’s formulation and the selected optimization solver, available at https://github.com/tsudalab/fmqa. To in-
regardless of the decision variables, polynomial order, con- stall the package, obtain fmqa-master.zip from the GitHub
straints, and hardware specifications. For this, the Amplify page and execute the following commands:
SDK performs conversions such as variable conversion, or-
$ python setup.py install
der reduction, penalization of constraints, and graph embed-
ding. Thus, FMQA is no longer limited to binary-variable To train the FM, xs and ys are prepared as the training data.
BBO. It can also solve problems with integer and real vari- The former is a matrix containing bit sequences of the training
ables and explicit constraints without additional implemen- data. The latter is a list of the values of the objective function.
tation of variable encoding or constraint penalization in the The training was performed using the following commands:
resulting QUBO equation.
1 model = fmqa.FMBQM.from_data(xs, ys)

The training was performed using Adam85 . The next candi-

3. Digital annealer date was obtained by solving this model using an IM. For in-
stance, when solving by SA, the dimod package provided by
The DA developed by Fujitsu Limited, is a machine that D-Wave can be used, and the following is an example to ob-
searches for the optimal solution to the Ising model through tain the next candidate.
simulation in a digital circuit. DA is implemented using com-
1 import dimod
plementary metal-oxide-semiconductor (CMOS) technology, 2 sa_sampler = dimod.samplers.
and several hundred thousand (∼ 100, 000) bits have already SimulatedAnnealingSampler()
been realized with all couplings78–80 . Specifically, the hard- 3 res = sa_sampler.sample(model)
ware minimizes the evaluation function while updating the bit
states using a stochastic search based on the Markov chain res contains the bit sequence as the next candidate.
Monte Carlo (MCMC) method. The basic operation is based The coefficients of the trained FM are obtained using the
on the simulated annealing (SA) method81,82 . An acceleration following commands:
technique using digital circuit technology is applied. Typical
1 model.linear[k] for k in model.linear
innovations for quickly obtaining the optimal solution in DA 2 model.quadratic[(k, l)] for (k, l) in model.
are as follows: quadratic
Parallel bit manipulation: A Monte Carlo simulation using
the Metropolis method is effective for determining the opti- The QUBO model was constructed using these coefficients.
mal solution to the Ising model. With this method, the change The QUBO model can be solved using various IMs, and the
in the total energy of the system is calculated when a bit is FMQA algorithm can be applied.
flipped, and the algorithm determines whether the bit should
be flipped. DA attempts to flip all bits in parallel, which im-
proves the speed of the optimization search. 2. FMQA implementation using the Amplify SDK and
Rejection-free search: All the spin states are generated PyTorch
from the transition probabilities calculated from the energy
change of the flip operation for each bit. This improves the Using the aforementioned Amplify SDK, one can im-
search speed because the state is always updated in a single plement FMQA using Python in a straightforward manner.
operation without rejection. Several such examples are available at https://amplify.
Replica exchange method: A replica exchange method is fixstars.com/en/demo#blackbox for a range of BBO
used to reduce the risk of becoming trapped in a local mini- problems such as optimization of material design, chemical
mum and not reaching a global minimum83,84 . This technique plant operation conditions, airfoil design and city traffic sig-
overcomes trapping at the local minimum by simultaneously nal control. These example programs utilize PyTorch86 to
running Monte Carlo calculations for different temperatures train the FM models on par with the original FMQA pack-
(replicas) and implementing an exchange of replicas between age. The programs convert the trained FM model using a Py-
each temperature. Torch class to QUBO formulation, which is then optimized
 Black-box optimization using factorization and Ising machines 8

using an IM such as Amplify AE. Additionally, this imple- Before setting up the optimizer, the solver client must be
mentation allows the use of “hard” constraints such as one- specified for use during the annealing process. Similar to Am-
hot and domain-wall, as well as other equality and inequality plify SDK, Amplify-BBOpt supports various commercially
constraints. Therefore, non-binary variables are implemented available solvers; howeber, here, we use Amplify AE.
in a straightforward manner, and such examples are included
1 from datetime import timedelta
in the above example programs. In addition, the user can con- 2 from amplify import FixstarsClient
sider additional constraints of variables in FMQA using this 3
functionality of the Amplify SDK, regardless of the target op- 4 client = FixstarsClient()
timization solvers. 5 client.parameters.timeout = timedelta(
milliseconds=2000) # annealing timeout
for 2 seconds
6 # client.token = "xxxxxxxxxxx" # Enter your
3. Python library: Amplify-BBOpt Amplify AE API token.

Fixstars Amplify has been developing a Python library Finally, we instantiate an optimizer known as
known as Amplify-BBOpt87 , which facilitates QA-based FMQAOptimizer and start the optimization for five iter-
BBO, such as FMQA. Amplify-BBOpt leverages Fixstars ations.
Amplify’s features explained in the previous section “Fixstars 1 from amplify_bbopt import FMQAOptimizer
Amplify Annealing Engine”. In addition, the above afore- 2
mentioned encoding from non-binary to binary variables is 3 optimizer = FMQAOptimizer(data=data, client=
included. Thus, the user can exploit the fast optimization of client, objective=my_blackbox_func)
various IMs while considering non-binary decision variables 4 optimizer.optimize(num_cycles=5)
without additional encoder/decoder implementation. The user
can install the library as follows: As QA-based BBO methods, such FMQA are relatively
new and are still under maturity. Therefore, Amplify-BBOpt
$ pip install amplify_bbopt plans to implement features along with the research and de-
velopment of relevant studies.
A typical usage of the library for a black-box function com-
prising of “fluid dynamics simulation” of a flow around an
airfoil is shown below. The optimization results and the prob- IV. APPLICATION EXAMPLES
lem detail are introduced in the Sec. 3 as “Wing shape opti-
mization”. The simulator, wing_simulator, takes “width,” A. Binary optimization problems without constraint
“height,” and “attack angle” of an airfoil and computes and
returns the drag and lift forces. The return value from the
If all the inputs, that is, the binary variables x, to the FM
black-box function is the corresponding negative lift-to-drag
are the same as the inputs to the BBF, the FM is minimized
ratio.
by IMs, and the optimal solution to the FM can be used as the
1 from amplify_bbopt import RealVariable, input for the BBF. In this simple case, there is no need to add
blackbox constraints to solve the FM in FMQA algorithm. This type of
2
3 @blackbox problem has been applied to the design of metamaterials41,88 ,
4 def my_blackbox_func( layered materials48,89–91 , and nanoparticles92 . In addition, it
5 wing_width: float = RealVariable(bounds can be applied to the lossy compression of matrices54 and fea-
=(1, 20), nbins=100), ture selection to obtain an ML model with better prediction
6 wing_height: float = RealVariable(bounds accuracy93 . In the following sections, we present some appli-
=(1, 5), nbins=20), cations of metamaterial design, which is the first application
7 wing_angle: float = RealVariable(bounds
=(0, 45), nbins=20), of FMQA, layered material design, and feature selection, as
8 ) -> float: representative examples.
9 lift, drag = wing_simulator(wing_width,
wing_height, wing_angle)
10 return -lift / drag # value to minimize 1. Metamaterial design
The initial training data samples are generated based on the
randomly sampled inputs. In the BBO context, the training Metamaterials exhibit characteristic properties when artifi-
dataset does not need to be large. Here, we start with only ten cially and intricately combined with multiple materials. Be-
samples. cause the material properties are signifiantly dependent on
the elements, compounds, geometry, and their arrangements,
1 from amplify_bbopt import DatasetGenerator BBO is effective when the structure of the metamaterial is the
2
3 num_init_data = 10 input and the properties obtained from simulations or exper-
4 data = DatasetGenerator(objective= iments are the output. However, if the types of elements and
my_blackbox_func).generate(num_samples= compounds increase, and a complex arrangement is consid-
num_init_data) ered, a combinatorial explosion will occur in the candidate
Black-box optimization using factorization and Ising machines 9

materials. Thus, FMQA algorithm is useful for designing resented by 2Nl bits. Optimization calculations were per-
metamaterials. formed up to the case Nl = 24; therefore, a 48-bit problem
As an application of FMQA, a metamaterial design for ra- was treated at most. When the materials to be used as layered
diative cooling was reported by Kitai et al41 . The arrange- materials and their arrangements are known, the wavelength-
ments of three types of rod-like compounds (SiO2 , SiC, and dependent optical properties are calculated using the transfer
PMMA) were optimized to obtain better radiative cooling matrix method (TMM). The thickness of each layer is fixed
properties, as shown in Fig. 3 (a). In this problem setting, at 50 nm. The difference between the calculated and desired
three different materials are present, and it is impossible to optical properties was defined as FOM. This FOM was used
represent them as binary variables. However, by considering as the output of BBF. Using FMQA, the structures of lay-
the constraint that each layer can contain only SiO2 or SiC, it ered materials with a small FOM were explored. The results
is possible to represent the arrangement of the rods as binary of the optimization processes are shown in Fig. 3 (e), when
variables. In other words, as shown in Fig. 3 (a), a binary vari- Nl = 6, 8, 10, 12, 16, 20, 24. For Nl = 16, 20, 24, a total of 5000
able representing Si compounds or PMMA for each position iterations of FMQA optimization were performed. D-Wave
and another binary variable representing SiC = 0 and SiO2 = Advantage 4.1 was used as the IM. For the Nl = 24 case, the
1 in each layer are prepared. Let C denote the number of com- best structure with the smallest FOM was determined, and the
pounds in the horizontal direction and L denote the number layered material was fabricated using a thin-film deposition
of compounds in the vertical direction. The arrangement can process. In the 24 layers in the best structure, some of the
then be expressed as (C + 1) × L bits. When learning the FM, neighboring layers are made of the same materials, and a ma-
we flattened these bits into one dimension and used them as terial with ten layers was proposed, that is, TiO2 (100 nm) /
inputs. SiO2 (150 nm) / TiO2 (100 nm) / SiO2 (150 nm) / TiO2 (100
The emissive power was evaluated using rigorous coupled- nm) / SiO2 (50 nm) / Al2 O3 (100 nm) / TiO2 (200 nm) / Al2 O3
wave analysis (RCWA) calculations. The difference between (150 nm) / TiO2 (100 nm). In this material, Si3 N4 was not
the calculated properties and atmospheric window was de- used. This material was fabricated, and its energy-dispersive
fined as the figure of merit (FOM). This FOM was used as X-ray spectroscopy map is shown in Fig. 3 (f), resulting in
the output of BBF, and a structure with a large value of FOM a 10 layered material being correctly obtained. It was con-
is required. To maximize the properties in FMQA, the output firmed using transmitted irradiance measurements that this
of BBF is defined as y = −FOM and the FMQA algorithm material had the lowest FOM among the materials reported
was applied. In Fig. 3 (b), the results of the optimization pro- in the literature48 .
cesses are shown when L = 6 and C = 3, that is, 24 bits prob- The same technique is used to design a wide-angle spectral
lem was considered. As an IM, the D-Wave 2000Q was used, filter for energy saving windows. In this case, the solar inci-
and more than 2000 BBO cycles were performed. For com- dent angle dependence of the optical properties was evaluated
parison, random search results are also shown. If the FMQA using the TMM, and the averaged FOM for the angle depen-
algorithm is used, structures with large FOM values are found dence was used as the output of BBF. The same four materials
even if the number of calculated structures is smaller than that were used in addition to the high-performance transparent ra-
obtained by random sampling (RS). In this study, optimiza- diative cooling design, and the structural design using FMQA
tions for larger systems with a maximum bit size of 60 have up to Nl = 20 was performed. When Nl = 20, 8000 opti-
also been attempted. The metamaterials suitable for radiative mization cycles were performed, and a world record material
cooling have also been developed. In addition, the computa- among other materials in the literature was also obtained90 .
tion times of BBO with and without D-Wave 2000Q are sum-
marized in Fig. 3 (c). Without using IMs, the computation
time required to select candidates increases rapidly depending 3. Polymer design
on the problem size, whereas the search time remains almost
constant as the problem size is increased when IM is used.
The polymer was designed using the FMQA algorithm94 .
The results indicate that complex metamaterial structures can
A total of 32 polymer motifs were prepared in advance and
be optimized within a short time using the FMQA algorithm.
represented as 5-bit binary variables (215 = 32), that is, the
concept is similar to that in the case where the number of ma-
terials in a layered material is large. Three polymer motifs
2. Layered material design were selected to create triblock polymers (see Fig. 4 (a)) to op-
timize structures with high thermal conductivity (TC). Thus,
Materials exploration has been performed using FMQA, the triblock polymer could be represented in 15 bits. The TC
which focuses on layered materials48 . Layered materials were was obtained from a trained deep neural network (DNN) us-
designed by layering multiple types of films (see Fig. 3 (d)). ing Morgan fingerprints95 with frequency as a feature. A total
Layered materials for high-performance transparent radiative of 1144 polymer datasets were used to train the DNN. In the
cooling were designed by combining four materials: SiO2 , QUBO, there is no need to consider the additional constraints,
Si3 N4 , Al2 O3 , and TiO2 . Each material can be represented and only the FM model is included to predict the TC. A solver
by two bits, that is, the materials are defined as SiO2 = 00, provided by D-Wave, dimod96 , was used to solve QUBO. The
Si3 N4 = 01, Al2 O3 = 10, and TiO2 = 11. If the number iterative dependence of the largest TC is shown in the Fig. 4
of layers is Nl , then the one layered material can be rep- (b), and the TC gradually increases with the number of gen-
Black-box optimization using factorization and Ising machines 10

FIG. 3. (a) Optimization target of metamaterial structure constructed by three types rod-like materials (SiO2 , SiC, and PMMA) to obtain better
radiative cooling property. This structure can be translated to binary variables, and emissive power was calculated by RCWA calculations.
The FOM can evaluate the shape of emissive power, and when all finite emissive power is located within the atmospheric window, FOM is
maximized. (b) FOM results depending on the number of calculated structures. The results by FMQA and RS were compared. (c) Computation
time of BBO with and without using the D-Wave 2000Q depending on the system size (number of encoding bits). (d) Structure of layered
materials expressed by binary variables. (e) FOM results depending on the number of calculated structures depending on the number of layers
Nl in the layered materials optimization. (f) Energy-dispersive X-ray spectroscopy map for the fabricated layered materials. (For (a), (b), and
(c) reprinted from ref.41 under a Creative Commons Attribution 4.0 International license. For (d), (e), and (f), adapted with permission from
ref.48 . Copyright 2022 American Chemical Society.)

erations. However, FMQA did not exhibit better optimization 4. Analog circuit design
performance than Bayesian optimization or the genetic algo-
rithm (GA). Optimization performance can be improved by
An analog circuit was designed using FMQA algorithm97 .
using IMs and improved encoding methods.
The objective was to optimize the conductive structure of the
filter used in the power amplifier, as shown in Fig. 4 (c). In a
conductive structure with input and output ports, optimization
is performed by targeting a specific frequency response. Two-
Black-box optimization using factorization and Ising machines 11

dimensional (2D) square blocks are used to define the con- measure prediction accuracy. Fig. 4 (f) shows the cycle de-
ductive structure. For each block, binary variables are defined pendence of the mean squared error by FMQA algorithm and
such that if the conductor is present, it is 1; if the conductor RS when the tensile modulus is a materials property. Even if
is absent, it is 0. A 22 × 22 arrangement of blocks is con- the number of cycles is small, the FMQA algorithm can iden-
sidered with input and output ports on both sides. Given this tify a better combination of features that exhibits a small mean
structure, the finite element method (FEM) was used to sim- squared error. The comparison plots between the real tensile
ulate the frequency response of the conductive structure. The modulus and the prediction results obtained using the optimal
FOM is defined as the value that decreases when the desired model are summarized in Fig. 4 (g). FMQA can successfully
frequency is reached, and optimization is performed to reduce select features with a higher prediction accuracy than LASSO.
it. In the QUBO for FMQA, there is no need to consider ad- Therefore, this study shows that the FMQA algorithm can be
ditional constraints, and only the FM model is considered. In applied as a feature selection method.
this study, a DA was used as an Ising model solver. In addi-
tion, an optimization method that alternates between FMQA
and GA is proposed to sample various structures. The itera-
tive dependence of the FOM is illustrated in Fig. 4 (d). The B. Binary optimization problems with constraint
results are compared with those of particle swarm optimiza-
tion (PSO). The method combining FMQA and GA obtained Here, we introduce some applications of optimization prob-
the desired structure with fewer computations than PSO. lems that require constraints in addition to the FM model. This
is the case when not all the inputs to the FM are valid, and
the solution must be selected with an appropriate constraint.
5. Feature selection By minimizing the QUBO, which adds constraint terms to the
FM, only valid solutions satisfying the constraint can be se-
Feature selection is essential for improving the prediction lected. In the following section, we demonstrate the assign-
accuracy of ML. It has been reported that the FMQA algo- ment of constraint terms to specific problem examples.
rithm can also be applied to feature selection93 . The least ab-
solute shrinkage and selection operator (LASSO) method is a
well-known feature selection method used to avoid overfitting
in linear regression98 . For a more precise implementation, the 1. Magnetic tunnel junction structure design
exhaustive search method99,100 is useful for determining bet-
ter combinations of features with the highest prediction ac- Magnetic tunnel junctions (MTJs) have garnered consider-
curacy. However, the total number of combinations becomes able attention in the field of spintronics, because they are ex-
2N − 1 when the number of features is N. Thus, as the num- pected to be applicable in non-volatile memory devices. The
ber of features is increased, a combinatorial explosion occurs. MTJ structure comprises two layers of magnetic material with
This feature selection can be considered as a BBO with binary a very thin insulating layer between them. An optimization
variables. In other words, the input of BBF is a combination study combining FMQA and first-principles calculations was
of features, and the prediction accuracy of the ML models conducted to understand the suitable structure of the insulat-
is used as the output. To solve this problem using BBO, N- ing layer49 . The target material was MgGa2 O4 103 , which has
dimensional binary variables {x1 , ..., xN } with xi = 0 or 1 are an inverse spinel structure as shown in Fig. 5 (a). The op-
prepared. If xi = 0, the ith feature is not used, and if xi = 1, timization target is the arrangement of atoms to improve the
the ith feature is used. The feature selection was performed physical properties when Mg or Ge ions are placed at the dis-
using by the FMQA algorithm. ordering sites. A structure with 10 disordered sites was con-
In Ref.93 , feature selection using the FMQA algorithm was sidered: Mg ions were placed in five of them, and Ge ions
performed to prepare the ML model with a high prediction were placed in the remaining five. In other words, there are
accuracy for the mechanical properties of polymer materials. 10 C5 = 252 possible configurations. By providing a binary
A total of 75 samples of homo polypropylene materials were variable xi = 0 or 1 for each site, xi = 0 when Mg ions are
collected. The target features were the parameters extracted placed and xi = 1 when Ge ions are placed. Considering the
from X-ray diffraction (XRD) patterns using Bayesian spec- case where there are N disordered sites in total, the arrange-
tral deconvolution101,102 . Examples of peak deconvolution for ment of the ions can be expressed by the binary variables:
XRD measurements in the machine direction (MD; flow di- x = {x1 , x2 , ..., xN }. By training the FM defined in Eq. (3) with
rection of the resin) and transverse direction (TD; width di- known y(x) as the target physical property, a surrogate model
rection of the resin) are shown in Fig. 4 (e). Five feature for predicting this property was trained. However, there is
types were generated for each peak, with 24 peaks for each no guarantee that the solution minimizing the FM will have
sample. Thus, the number of dimensions of the prepared fea- five ions each. In other words, when the FM is solved by
tures was 120, and it was difficult to exhaustively search for the IM, the solution is selected from the 2N possible struc-
2120 − 1 candidates. In this study, the mean squared error be- tures expressed by x instead of from the 10 C5 configurations.
tween the predicted and real mechanical properties was used Therefore, to ensure that the obtained solution is correct for
as the prediction accuracy, which is the output of the BBF. the number of ions, that is, the numbers of xi = 1 and xi = 0
A linear regression model and cross-validation were used to are the same, a QUBO with the following constraint was pre-
Black-box optimization using factorization and Ising machines 12

FIG. 4. (a) Triblock polymers defined by 15 bits. (b) Cycle dependence of the best TC for 20 independent runs. (c) Example of conductive
structures where 1 is input and 2 is output ports, respectively. (d) Cycle dependence of the objective function defined by frequency response.
Present results were obtained by the combination method of FMQA and GA, and benchmark results were by PSO. (e) Examples of peak
deconvolution results for XRD measurements in the machine direction (MD; flow direction of the resin) and transverse direction (TD; width
direction of the resin) of the homo polypropylene materials. (f) Cycle dependence of the mean squared error between the predicted tensile
modulus and the real one by FMQA algorithm and RS for the homo-polypropylene materials. (g) Prediction results of the optimal model for
the tensile modulus. For comparison, the results by LASSO are also plotted. (For (a) and (b), reprinted from ref.94 under a Creative Commons
Attribution (CC BY) license. For (c) and (d), reprinted from ref.97 under a Creative Commons Attribution-NonCommercial-NoDerivs license.
For (e), (f), and (g), reprinted from ref.93 under a Creative Commons Attribution license.)

pared: must be set appropriately, and used to optimize the MTJ struc-
!2 ture, where α = 1 is assumed.
N
N
HQUBO = HFM + α ∑ xi − 2 , (19) In this study, three types of MTJ properties were con-
i=1
sidered. The target properties were the total energy differ-
where α denotes a positive hyperparameter. The second term ence (∆Etotal ), tunneling magnetoresistance (TMR), and the
is a constraint that is minimized when ∑Ni=1 xi = N/2. For resistance area product, which were calculated using first-
this Hamiltonian, the optimal solution is one in which the out- principles calculations. Small values for these parameters are
put value of the FM is sufficiently small and ∑Ni=1 xi = N/2. required for spintronics devices. Therefore, to train the FM,
Parameter α determines the constraint strength. If it is ex- the properties multiplied by a negative sign were used as ob-
cessively small, a solution that does not satisfy the constraint jective functions. The optimization results obtained by FMQA
is obtained, whereas if it is excessively large, a solution with are shown in Fig. 5 (b), which shows the number of cycles re-
better properties is not obtained. Therefore, these parameters quired to determine the best MTJ structure. The results were
Black-box optimization using factorization and Ising machines 13

compared with those obtained using the D-Wave quantum an- 3. Structural design for resonance avoidance
nealer, SA, and random sampling. In addition, the results
were obtained by the strategy in which the minimization of Printed circuit boards are used in the power control units of
the FM was performed by an exhaustive search (ES) without hybrid vehicles and require a proper structural design to avoid
using IMs and Bayesian optimization, because the size of the the resonance. For instance, if the location of the mounting
problem was small. The FMQA algorithm is more efficient holes can be optimized to increase the natural frequency, low-
than random search, and the optimization of ∆Etotal and TMR frequency resonance can be avoided. To solve the optimiza-
requires fewer computations than the Bayesian optimization tion problem of increasing the natural frequency, the FMQA
strategy. This research shows that the FMQA algorithm is algorithm was adopted105 . In this study, the candidate posi-
also a powerful tool for optimizing the arrangement of ions tions of mounting holes were determined on a 2D substrate. If
with the desired physical properties, if a constraint is required. N candidate positions exist, then N-dimensional binary vari-
ables x = {x1 , x2 , ..., xN } were used to determine the arrange-
ment of the mounting holes. If xi = 1, a mounting hole is
2. Atomistic configuration optimization in molecules formed at that position; when xi = 0, no mounting hole is
formed at that position. After the number of mounting holes
The FMQA algorithm was used to optimize the atomistic is determined, the natural frequency is calculated using the
configuration of carbon materials104 . Carbon materials are FEM. To obtain a higher natural frequency, the FM is trained
used such as catalysts in fuel cells and electrode materials in using the natural frequency multiplied by the negative sign as
Li ion batteries, and their industrial applications are highly the objective function. When the number of mounting holes is
active. This study focused on perylene (C20 H12 ), and sub- fixed, additional constraints are required. For instance, when
stitution of carbon atoms was considered. Substitution was the number of mounting holes to be used is Nc , the QUBO
performed by replacing N or B atoms with C atoms. For each model constructed by the FM and constraint is as follows:
carbon position, two binary variables are used to determine
!2
the atomistic configuration. In other words, if the first bit is N
0, no substitution is performed, and if it is 1, substitution is HQUBO = HFM + α ∑ xi − Nc , (21)
performed for the N or B atoms. If the second bit is 0, the i=1
C atom is substituted with a B atom; if the second bit is 1,
the C atom is substituted with a N atom. Thus, for each C where α denotes a positive hyperparameter.
atom position, ‘00’ and ‘01’ mean C atom, ‘10’ is B atom, Fig. 5 (e) shows the results of the optimization for prob-
and ‘11’ is C atom. In perylene, N = 20 carbons are present; lems with N = 20 and Nc = 6. Here, the results of FMQA,
thus, the substituted structure is represented by 40 binary vari- BOCS, and RS are compared, and the D-Wave Hybrid Solver
ables ({xi }i=1,....,2N ). The polarizability calculated using the was used as the IM. The FMQA algorithm can determine a
molecular orbital method was used as the objective function, better mounting hole location with a higher natural frequency
and maximization was desirable. Therefore, when training the and fewer calculations than RS. The results of the frequency
FM, polarizability multiplied by a negative sign was used as analysis obtained at the optimal mounting hole position are
the objective function. In addition to the MTJ structure de- shown in Fig. 5 (f).
sign, better structures were optimized for polarizability when
the number of substituted atoms was fixed. To ensure that
the number of substituted atoms was Nc , the following QUBO C. Integer variable problems
model was developed:
!2 If the integer variable problems are solved using the FMQA
N algorithm, then encoding from integers to binary variables
HQUBO = HFM + α ∑ x2i−1 − Nc , (20) must be used. Some applications were reported55,67,106,107 ,
i=1
and the following representative examples, the structural
where α denotes a positive hyperparameter. This is because optimization of photonic crystal lasers108 , wing shape
the first bit of each carbon position expresses substitution. optimization106 , signal control optimization55 , vehicle body
The results for the case of 6-site substitution (Nc = 6) are structures optimization109 are introduced in this review.
shown in Fig. 5 (c). The number of cycles was 96, and the op-
timized structure and its polarizability values were compared
using the FMQA algorithm and RS. A DA was used as an IM. 1. Structural optimization of photonic crystal lasers
The structure obtained by FMQA had a larger polarizability
value than that obtained by RS. In addition, Fig. 5 (d) shows The demand for semiconductor lasers110 with high output
the obtained objective function values against the computa- power and beam quality is rapidly increasing for various ap-
tion time using the FMQA and GA. The results demonstrate plications, including optical detection in smart mobility and
that the FMQA algorithm can find a better molecular struc- high-precision laser processing in smart manufacturing. The
ture with large polarizability within a short computation time. FMQA algorithm is used to design the next generation of
It has been reported that FMQA exhibits high performance in semiconductor lasers, known as photonic crystal lasers, which
solving molecular replacement problems. are characterized by their ability to operate at high power
Black-box optimization using factorization and Ising machines 14

FIG. 5. (a) Structure of MgGa2 O4 . The disordered site is the orange site. (b) Optimization results by FMQA in MTJ structure design depending
on the algorithms for each physical property. FM + QA, FM + SA, and FM + ES are BBOs using FM with D-Wave quantum annealer, SA,
and ES, respectively. BO and RS are the Bayesian optimization and random sampling. (c) Optimization results of atomistic configuration in
perylene by RS and FMQA when 6-site substitution. The polarizability values are also shown. The green, gray, and brown atoms are B, N, and
C atoms, respectively. (d) Objective function values against the computation time by FMQA and GA for the molecular structure optimization.
(e) Optimization results of the structural design for resonance avoidance by FMQA, BOCS, and RS. (f) Result of frequency analysis obtained
at the optimal mounting hole position. (For (a) and (b), reprinted from ref.49 under a Creative Commons Attribution 4.0 International license.
For (c) and (d), reprinted from ref.104 with permission. For (e) and (f), reprinted from ref.105 under a Creative Commons Attribution 4.0
International license.)

and high beam quality108 . In this study, 11 continuous vari- defined below is used as the objective function:
ables were considered to determine the structure of photonic
P×η
lasers. Continuous variables cannot be handled by FMQA as Q= . (22)
is. Thus, the upper and lower bounds of each variable are θx θy
determined, and the discrete variables are defined by divid- The value of Q increases as P and η increase, and θx and θy
ing them into 16 equal distances between the upper and lower decrease. Because photonic crystal lasers with large Q val-
bounds. Using binary encoding, the 16 values are represented ues are desirable, the properties multiplied by a negative sign
by four binary variables, and 11 discrete parameters can be were used as objective functions to train the FM. In this case,
represented by 44 binary variables. To simultaneously con- all inputs to the FM are valid, that is, all states defined by
sider the three properties of output power (P), polarization the 44 binary variables can be converted into 11-dimensional
ratio (η), and divergence angle (θx and θy ), the quantity Q continuous variables for photonic crystal lasers.
Fig. 6 (a) shows the cycle dependence of Q. The IM used
was a D-Wave Advantage quantum annealer. It is shown that
the FMQA algorithm can be used to find parameters with
Black-box optimization using factorization and Ising machines 15

larger Q values and fewer computations than other optimiza- 3. Signal control optimization
tion methods. Other methods include PSO111,112 and GAs113 .
The band edge frequency of the structure determined using Signal control using the FMQA algorithm was considered
the FMQA algorithm is shown in Fig. 6 (b). For compar- as an example of an application to a social science problem55 .
ison, the results for the initial structure are also presented. In this application, the control of traffic signals was optimized
The band edge frequency in the structure determined by the to maximize the average speed of all vehicles traveling in a
FMQA algorithm is larger than that of the initial structure. city. The objective function was the average vehicle speed,
This study demonstrated the effectiveness of FMQA for opti- which was evaluated using a multi-agent simulation (MAS).
mization problems in which continuous values are discretized Based on a simulation model known as the optimal speed
using binary encoding. model115,116 , the driver of each vehicle accelerated and de-
celerated according to the optimal speed determined from the
distance between the vehicles. The city, as an optimization
target, is shown in Fig. 6 (e). The city has 4 × 4 sections di-
2. Wing shape optimization vided by roads and two shopping malls. There are 16 intersec-
tions in the city. Traffic signals were present at each intersec-
The wing shape of an aircraft is optimized to maximize the tion. The traffic signal can be controlled by three parameters:
lift and minimize the drag. These wing forces can be evalu- the length of the red light Lr , length of the green light Lg , and
ated using actual experiments and computational fluid dynam- the start time of the red light Ls , in seconds. For instance, if
ics simulations. For this optimization problem using the fluid the control parameters for a traffic signal at an intersection are
dynamics simulations, the FMQA algorithm was applied106 . Lr = 15.0, Lg = 12.0, and Ls = 5.0, the north–south signal at
The Joukowski transform was used to determine the shape the intersection will be red 5.0 s after the simulation start time,
of the wing, which is known as Joukowski airfoil. Three con- red for 15.0 s, and green for 12.0 s. When using the FMQA
tinuous parameters (ξ , η, α) are used to determine shape. algorithm, continuous variables must be expressed as binary
Thus, the task of the optimization problem is to find a better variables. One-hot encoding was used. The lower and upper
set of continuous parameters that yields the maximum lift and limits of the three parameters were set to 1 s and 20 s, respec-
minimum drag forces. To evaluate these forces, a 2D fluid tively, and the discrete values were divided into 20 equal parts,
analysis method based on the lattice Boltzmann method114 that is, every second. Thus, when the one-hot encoding was
was used. Here, FL and FD denote the lift and drag forces, re- considered, 20 bits were prepared as binary variables for each
spectively. Using these forces, the goal is to maximize the lift– continuous parameter. The control parameters for each inter-
drag ratio defined by rLD = FL /FD . When training the FM, section can be expressed using 60 binary variables, and 960
−rLD was used as the objective function. When the FMQA al- binary variables ({xi }i=1,...,960 ) are required for all the control
gorithm is applied, continuous variables must be discretized, signals in the city.
and ξ is divided into 20 equal parts between 1.0 and 9.55, η In MAS, a problem setting in which the maximum num-
is divided into 40 equal parts between 0 and 9.75, and α is ber of vehicles is 20 and 50% of these vehicles have des-
divided into 40 equal parts between 0 and 39. In this study, tinations in one of the two shopping malls was considered.
one-hot encoding was used, and a binary variable was pre- The destinations of the remaining 50% of vehicles were ran-
pared for each part. Thus, 20 bits were used for ξ , whereas 40 domly selected from the city. The homes for each vehicle
bits were used for η and α. In other words, 100 bits were used were randomly located in the city, and a simple simulation
to determine the shape of the wing as {xi }i=1,...,100 . To imple- was performed in which the vehicles continued to travel back
ment one-hot encoding, a QUBO model with the following and forth between the shopping mall and their homes or be-
constraint was prepared: tween their destinations and their homes. When the control
parameters of the signals are determined, the average speed
 !2 of vehicles, v̄, was evaluated using MAS. Thus, the control
20
parameters that maximize v̄ are searched for using the FMQA
HQUBO = HFM + α  ∑ xi − 1 algorithm. When training the FM, −v̄ is used as the objec-
i=1
tive function. To implement one-hot encoding, we prepared a
!2 !2  QUBO model with the following constraints:
60 100
+ ∑ xi − 1 + ∑ xi − 1  , (23)
16 3 20+20( j−1)+60(l−1)
!2
i=21 i=41
HQUBO = HFM + α ∑ ∑ ∑ xi − 1 (24)
l=1 j=1 i=1+20( j−1)+60(l−1)
where α denotes positive hyperparameter. Fixstars Amplify
AE is used to solve this QUBO model. The values of −rLD In this study, Fixstars Amplify AE was used as the IM. The
in the optimization process are shown in Fig. 6 (c). In this simulation was performed using 10 random parameters as the
process, the first 20 cycles were performed randomly and the initial data. BBO with the FMQA algorithm was executed
remaining 20 cycles were performed using the FMQA algo- for 90 cycles. The −v̄ values for the optimization cycles are
rithm. Fig. 6 (d) shows the optimized wing shape after 40 shown in Fig. 6 (f). This indicates that the use of FMQA im-
cycles, which improved the lift–drag ratio by 17% compared proved the search for better control parameters for traffic sig-
with that obtained after 20 cycles by RS. nals. The MAS simulation results for the optimized parame-
Black-box optimization using factorization and Ising machines 16

FIG. 6. (a) Cycle dependence of the Q value in the structural optimization of photonic crystal lasers. The results by FMQA, PSO, and GAs are
compared. (b) Band edge frequency difference between initial structure and the optimized structure obtained by FMQA. (c) Cycle dependence
of −rLD for the wing shape optimization by FMQA. (d) Optimized wing shape by FMQA. (e) Optimization target city for signal control
optimization problem. Squares represent the vehicles, whereas white circles represent the homes. Red and green circles represent the signals,
and two malls are present. (f) Cycle dependence of the average of vehicle speeds v̄ by FMQA. (For (a) and (b), reprinted from ref.108 under a
Creative Commons license. For (c) and (d), reprinted from ref.106 with permission. For (e) and (f), reprinted from ref.55 with permission.)

ters are available at https://amplify.fixstars.com/en/ most difficult one is the multi-objective optimization of ve-
demo/fmqa_4_traffic. hicle body structures for multiple vehicles, which is based on
real-world problem settings and complexity. The problem in-
volves:
4. Design optimization of vehicle body structures for • A black-box objective function, which is the total
multiple vehicles weight of the three vehicle models (to be minimized),
• An objective function that can be expressed in a
The benchmark problem for the optimization of vehicle quadratic manner, representing the number of common
body structures for multiple vehicles was based on Mazda’s parts across the three models (to be maximized) and
production models proposed in Ref.117 and is available at
https://ladse.eng.isas.jaxa.jp/benchmark/. The • A total of 18 constraints among which there are 14
benchmark comprises three problem classes; however, the black-box constraint functions per vehicle (total of 42
Black-box optimization using factorization and Ising machines 17

black-box constraints). These black-box constraints are optimized for use in a thermophotovoltaics (TPV) engines.
related to crashworthiness, eigenvalues, body rigidity, The TPV engine generates electricity through radiative heat
and manufacturing constraints. The values of these transfer between a heater and an array of photovoltaic cells.
black-box functions are required to be zero or positive The design of an efficient thermal radiator determines the ef-
values such that the constraints are met. ficiency of the TPV engine. Thus, an efficient design is es-
sential. The target topological shape of the thermal radia-
A total of 222 real variables were used to control the ob- tor is 2D, as shown in Fig. 7 (a). The white colored area is
jectives and constraints. Various groups have attempted to TiN, whereas the black area is filled with air. This topological
solve this benchmark problem based on various GAs includ- shape is treated as a 64 × 64 binarized figure, and 5,000 train-
ing elitist non-dominated sorting genetic algorithm II (NSGA- ing data points were available. When the topological shape is
II)118,119 and Bayesian optimization algorithms120 . All these converted into a binary variables, 4,096 bits are needed. To
experiments consider the number of cycles (number of eval- reduce the dimensions of the explanatory variables, the bVAE
uations of black-box functions) of 10,000 or 30,000, which is useful.
is excessively large for the typical real-world design process Using the bVAE, a 500-dimensional binarized latent vari-
standards. able space was trained. Thus, each topological shape could be
Kondo et al.109 have tackled this problem with only less expressed by 500-dimensional binary variables. The objective
than 1,000 evaluations of black-box functions by adding the function in the optimization problem was defined by the effi-
following two improvements to the original FMQA algorithm: ciency of the thermal radiator. This efficiency was determined
as the product of the in-band and out-of-band efficiencies, and
• Weights to the objectives and constraints are determined
higher values are desirable for the TPV engines. When train-
based on the distances between the objective (con-
ing the FM, efficiency multiplied by the negative sign was
straint) values and corresponding target values. The
used as the objective function. In this case, all states defined
target values are zero for the constraints, whereas for
by the 500-dimensional binary variables can be decoded to
the objective functions, the target values are the aver-
the topological shape of the thermal radiator using the trained
age values in the dataset minus their standard deviation.
bVAE. Thus, there is no need to consider constraints to solve
Thus, the weights for the objectives change over cycles.
the QUBO model. SA and a D-Wave Advantage hybrid sam-
• Multiple searches per cycle. In each optimization cycle, pler were used as IMs. Fig. 7 (b) shows the efficiency as a
the following three FMs are constructed. The first FM function of the iterations using a D-Wave Advantage hybrid
is constructed using a relatively large K (see Eq. (3)) to sampler. It was shown that a structure with better efficiency
make the resulting model over-fit. The second FM is could be obtained using the FMQA algorithm. The results of
based on a relatively small K, expected to be an “under- the top 100 thermal emitter designs for both IMs are shown
fitted” model. The third model is constructed by lin- in Fig. 7 (c). The optimal structures are shown in Fig. 7 (a).
early combining the first and the second FM models. This study reported that better topological structures could be
Then, annealing is performed for each of these three effectively identified by combining bVAE and the FMQA al-
FM models to obtain three input vectors. gorithm.

With these additions, Kondo et al. have performed three in-

dependent evaluations of black-box functions per cycle, mul- 2. Molecular structure design
tiplied by the number of cycles, which was 300, resulting
in a total of 900 function evaluations. The Pareto solutions
Molecular structures can be represented by strings, known
obtained exhibited improved solutions by the NSGA-II with
as SMILES (simplified molecular input line entry system)124 .
30,000 function evaluations.
For instance, benzene is represented as C1=CC=CC=C1 and
acetic acid is represented as CC(=O)O. The length of the
SMILES varies depending on the molecule. To convert
D. Others
SMILES to binary variables, junction tree variational autoen-
coder (JTVAE) is useful for treating molecular structures as
Using the bVAE, various optimization problems can be graphs125 . JTVAE was modified to the bJTVAE using Gum-
solved even if they are not binary or integer optimization bel softmax reparameterization126 to obtain the latent space
problems121 . Applications in thermal emitter design122 , with binary variables. A total of 250,000 molecules were used
molecular structure design123 , and peptide design53 are intro- to train the bJTVAE. The dimension of the latent space was
duced. 300123 .
The three types of molecular properties (penalized LogP,
TPSA, GSK3β + JNK3 + QED + SA) were optimized to
1. Thermal emitter design verify the performance of the combined bJTVAE and the
FMQA algorithm123 . These properties were calculated from
The first study that used the bVAE and FMQA algorithms the molecular structures, and the larger values were required.
was the topology structure optimization of a thermal radia- Thus, to train the FM, the negative sign was multiplied by
tor by Wilson et al.122 . In this study, a thermal radiator was these properties. Fig. 7 (d) compares the results of the op-
Black-box optimization using factorization and Ising machines 18

FIG. 7. (a) Target topological shape of the thermal radiator for the thermal emitter design. (b) Cycle dependence on efficiencies by FMQA.
(c) Results of the top 100 thermal emitter designs by D-Wave quantum annealer with hybrid solvers (QCH) and SA. (d) Optimization results
in molecular structure design. Optimization results by combining bVAE and FMQA (bVAE + IM), and randomly generated in the latent space
of bVAE (bVAE-Random) are shown. For comparison, the results using Bayesian optimization (VAE-BO) and RS (VAE-Random) using a
non-binary VAE are shown. The distribution of initial data is shown. (e) Ranking of Pareto solutions for multiobjective optimizations. (f)
Peptide optimization results that were experimentally synthesized and evaluated. Sequences, optical density, and hemolytic activity for four
peptides are shown. (For (a), (b), and (c), reprinted from ref.122 with permission. For (d), reprinted from ref.123 under a Creative Commons
Attribution 3.0 Unported license. For (e) and (f), reprinted from ref.53 under a Creative Commons Attribution-NonCommercial-NoDerivatives
4.0 International license.)

timizations combining bVAE and FMQA and randomly gen- 3. Peptide design
erated in the latent space of bVAE. For comparison, it also
shows the results using a non-binary VAE. Bayesian optimiza-
tion and RS were performed using a non-binary VAE. In all Peptide design was performed using bVAE and the FMQA
cases, the combined search for bJTVAE and FMQA generated algorithm53 . A peptide is a combination of amino acids, and
several molecules with excellent properties. For this verifica- the 20 amino acids can be represented by strings such as ‘A,’
tion, Fixstars Amplify AE was used as the IM. In addition, ‘C,’ and ‘D.’ Thus, the peptide sequence was represented by
this study reported that increasing the number of binary vari- a string of letters. To convert a sequence into binary vari-
ables in the latent space is expected to improve the accuracy ables, bVAE is useful for text modeling127 . The training data
of bVAE, and IM with several binary variables would be ef- for the sequences included 19,530 antimicrobial and 5,583
fective for molecular structure optimization. non-antimicrobial peptides. The latent space has 64 dimen-
sions; thus, each peptide sequence can be represented by 64-
dimensional binary variables.
In this study, peptide sequence optimization was performed
to simultaneously increase antimicrobial and non-hemolytic
Black-box optimization using factorization and Ising machines 19

properties. These are generally known to have trade-offs, and computers increases, large-scale BBO problems can be han-
designing peptides that simultaneously achieve both is an im- dled using FMQA and quantum computers instead of IMs.
portant issue in the development of therapeutic peptides. To Optimization of the experimental results: The applications
perform this optimization, it is necessary to set up a single introduced in this review are all optimization results when
objective function that handles multi-objective optimization. the objective functions are obtained through simulations. Us-
Multi-objective optimization is performed to find the maxi- ing the results obtained from real experiments as the objec-
mum number of Pareto solutions when two or more objective tive function, they can contribute directly to the basic science,
functions with a trade-off are present. The Pareto solution is such as materials development and drug discovery. However,
the optimal solution after the balance of the objective func- although FMQA can handle a large search space, it requires
tion is determined. Depending on this balance, a different several BBF readings, which is time-consuming. Therefore,
solution becomes a Pareto solution. This set of Pareto solu- when real experiments are conducted in optimization cycles,
tions is known as the Pareto front and is a trade-off curve. FMQA algorithm is expected to be more powerful in com-
To treat the multi-objective optimization problem with a sin- bination with automated experimental systems such as self-
gle objective function, a ranking variable rm is introduced into driving laboratories. Recently, self-driving laboratories have
the currently available data, which indicates the distance be- been developed for various fields130–135 . The linkage between
tween each data point and the Pareto front. If the solution is these factors and FMQA may lead to the development of in-
on the Pareto front, rm = 1 is assigned. After removing all the novative materials and drug discovery strategies.
solutions with rm = 1, a new Pareto front is defined in the re- Applications to larger-scale problems: Up to 1,000,000 bi-
maining samples, and rm = 2 is assigned to the sample on the nary variables can be used in the current IMs. In contrast,
new Pareto front. This procedure is repeated until a ranking the application problems discussed in this review can all be
is obtained for all samples (see Fig. 7 (e)). To train the FM, regarded as small scale problems. Therefore, IMs can poten-
ym = −1/rm is used as the objective function; thus, a small tially be applied to larger-scale problems using the FMQA al-
value of ym is desired. In this way, multi-objective optimiza- gorithm. To address the lager-scale problem, room for further
tion can be treated as single-objective optimization and solved innovation exists in the FMQA algorithm. One is the develop-
using the FMQA algorithm. The BBO method using FMQA ment of effective batch processing methods, in which multiple
that can handle this multi-objective optimization problem is suggestions of inputs for BBF are effectively generated. For
known as multi-objective optimization by quantum annealing instance, batch versions of Bayesian optimization have been
(MOQA). MOQA was used to design 200,000 peptides. The widely developed in the field of informatics136–138 . However,
D-Wave Advantage System 4.1 was used as IM. Four peptides no definitive method has been developed yet, and this is still
were synthesized and evaluated (Fig. 7 (f)) from peptides de- the subject of research. Similarly, effective batch techniques
signed using MOQA. Two of these peptides were experimen- for the FMQA algorithm will definitely be needed in the fu-
tally determined to be both antimicrobial and non-hemolytic. ture. In particular, IM-specific batch techniques, such as us-
This study showed that peptides were successfully designed ing uncertainty in QA and sampling methods using IM, may
by multi-objective optimization using the FMQA algorithm. be effective and will be an important technological develop-
ments. The second is the development of a surrogate model
suitable for BBO using IMs, which is a more advanced ver-
V. CLOSING STATEMENTS AND OUTLOOK sion of the FM. Although not directly solvable using the cur-
rent IM, the introduction of higher order binary optimization
In this review, we provide an overview of the BBO algo- (HUBO)139,140 , which has a higher order than QUBO, is use-
rithm using IMs, known as FMQA, and its applications. In ful for ensuring the prediction accuracy. Furthermore, model
this algorithm, the FM, which can be expressed as a QUBO fitting to the graph structure of IMs can be expected to im-
model, is used as a surrogate model. The FM can then be prove the performance of IMs as well as the FMQA.
solved efficiently using the IMs, and a fast BBO is realized.
In particular, using IMs, the FMQA algorithm has been ex-
plored in a large input space. To make the FMQA algorithm ACKNOWLEDGMENTS
easy to run, a Python package will be developed and released;
thus, R & D using FMQA will become increasingly popular The authors thank Taiga Hayashi, Ryo Ogawa, Tetsuro
in both academia and industry. We strongly believe that BBO Abe, Mayumi Nakano, Tokiya Fukuda, Masashi Yamashita,
using FMQA is a key technique for IMs. Finally, outlook of Hyakka Nakada, Shuta Kikuchi, Kotaro Terada, Yosuke
the FMQA algorithm is briefly described. Mukasa, Kotaro Tanahashi, Tadashi Kadowaki, Tengfei Luo,
Using a quantum computer instead of IM: Gao et al. de- and Junichiro Shiomi for the discussions about FMQA algo-
signed a highly efficient organic light-emitting diode emitter rithm.
using the FMQA algorithm128 . This study used the QAOA40 This work was partially supported by Japan Science and
to solve the QUBO. QAOA can be treated using quantum Technology Agency (JST) PRESTO (JPMJPR24T8), JST
computers. Although the small-sized problems with six bits CREST (JPMJCR21O2), the Japan Society for the Promotion
were attempted, the QAOA calculations were performed us- of Science (JSPS) KAKENHI (Grant Numbers JP23H05447,
ing a state vector simulator and the IBM Quantum System 25K07172), the Council for Science, Technology, and In-
One129 . In the future, when the number of qubits in quantum novation (CSTI) through the Cross-ministerial Strategic In-
Black-box optimization using factorization and Ising machines 20

novation Promotion Program (SIP), “Promoting the appli- 6 N. Wiener, Cybernetics or Control and Communication in the Animal and
cation of advanced quantum technology platforms to social the Machine (The MIT Press, 1961).
7 P. I. Frazier, “A tutorial on Bayesian optimization,” (2018),
issues” (Funding agency: QST), JST (Grant Number JP-
arXiv:1807.02811 [stat.ML].
MJPF2221). In addition, this paper is partially based on re- 8 J. Snoek, H. Larochelle, and R. P. Adams, “Practical Bayesian optimiza-
sults obtained from a project, JPNP23003, commissioned by tion of machine learning algorithms,” in Advances in Neural Information
the New Energy and Industrial Technology Development Or- Processing Systems, Vol. 25, edited by F. Pereira, C. Burges, L. Bottou,
ganization (NEDO). S.T. wishes to express their gratitude to and K. Weinberger (Curran Associates, Inc., 2012).
9 B. Shahriari, K. Swersky, Z. Wang, R. P. Adams, and N. de Freitas, “Tak-
the World Premier International Research Center Initiative ing the human out of the loop: A review of Bayesian optimization,” Pro-
(WPI), MEXT, Japan, for their support of the Human Biology- ceedings of the IEEE 104, 148–175 (2016).
Microbiome-Quantum Research Center (Bio2Q). 10 X. Wang, Y. Jin, S. Schmitt, and M. Olhofer, “Recent advances in

Bayesian optimization,” ACM Comput. Surv. 55 (2023).

11 A. Seko, A. Togo, H. Hayashi, K. Tsuda, L. Chaput, and I. Tanaka, “Pre-

diction of low-thermal-conductivity compounds with first-principles an-

AUTHOR DECLARATIONS harmonic lattice-dynamics calculations and Bayesian optimization,” Phys-
ical Review Letters 115, 205901 (2015).
12 S. Ju, T. Shiga, L. Feng, Z. Hou, K. Tsuda, and J. Shiomi, “Designing
Conflict of Interest
nanostructures for phonon transport via Bayesian optimization,” Physical
Review X 7, 021024 (2017).
13 S. Jalas, M. Kirchen, P. Messner, P. Winkler, L. Hübner, J. Dirkwinkel,
The authors have no conflicts to disclose.
M. Schnepp, R. Lehe, and A. R. Maier, “Bayesian optimization of a laser-
plasma accelerator,” Physical Review Letters 126, 104801 (2021).
14 M. Yu, S. Yang, C. Wu, and N. Marom, “Machine learning the hubbard
Author Contributions U parameter in DFT+U using Bayesian optimization,” npj Computational
Materials 6, 180 (2020).
15 F. Häse, L. M. Roch, C. Kreisbeck, and A. Aspuru-Guzik, “Phoenics:
Ryo Tamura: Conceptualization (lead); Writing – origi- A Bayesian optimizer for chemistry,” ACS Central Science 4, 1134–1145
nal draft (lead); Supervision (equal). Yuya Seki: Writing – (2018).
16 K. Homma, Y. Liu, M. Sumita, R. Tamura, N. Fushimi, J. Iwata,
review & editing (supporting). Yuki Minamoto: Writing –
original draft (supporting); Writing – review & editing (sup- K. Tsuda, and C. Kaneta, “Optimization of a heterogeneous ternary
Li3 PO4 –Li3 BO3 –Li2 SO4 mixture for Li-ion conductivity by machine
porting). Koki Kitai: Writing – review & editing (support-
learning,” The Journal of Physical Chemistry C 124, 12865–12870 (2020).
ing). Yoshiki Matsuda: Supervision (supporting); Writing 17 L. Fang, E. Makkonen, M. Todorović, P. Rinke, and X. Chen, “Effi-
– review & editing (supporting). Shu Tanaka: Conceptual- cient amino acid conformer search with Bayesian optimization,” Journal
ization (supporting); Supervision (equal); Writing – review & of Chemical Theory and Computation 17, 1955–1966 (2021).
18 K. Wang and A. W. Dowling, “Bayesian optimization for chemical prod-
editing (supporting). Koji Tsuda: Conceptualization (sup-
ucts and functional materials,” Current Opinion in Chemical Engineering
porting); Supervision (equal); Writing – review & editing 36, 100728 (2022).
(supporting). 19 Y. K. Wakabayashi, T. Otsuka, Y. Krockenberger, H. Sawada, Y. Taniyasu,

and H. Yamamoto, “Machine-learning-assisted thin-film growth: Bayesian

optimization in molecular beam epitaxy of SrRuO3 thin films,” APL Ma-
terials 7, 101114 (2019).
DATA AVAILABILITY STATEMENT 20 Y. Zhang, D. W. Apley, and W. Chen, “Bayesian optimization for mate-

rials design with mixed quantitative and qualitative variables,” Scientific

Data sharing is not applicable to this article as no new data Reports 10, 4924 (2020).
21 R. Tamura, T. Osada, K. Minagawa, T. Kohata, M. Hirosawa, K. Tsuda,
were created or analyzed in this study.
and K. Kawagishi, “Machine learning-driven optimization in powder man-
ufacturing of Ni-Co based superalloy,” Materials and Design 198, 109290
(2021).
22 Q. Liang, A. E. Gongora, Z. Ren, A. Tiihonen, Z. Liu, S. Sun, J. R.
REFERENCES
Deneault, D. Bash, F. Mekki-Berrada, S. A. Khan, K. Hippalgaonkar,
B. Maruyama, K. A. Brown, J. Fisher III, and T. Buonassisi, “Bench-
1 D.
marking the performance of Bayesian optimization across multiple exper-
Golovin, B. Solnik, S. Moitra, G. Kochanski, J. Karro, and D. Sculley, imental materials science domains,” npj Computational Materials 7, 188
“Google vizier: A service for black-box optimization,” in Proceedings of (2021).
the 23rd ACM SIGKDD International Conference on Knowledge Discov- 23 A. Lucas, “Ising formulations of many NP problems,” Frontiers in Physics
ery and Data Mining, KDD ’17 (Association for Computing Machinery, 2 (2014).
New York, NY, USA, 2017) pp. 1487–1495. 24 W. J. Cook, W. H. Cunningham, W. R. Pulleyblank, and A. Schrijver,
2 S. Alarie, C. Audet, A. E. Gheribi, M. Kokkolaras, and S. Le Digabel,
Combinatorial Optimization (John Wiley & Sons, Inc., 1998).
“Two decades of blackbox optimization applications,” EURO Journal on 25 N. Mohseni, P. L. McMahon, and T. Byrnes, “Ising machines as hardware
Computational Optimization 9, 100011 (2021). solvers of combinatorial optimization problems,” Nature Reviews Physics
3 K. Terayama, M. Sumita, R. Tamura, and K. Tsuda, “Black-box opti-
4, 363–379 (2022).
mization for automated discovery,” Accounts of Chemical Research 54, 26 S. Tanaka, R. Tamura, and B. K. Chakrabarti, Quantum Spin Glasses,
1334–1346 (2021). Annealing and Computation, 1st ed. (Cambridge University Press, USA,
4 W. Kumagai and K. Yasuda, “Black-box optimization and its applications,”
2017).
in Innovative Systems Approach for Facilitating Smarter World, edited by 27 K. Tanahashi, S. Takayanagi, T. Motohashi, and S. Tanaka, “Application
T. Kaihara, H. Kita, S. Takahashi, and M. Funabashi (Springer Nature of Ising machines and a software development for Ising machines,” Journal
Singapore, Singapore, 2023) pp. 81–100. of the Physical Society of Japan 88, 061010 (2019).
5 V. Belevitch, “Summary of the history of circuit theory,” Proceedings of

the IRE 50, 848–855 (1962).

Black-box optimization using factorization and Ising machines 21

28 S. Yarkoni, E. Raponi, T. Bäck, and S. Schmitt, “Quantum annealing for 50 T. Luo, Z. Xu, W. Shang, S. Kim, and E. Lee, “QALO: Quantum
industry applications: Introduction and review,” Reports on Progress in annealing-assisted lattice optimization,” (2024).
Physics 85, 104001 (2022). 51 Y. Couzinié, Y. Seki, Y. Nishiya, H. Nishi, T. Kosugi, S. Tanaka, and Y.-i.
29 K. G. Wilson, “The renormalization group: Critical phenomena and the Matsushita, “Machine learning supported annealing for prediction of grand
kondo problem,” Rev. Mod. Phys. 47, 773–840 (1975). canonical crystal structures,” Journal of the Physical Society of Japan 94,
30 A. Pelissetto and E. Vicari, “Critical phenomena and renormalization- 044802 (2025).
group theory,” Physics Reports 368, 549–727 (2002). 52 J. Lin, T. Tada, A. Koizumi, M. Sumita, K. Tsuda, and R. Tamura, “Deter-
31 G. Kochenberger, J.-K. Hao, F. Glover, M. Lewis, Z. Lü, H. Wang, and mination of stable proton configurations by black-box optimization using
Y. Wang, “The unconstrained binary quadratic programming problem: a an Ising machine,” The Journal of Physical Chemistry C 129, 2332–2340
survey,” Journal of Combinatorial Optimization 28, 58–81 (2014). (2025).
32 M. Yamaoka, C. Yoshimura, M. Hayashi, T. Okuyama, H. Aoki, and 53 A. Tučs, F. Berenger, A. Yumoto, R. Tamura, T. Uzawa, and K. Tsuda,

H. Mizuno, “A 20k-spin Ising chip to solve combinatorial optimization “Quantum annealing designs nonhemolytic antimicrobial peptides in a dis-
problems with CMOS annealing,” IEEE Journal of Solid-State Circuits crete latent space,” ACS Medicinal Chemistry Letters 14, 577–582 (2023).
51, 303–309 (2016). 54 T. Kadowaki and M. Ambai, “Lossy compression of matrices by black box
33 S. Matsubara, H. Tamura, M. Takatsu, D. Yoo, B. Vatankhahghadim, optimisation of mixed integer nonlinear programming,” Scientific Reports
H. Yamasaki, T. Miyazawa, S. Tsukamoto, Y. Watanabe, K. Takemoto, 12, 15482 (2022).
and A. Sheikholeslami, “Ising-model optimizer with parallel-trial bit-sieve 55 “Black-box optimiztion of traffic light control,” (), https://amplify.

engine,” in Complex, Intelligent, and Software Intensive Systems, edited by fixstars.com/en/demo/fmqa_4_traffic.

L. Barolli and O. Terzo (Springer International Publishing, Cham, 2018) 56 F. Mena and R. Ñanculef, “A binary variational autoencoder for hashing,”

pp. 432–438. in Progress in Pattern Recognition, Image Analysis, Computer Vision, and
34 H. Goto, “Bifurcation-based adiabatic quantum computation with a non- Applications, edited by I. Nyström, Y. Hernández Heredia, and V. Mil-
linear oscillator network,” Scientific Reports 6, 21686 (2016). ián Núñez (Springer International Publishing, Cham, 2019) pp. 131–141.
35 H. Goto, K. Tatsumura, and A. R. Dixon, “Combinatorial optimization by 57 T. Akiba, S. Sano, T. Yanase, T. Ohta, and M. Koyama, “Optuna: A

simulating adiabatic bifurcations in nonlinear hamiltonian systems,” Sci- next-generation hyperparameter optimization framework,” in Proceedings
ence Advances 5, eaav2372 (2019). of the 25th ACM SIGKDD International Conference on Knowledge Dis-
36 M. W. Johnson, M. H. S. Amin, S. Gildert, T. Lanting, F. Hamze, N. Dick- covery and Data Mining (2019).
son, R. Harris, A. J. Berkley, J. Johansson, P. Bunyk, E. M. Chapple, C. En- 58 “Gpyopt,” https://sheffieldml.github.io/GPyOpt/.

derud, J. P. Hilton, K. Karimi, E. Ladizinsky, N. Ladizinsky, T. Oh, I. Per- 59 K. C. Felton, J. G. Rittig, and A. A. Lapkin, “Summit: Benchmarking

minov, C. Rich, M. C. Thom, E. Tolkacheva, C. J. S. Truncik, S. Uchaikin, machine learning methods for reaction optimisation,” Chemistry–Methods
J. Wang, B. Wilson, and G. Rose, “Quantum annealing with manufactured 1, 116–122 (2021).
spins,” Nature 473, 194–198 (2011). 60 M. Balandat, B. Karrer, D. R. Jiang, S. Daulton, B. Letham, A. G. Wil-
37 C. Easttom, “Quantum annealing,” in Hardware for Quantum Computing son, and E. Bakshy, “BoTorch: A framework for efficient Monte-Carlo
(Springer Nature Switzerland, Cham, 2024) pp. 101–112. Bayesian optimization,” (2020), arXiv:1910.06403 [cs.LG].
38 T. Inagaki, Y. Haribara, K. Igarashi, T. Sonobe, S. Tamate, T. Honjo, 61 T. Ueno, T. D. Rhone, Z. Hou, T. Mizoguchi, and K. Tsuda, “COMBO:

A. Marandi, P. L. McMahon, T. Umeki, K. Enbutsu, O. Tadanaga, H. Take- An efficient bayesian optimization library for materials science,” Materials
nouchi, K. Aihara, K. Kawarabayashi, K. Inoue, S. Utsunomiya, and Discovery 4, 18–21 (2016).
H. Takesue, “A coherent Ising machine for 2000-node optimization prob- 62 Y. Motoyama, R. Tamura, K. Yoshimi, K. Terayama, T. Ueno, and

lems,” Science 354, 603–606 (2016). K. Tsuda, “Bayesian optimization package: PHYSBO,” Computer Physics
39 T. Honjo, T. Sonobe, K. Inaba, T. Inagaki, T. Ikuta, Y. Yamada, T. Kazama, Communications 278, 108405 (2022).
K. Enbutsu, T. Umeki, R. Kasahara, K. Kawarabayashi, and H. Takesue, 63 S. Rendle, “Factorization machines with libFM,” ACM Trans. Intell. Syst.

“100,000-spin coherent Ising machine,” Science Advances 7, eabh0952 Technol. 3, 57:1–57:22 (2012).
(2021). 64 I. Bayer, “fastfm: A library for factorization machines,” Journal of Ma-
40 E. Farhi, J. Goldstone, and S. Gutmann, “A Quantum Approximate Opti- chine Learning Research 17, 1–5 (2016).
mization Algorithm,” (2017), arXiv:1411.4028 [quant-ph]. 65 A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan,
41 K. Kitai, J. Guo, S. Ju, S. Tanaka, K. Tsuda, J. Shiomi, and R. Tamura, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf,
“Designing metamaterials with quantum annealing and factorization ma- E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner,
chines,” Physical Review Research 2, 013319 (2020). L. Fang, J. Bai, and S. Chintala, “Pytorch: An imperative style, high-
42 S. Rendle, “Factorization machines,” in 2010 IEEE International Confer- performance deep learning library,” in Advances in Neural Information
ence on Data Mining (2010) pp. 995–1000. Processing Systems 32 (Curran Associates, Inc., 2019) pp. 8024–8035.
43 R. Baptista and M. Poloczek, “Bayesian optimization of combinatorial 66 “GitHub, fmqa,” https://github.com/tsudalab/fmqa.

structures,” (2018), arXiv:1806.08838 [stat.ML]. 67 Y. Seki, R. Tamura, and S. Tanaka, “Black-box optimization for integer-
44 A. S. Koshikawa, M. Ohzeki, T. Kadowaki, and K. Tanaka, “Benchmark variable problems using Ising machines and factorization machines,”
test of black-box optimization using D-Wave quantum annealer,” Journal (2022), arXiv:2209.01016 [cs.LG].
of the Physical Society of Japan 90, 064001 (2021). 68 E. Agrell, J. Lassing, E. Strom, and T. Ottosson, “On the optimality of the
45 A. Okada, H. Yoshida, K. Kidono, T. Matsumori, T. Takeno, and T. Kad- binary reflected gray code,” IEEE Transactions on Information Theory 50,
owaki, “Design optimization of noise filter using quantum annealer,” IEEE 3170–3182 (2004).
Access 11, 44343–44349 (2023). 69 K. Endo and K. Z. Takahashi, “Function smoothing regularization for pre-
46 K. Morita, Y. Nishikawa, and M. Ohzeki, “Random postprocessing for cision factorization machine annealing in continuous variable optimization
combinatorial Bayesian optimization,” Journal of the Physical Society of problems,” Physical Review Research 7, 013149 (2025).
Japan 92, 123801 (2023). 70 N. Chancellor, “Domain wall encoding of discrete variables for quan-
47 K. Hatakeyama-Sato, T. Kashikawa, K. Kimura, and K. Oyaizu, “Tack- tum annealing and QAOA,” Quantum Science and Technology 4, 045004
ling the challenge of a huge materials science search space with quantum- (2019).
inspired annealing,” Advanced Intelligent Systems 3, 2000209 (2021). 71 J. Berwald, N. Chancellor, and R. Dridi, “Understanding domain-wall
48 S. Kim, W. Shang, S. Moon, T. Pastega, E. Lee, and T. Luo, “High- encoding theoretically and experimentally,” Philosophical Transactions of
performance transparent radiative cooler designed by quantum comput- the Royal Society A: Mathematical, Physical and Engineering Sciences
ing,” ACS Energy Letters 7, 4134–4141 (2022). 381, 20210410 (2023).
49 K. Nawa, T. Suzuki, K. Masuda, S. Tanaka, and Y. Miura, “Quantum an- 72 T. Kadowaki and H. Nishimori, “Quantum annealing in the transverse

nealing optimization method for the design of barrier materials in magnetic Ising model,” Physical Review E 58, 5355–5363 (1998).
tunnel junctions,” Physical Review Applied 20, 024044 (2023). 73 G. E. Santoro and E. Tosatti, “Optimization using quantum mechanics:

Quantum annealing through adiabatic evolution,” Journal of Physics A:

Black-box optimization using factorization and Ising machines 22

Mathematical and General 39, R393 (2006). 94 X. Huang and S. Ju, “Tutorial: AI-assisted exploration and active design
74 “D-Wave, Advantage2,” https://www.dwavesys.com/ of polymers with high intrinsic thermal conductivity,” Journal of Applied
solutions-and-products/systems/. Physics 135 (2024).
75 A. D. King, S. Suzuki, J. Raymond, A. Zucca, T. Lanting, F. Altomare, 95 H. L. Morgan, “The generation of a unique machine description for chem-

A. J. Berkley, S. Ejtemaee, E. Hoskinson, S. Huang, E. Ladizinsky, A. J. R. ical structures-A technique developed at chemical abstracts service.” Jour-
MacDonald, G. Marsden, T. Oh, G. Poulin-Lamarre, M. Reis, C. Rich, nal of Chemical Documentation 5, 107–113 (1965).
Y. Sato, J. D. Whittaker, J. Yao, R. Harris, D. A. Lidar, H. Nishimori, and 96 “GitHub, dimod,” https://github.com/dwavesystems/dimod.

M. H. Amin, “Coherent quantum annealing in a programmable 2,000 qubit 97 M. Hida, H. Ikeda, A. Maruo, M. Sato, and T. Yamazaki, “Topology op-

Ising chain,” Nature Physics 18, 1324–1328 (2022). timization of analog circuit design via global optimization using factor-
76 A. D. King, J. Raymond, T. Lanting, R. Harris, A. Zucca, F. Altomare, A. J. ization machines with digital annealer,” Journal of Advanced Mechanical
Berkley, K. Boothby, S. Ejtemaee, C. Enderud, E. Hoskinson, S. Huang, Design, Systems, and Manufacturing 18, JAMDSM0076 (2024).
E. Ladizinsky, A. J. R. MacDonald, G. Marsden, R. Molavi, T. Oh, 98 R. Tibshirani, “Regression Shrinkage and Selection via the Lasso,” Journal

G. Poulin-Lamarre, M. Reis, C. Rich, Y. Sato, N. Tsai, M. Volkmann, J. D. of the Royal Statistical Society. Series B (Methodological) 58 (1996).
Whittaker, J. Yao, A. W. Sandvik, and M. H. Amin, “Quantum critical 99 K. Nagata, J. Kitazono, S. Nakajima, S. Eifuku, R. Tamura, and M. Okada,

dynamics in a 5,000-qubit programmable spin glass,” Nature 617, 61–66 “An exhaustive search and stability of sparse estimation for feature selec-
(2023). tion problem,” IPSJ Online Transactions 8, 25–32 (2015).
77 “Fixstars Amplify,” (), https://amplify.fixstars.com/en/. 100 Y. Igarashi, H. Takenaka, Y. Nakanishi-Ohno, M. Uemura, S. Ikeda, and
78 “Fujitsu digital annealer,” https://www.fujitsu.com/global/ M. Okada, “Exhaustive search for sparse variable selection in linear re-
services/business-services/digital-annealer/. gression,” Journal of the Physical Society of Japan 87, 044802 (2018).
79 S. Matsubara, M. Takatsu, T. Miyazawa, T. Shibasaki, Y. Watanabe, 101 K. Nagata, S. Sugita, and M. Okada, “Bayesian spectral deconvolution

K. Takemoto, and H. Tamura, “Digital annealer for high-speed solving with the exchange Monte Carlo method,” Neural Networks 28, 82–89
of combinatorial optimization problems and its applications,” in 2020 25th (2012).
Asia and South Pacific Design Automation Conference (ASP-DAC) (2020) 102 T. Matsumura, N. Nagamura, S. Akaho, K. Nagata, and Y. Ando, “High-

pp. 667–672. throughput XPS spectrum modeling with autonomous background sub-
80 M. Aramon, G. Rosenberg, E. Valiante, T. Miyazawa, H. Tamura, and traction for 3d5/2 peak mapping of SnS,” Science and Technology of Ad-
H. G. Katzgraber, “Physics-inspired optimization for quadratic uncon- vanced Materials: Methods 3, 2159753 (2023).
strained problems using a digital annealer,” Frontiers in Physics 7 (2019). 103 H. Sukegawa, Y. Kato, M. Belmoubarik, P.-H. Cheng, T. Daibou, N. Shi-
81 S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi, “Optimization by simulated momura, Y. Kamiguchi, J. Ito, H. Yoda, T. Ohkubo, S. Mitani, and
annealing,” Science 220, 671–680 (1983). K. Hono, “MgGa2 O4 spinel barrier for magnetic tunnel junctions: Coher-
82 V. Černý, “Thermodynamical approach to the traveling salesman problem: ent tunneling and low barrier height,” Applied Physics Letters 110, 122404
An efficient simulation algorithm,” Journal of Optimization Theory and (2017).
Applications 45, 41–51 (1985). 104 Y. Suga, A. Maruo, and H. Jippo, “A feasibility study for quantum
83 R. H. Swendsen and J.-S. Wang, “Replica Monte Carlo simulation of spin- computing methodologies in automotive advanced material investigation,”
glasses,” Physical Review Letters 57, 2607–2609 (1986). Transactions of Society of Automotive Engineers of Japan 55, 621–627
84 K. Hukushima and K. Nemoto, “Exchange Monte Carlo method and appli- (2024).
cation to spin glass simulations,” Journal of the Physical Society of Japan 105 T. Matsumori, M. Taki, and T. Kadowaki, “Application of qubo solver us-

65, 1604–1608 (1996). ing black-box optimization to structural design for resonance avoidance,”
85 D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” Scientific Reports 12, 12143 (2022).
(2014), arXiv:1412.6980. 106 “Black-box optimization of the airfoil geometry with fluid flow sim-
86 A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, ulation,” (), https://amplify.fixstars.com/en/demo/fmqa_3_
T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Köpf, aerofoil.
E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, 107 J. Lin, T. Tada, A. Koizumi, M. Sumita, K. Tsuda, and R. Tamura, “Deter-

L. Fang, J. Bai, and S. Chintala, “Pytorch: an imperative style, high- mination of stable proton configurations by black-box optimization using
performance deep learning library,” in Proceedings of the 33rd Interna- an Ising machine,” The Journal of Physical Chemistry C 129, 2332–2340
tional Conference on Neural Information Processing Systems (Curran As- (2025).
sociates Inc., Red Hook, NY, USA, 2019). 108 T. Inoue, Y. Seki, S. Tanaka, N. Togawa, K. Ishizaki, and S. Noda, “To-
87 “Amplify-BBOpt,” (), https://amplify.fixstars.com/en/docs/ wards optimization of photonic-crystal surface-emitting lasers via quan-
amplify-bbopt/v0/. tum annealing,” Opt. Express 30, 43503–43512 (2022).
88 S. Kim, S.-J. Park, S. Moon, Q. Zhang, S. Hwang, S.-K. Kim, T. Luo, 109 T. Kondo, T. Kohira, and Y. Minamoto, “Simultaneous structure design

and E. Lee, “Quantum annealing-aided design of an ultrathin-metamaterial optimization of multiple car models using FMQA,” Transactions of the
optical diode,” Nano Convergence 11 (2024). Society of Automotive Engineers of Japan 56, 229–236 (2025).
89 S. Kim, S. Wu, R. Jian, G. Xiong, and T. Luo, “Design of a high- 110 I. D. W. Samuel and G. A. Turnbull, “Organic semiconductor lasers,”

performance titanium nitride metastructure-based solar absorber using Chemical Reviews 107, 1272–1295 (2007).
quantum computing-assisted optimization,” ACS Applied Materials & In- 111 J. Kennedy and R. Eberhart, “Particle swarm optimization,” in Proceedings

terfaces 15, 40606–40613 (2023). of ICNN’95 - International Conference on Neural Networks, Vol. 4 (1995)
90 S. Kim, S. Jung, A. Bobbitt, E. Lee, and T. Luo, “Wide-angle spectral fil- pp. 1942–1948 vol.4.
ter for energy-saving windows designed by quantum annealing-enhanced 112 D. Wang, D. Tan, and L. Liu, “Particle swarm optimization algorithm: an

active learning,” Cell Reports Physical Science 5, 101847 (2024). overview,” Soft Computing 22, 387–408 (2018).
91 J. Guo, K. Kitai, H. Jippo, and J. Shiomi, “Boosting the quality factor of 113 D. Whitley, “A genetic algorithm tutorial,” Statistics and Computing 4,

tamm structures to millions by quantum inspired classical annealer with 65–85 (1994).
factorization machine,” arXiv:2408.05799 (2024). 114 K. Suzuki, K. Minami, and T. Inamuro, “Lift and thrust generation by
92 M. Urushihara, M. Karube, K. Yamaguchi, and R. Tamura, “Optimization a butterfly-like flapping wing–body model: immersed boundary–lattice
of core–shell nanoparticles using a combination of machine learning and boltzmann simulations,” Journal of Fluid Mechanics 767, 659–695 (2015).
ising machine,” Advanced Photonics Research 4, 2300226 (2023). 115 M. Bando, K. Hasebe, A. Nakayama, A. Shibata, and Y. Sugiyama,
93 R. Tamura, K. Nagata, K. Sodeyama, K. Nakamura, T. Tokuhira, S. Shi- “Structure stability of congestion in traffic dynamics,” Japan Journal of
bata, K. Hammura, H. Sugisawa, M. Kawamura, T. Tsurimoto, M. Naito, Industrial and Applied Mathematics 11, 203–223 (1994).
M. Demura, and T. Nakanishi, “Machine learning prediction of the 116 M. Bando, K. Hasebe, A. Nakayama, A. Shibata, and Y. Sugiyama, “Dy-

mechanical properties of injection-molded polypropylene through X-ray namical model of traffic congestion and numerical simulation,” Physical
diffraction analysis,” Science and Technology of Advanced Materials 25, Review E 51, 1035–1042 (1995).
2388016 (2024).
Black-box optimization using factorization and Ising machines 23

117 T. Kohira, H. Kemmotsu, O. Akira, and T. Tatsukawa, “Proposal of bench- 130 G. Schneider, “Automating drug discovery,” Nature Reviews Drug Discov-

mark problem based on real-world car structure design optimization,” Pro- ery 17, 97–113 (2018).
ceedings of the Genetic and Evolutionary Computation Conference Com- 131 B. Burger, P. M. Maffettone, V. V. Gusev, C. M. Aitchison, Y. Bai,
panion , 183–184 (2018). X. Wang, X. Li, B. M. Alston, B. Li, R. Clowes, N. Rankin, B. Harris,
118 S. Ootomo, T. Harada, and R. Thawonmas, “Proposal of optimization R. S. Sprick, and A. I. Cooper, “A mobile robotic chemist,” Nature 583,
method using common parts information and virtual parent in simultane- 237–241 (2020).
ous design optimization problem of multiple car structures,” Transaction 132 B. P. MacLeod, F. G. L. Parlane, T. D. Morrissey, F. Häse, L. M. Roch,

of the Japanese Society for Evolutionary Computation 9, 41–52 (2018). K. E. Dettelbach, R. Moreira, L. P. E. Yunker, M. B. Rooney, J. R. Deeth,
119 A. Oyama, “Report of evolutionary computation competition 2017,” V. Lai, G. J. Ng, H. Situ, R. H. Zhang, M. S. Elliott, T. H. Haley, D. J. Dvo-
Transaction of the Japanese Society for Evolutionary Computation 9, 86– rak, A. Aspuru-Guzik, J. E. Hein, and C. P. Berlinguette, “Self-driving
92 (2018). laboratory for accelerated discovery of thin-film materials,” Science Ad-
120 S. Daulton, D. Eriksson, M. Balandat, and E. Bakshy, “Multi-objective vances 6, eaaz8867 (2020).
Bayesian optimization over high-dimensional search spaces,” Proceedings 133 R. Tamura, K. Tsuda, and S. Matsuda, “NIMS-OS: an automation soft-

of Machine Learning Research 180, 507–517 (2022). ware to implement a closed loop between artificial intelligence and robotic
121 B. A. Wilson, Z. A. Kudyshev, A. V. Kildishev, S. Kais, V. M. Shalaev, experiments in materials science,” Science and Technology of Advanced
and A. Boltasseva, “Metasurface design optimization via D-Wave based Materials: Methods 3, 2232297 (2023).
sampling,” in Conference on Lasers and Electro-Optics (Optica Publishing 134 G. Tom, S. P. Schmid, S. G. Baird, Y. Cao, K. Darvish, H. Hao, S. Lo,

Group, 2021) p. FTh2M.2. S. Pablo-García, E. M. Rajaonson, M. Skreta, N. Yoshikawa, S. Corapi,

122 B. A. Wilson, Z. A. Kudyshev, A. V. Kildishev, S. Kais, V. M. Shalaev, G. D. Akkoc, F. Strieth-Kalthoff, M. Seifrid, and A. Aspuru-Guzik, “Self-
and A. Boltasseva, “Machine learning framework for quantum sampling of driving laboratories for chemistry and materials science,” Chemical Re-
highly constrained, continuous optimization problems,” Applied Physics views 124, 9633–9732 (2024).
Reviews 8, 041418 (2021). 135 N. Yoshikawa, Y. Asano, D. N. Futaba, K. Harada, T. Hitosugi, G. N.
123 Z. Mao, Y. Matsuda, R. Tamura, and K. Tsuda, “Chemical design with Kanda, S. Matsuda, Y. Nagata, K. Nagato, M. Naito, T. Natsume,
GPU-based Ising machines,” Digital Discovery 2, 1098–1103 (2023). K. Nishio, K. Ono, H. Ozaki, W. Shin, J. Shiomi, K. Shizume, K. Taka-
124 D. Weininger, “SMILES, a chemical language and information system. hashi, S. Takeda, I. Takeuchi, R. Tamura, K. Tsuda, and Y. Ushiku, “Self-
1. introduction to methodology and encoding rules,” Journal of Chemical driving laboratories in Japan,” Digital Discovery 4, 1384–1403 (2025).
Information and Computer Sciences 28, 31–36 (1988). 136 J. Gonzalez, Z. Dai, P. Hennig, and N. Lawrence, “Batch Bayesian opti-
125 W. Jin, R. Barzilay, and T. Jaakkola, “Junction tree variational autoen-
mization via local penalization,” in Proceedings of the 19th International
coder for molecular graph generation,” in Proceedings of the 35th Interna- Conference on Artificial Intelligence and Statistics, Proceedings of Ma-
tional Conference on Machine Learning, Proceedings of Machine Learn- chine Learning Research, Vol. 51, edited by A. Gretton and C. C. Robert
ing Research, Vol. 80, edited by J. Dy and A. Krause (PMLR, 2018) pp. (PMLR, Cadiz, Spain, 2016) pp. 648–657.
2323–2332. 137 J. Azimi, A. Fern, and X. Fern, “Batch Bayesian optimization via simu-
126 E. Jang, S. Gu, and B. Poole, “Categorical reparameterization with
lation matching,” in Advances in Neural Information Processing Systems,
Gumbel-softmax,” (2017), arXiv:1611.01144 [stat.ML]. Vol. 23, edited by J. Lafferty, C. Williams, J. Shawe-Taylor, R. Zemel, and
127 R. Baynazarov and I. Piontkovskaya, “Binary autoencoder for text model-
A. Culotta (Curran Associates, Inc., 2010).
ing,” in Artificial Intelligence and Natural Language, edited by D. Ustalov, 138 Z. Dai, Q. P. Nguyen, S. Tay, D. Urano, R. Leong, B. K. H. Low, and
A. Filchenkov, and L. Pivovarova (Springer International Publishing, P. Jaillet, “Batch Bayesian optimization for replicable experimental de-
Cham, 2019) pp. 139–150. sign,” in Advances in Neural Information Processing Systems, Vol. 36,
128 Q. Gao, G. O. Jones, T. Kobayashi, M. Sugawara, H. Yamashita,
edited by A. Oh, T. Naumann, A. Globerson, K. Saenko, M. Hardt, and
H. Kawaguchi, S. Tanaka, and N. Yamamoto, “Quantum-classical com- S. Levine (Curran Associates, Inc., 2023) pp. 36476–36506.
putational molecular design of deuterated high-efficiency OLED emitters,” 139 A. Glos, A. Krawiec, and Z. Zimborás, “Space-efficient binary optimiza-
Intelligent Computing 2, 0037 (2023). tion for variational quantum computing,” npj Quantum Information 8, 39
129 “IBMQ,” https://www.ibm.com/quantum/technology.
(2022).
140 K. Domino, A. Kundu, Ö. Salehi, and K. Krawiec, “Quadratic and higher-

order unconstrained binary optimization of railway rescheduling for quan-

tum computing,” Quantum Information Processing 21, 337 (2022).