# Black-Box Optimization of Model Superconducting Materials¶

To illustrate the effective use of black-box optimization, this sample code describes the optimization for superconducting materials composed of pseudo-materials as an example problem. Although the present sample code performs material searches based on nonlinear algebraic models, you can perform black-box optimization with the same steps based on high-precision simulations or experimental measurement results instead of model algebraic expressions. Even in such cases, you can use this example almost as is.

For a basic introduction to black-box optimization and FMQA, see "Black-Box Optimization with Quantum Annealing and Ising Machines".

Also, you can find example programs of FMQA applications in various fields such as material search, fluid engineering, chemical plant and urban transportation.

## Problem setting¶

### Search scenario for superconducting materials¶

Superconductivity technology is expected to be utilized in the fields of transportation, such as maglev trains, metrology, and energy. Various superconducting materials are currently being developed to realize superconductivity.

The temperature at which superconductivity is achieved (critical temperature) is generally around the absolute temperature of 0 K (Kelvin) for currently confirmed superconducting materials. Because of this, superconductivity requires costly cooling to be exploited, and its application in the real world is currently limited. Therefore, the search for high-temperature superconductors is a pressing issue.

Typically, the search for materials that realize superconductivity involves a trial-and-error process of selecting and synthesizing several materials, repeatedly evaluating the critical temperature of the synthesized materials by measurement, and identifying the material to be synthesized that achieves a higher critical temperature. This process of synthesis and critical temperature evaluation is considered time-consuming. For this search, a black box optimization method is used to find a combination of materials close to the optimal solution with a relatively small number of evaluations.

In this example, the search for superconducting materials consisting of pseudo materials is treated as an example to illustrate the material search by a black-box optimization method (FMQA), and a critical temperature model is used to evaluate the critical temperature. Note that the critical temperature models presented below and the combinations of materials obtained are not necessarily physically accurate, and thus the example serves for illustration purposes only.

## Solver client configuration¶

First, we will configure the solver client of the solver we will use during the optimization cycles. In this example program, we will use Amplify AE.

```
from amplify import FixstarsClient
from datetime import timedelta
# Set a solver client
client = FixstarsClient()
# If you use Amplify in a local environment, enter the Amplify API token.
# client.token = "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx"
# Set timeout to be 2 seconds
client.parameters.timeout = timedelta(milliseconds=2000)
```

## Definition of black-box function¶

This example code selects a combination of several materials from $D$ types of materials and performs an optimization to maximize the critical temperature of the superconducting material produced by their synthesis.

In general, the critical temperature can be evaluated by experimental measurement, which requires a
relatively large cost (time and money) each time the evaluation is performed. In this example code,
instead of measuring the critical temperature, the following critical temperature model
`supercon_temperature()`

is used for evaluation. However, this function is only a
substitute
for experimental measurement, and its contents and parameters are treated as unknown, and the number
of
calls to `supercon_temperature()`

is also treated as limited.

The following `make_blackbox_func`

function creates and returns a blackbox function
`blackbox`

, which is also the objective function in the present FMQA. The blackbox function
`blackbox`

executes `supercon_temperature()`

and returns the negative value of
the
obtained critical temperature.

In this example, optimization will progress to minimize the negative value of the obtained critical temperature.

```
import numpy as np
import math
from typing import Callable, Any
# Set the random seed
seed = 1234
rng = np.random.default_rng(seed)
def make_blackbox_func(d: int) -> Callable[[np.ndarray], float]:
"""Returns a function that takes a binary vector with size d as input and returns a float value"""
def set_properties(size: int) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
"""Returns randomly chosen material properties"""
mu, sigma, ratio = 0.0, 1.0, 0.2
table1 = rng.random(size) * 1e5 * (0.1 * math.log(size) - 0.23)
table2 = rng.lognormal(mu, sigma, size) * ratio
table3 = rng.lognormal(mu, sigma, size) * ratio
return table1, table2, table3
def supercon_temperature(
x: np.ndarray,
debye_table: np.ndarray,
state_table: np.ndarray,
interaction_table: np.ndarray,
) -> float:
"""For a given material choice (a binary variable vector of size d), compute and return the critical temperature. (substitute of simulation or experiment)"""
debye_temperature = np.sum(x * debye_table) / np.sum(x)
state_density = np.sum(x * state_table) / np.sum(x)
interaction = np.sum(x * interaction_table) / np.sum(x)
crit_temp = debye_temperature * math.exp(-1.0 / state_density / interaction)
return crit_temp
# Prepare property tables
debye_temperature_table, state_density_table, interaction_table = set_properties(d)
# Definition of the black-box function
def blackbox(x: np.ndarray) -> float:
"""For a given material choice (a binary vector of size d), returns the negative of the critical temperature."""
assert x.shape == (d,) # x is the 1D binary vector of size d
t_c = supercon_temperature(
x, debye_temperature_table, state_density_table, interaction_table
)
return -t_c
return blackbox
```

In the following, the black-box function `blackbox(x)`

for the critical temperature
defined
above is used to evaluate the critical temperature of superconducting materials synthesized from a
random
selection of materials. Here, `num_materials`

is the number of materials to be selected,
and
the input `x`

is a vector of size `num_materials`

consisting of 0 or 1.

For example, in the case of selecting the first and last of five materials to be combined, the input
vector would be `x = [1, 0, 0, 0, 0, 1]`

. In this case, there are $2^5-1=31$ possible
choices
(combinations).

For `num_materials = 100`

, the number of combinations is approximately $10^{30}$, and the
full-search method is considered difficult.

```
num_materials = (
100 # The size of the binary decision variable (number of materials to choose from)
)
blackbox_func = make_blackbox_func(num_materials)
# Evaluate the black-box function with the random inputs for n_cycle times, and print the minimum and average values of the blackbox function
n_cycle = 100
obj_min = 0.0 # The variable to save the minimum of the negative critical temperature
obj_mean = 0.0 # The variable to save the average of the negative critical temperature
for i in range(n_cycle):
x = rng.integers(0, 2, num_materials)
if np.sum(x) == 0:
continue
obj = blackbox_func(x)
if obj_min > obj:
obj_min = obj
obj_mean += obj
obj_mean /= n_cycle
print(f"Minimum objective function value: {obj_min:.2f} K")
print(f"Mean objective function value: {obj_mean:.2f} K")
```

## FMQA program implementation¶

### Model training by machine learning¶

This section implements the part of FMQA that learns the optimal parameters (weights and bias) of the
model by machine learning. First, the `TorchFM`

class representing the model by the
Factorization Machine is defined using PyTorch.

The following equation represents the Factorization Machine.

$$ \begin{aligned} f(\boldsymbol{x} | \boldsymbol{w}, \boldsymbol{v}) &= \underset{\color{red}{\mathtt{out\_linear}}}{\underline{ w_0 + \sum_{i=1}^d w_i x_i} } + \underset{\color{red}{\mathtt{out\_quadratic}}}{\underline{\frac{1}{2} \left[\underset{\color{red}{\mathtt{out\_1}}}{\underline{ \sum_{f=1}^k\left(\sum_{i=1}^d v_{i f} x_i\right)^2 }} - \underset{\color{red}{\mathtt{out\_2}}}{\underline{ \sum_{f=1}^k\sum_{i=1}^d v_{i f}^2 x_i^2 }} \right] }} \end{aligned} $$

The input $x$ of this model is a vector of the length $d$ as the input to the black-box function, with the following three parameters.

- $v$: 2-dimensional array of $d\times k$.
- $w$: 1D vector of length $d$.
- $w_0$: scalar

The only hyperparameter is $k$, which is given as a positive integer less than or equal to $d$.

The `TorchFM`

class defined below inherits from `torch.nn.Module`

and is
constructed from an input vector $x$ of size $d$ and a hyperparameter $k$. The hyperparameter $k$
controls
the number of parameters in the model; the larger the hyperparameter, the more parameters are, the
more
accurate the model becomes, but also the more prone to over-fitting.

The `TorchFM`

class has attributes $v$, $w$, and $w_0$ of the model parameters and updates
these parameters as the training proceeds. According to the above formula, the `forward`

method
also outputs an estimate of $y$ from the input $x$. Since the parameters $v$, $w$, and $w_0$ are
needed to
construct a QUBO model later, we also define a function `get_parameters`

to output these
parameters.

```
import torch
import torch.nn as nn
# Set a random seed
torch.manual_seed(seed)
class TorchFM(nn.Module):
def __init__(self, d: int, k: int):
"""Define a model
Args:
d (int): Size of input vector
k (int): Hyperparameter k
"""
super().__init__()
self.d = d
self.v = torch.randn((d, k), requires_grad=True)
self.w = torch.randn((d,), requires_grad=True)
self.w0 = torch.randn((), requires_grad=True)
def forward(self, x: torch.Tensor) -> torch.Tensor:
"""Take x as input and returns a predicted value y
Args:
x (torch.Tensor): 2D tensor of (number of samples × d)
Returns:
torch.Tensor: 1D tensor of predicted y (size is the number of samples)
"""
out_linear = torch.matmul(x, self.w) + self.w0
out_1 = torch.matmul(x, self.v).pow(2).sum(1)
out_2 = torch.matmul(x.pow(2), self.v.pow(2)).sum(1)
out_quadratic = 0.5 * (out_1 - out_2)
out = out_linear + out_quadratic
return out
def get_parameters(self) -> tuple[np.ndarray, np.ndarray, float]:
"""Returns the parameters (weights and bias) v, w, w0"""
np_v = self.v.detach().numpy().copy()
np_w = self.w.detach().numpy().copy()
np_w0 = self.w0.detach().numpy().copy()
return np_v, np_w, float(np_w0)
```

Next, we define a function `train()`

to perform machine learning of the model based on the
`TorchFM`

class. The input is the training data $x, y$ and an instance of the
`TorchFM`

model. By calling the `train()`

function, the `TorchFM`

parameters are trained and optimized.

In general machine learning, we split the data into training and validation data, optimize the parameters using the training data, and validate the model during training using the validation data. The model is validated at each epoch, and the parameters at the epoch with the highest prediction accuracy for the validation data are saved and used as the model after learning.

```
from torch.utils.data import TensorDataset, DataLoader, random_split
from tqdm.auto import tqdm, trange
import copy
def train(
x: np.ndarray,
y: np.ndarray,
model: TorchFM,
) -> None:
"""Perform training of a FM model
Args:
x (np.ndarray): Training data (input vectors)
y (np.ndarray): Training data (output values)
model (TorchFM): a TorchFM model
"""
# Number of iterations
epochs = 2000
# Model optimizer
optimizer = torch.optim.AdamW([model.v, model.w, model.w0], lr=0.1)
# Loss function
loss_func = nn.MSELoss()
# Prepare the dataset
x_tensor, y_tensor = (
torch.from_numpy(x).float(),
torch.from_numpy(y).float(),
)
dataset = TensorDataset(x_tensor, y_tensor)
train_set, valid_set = random_split(dataset, [0.8, 0.2])
train_loader = DataLoader(train_set, batch_size=8, shuffle=True)
valid_loader = DataLoader(valid_set, batch_size=8, shuffle=True)
# Execute the training
min_loss = 1e18 # Save the minimum loss value
best_state = model.state_dict() # Save the best model parameters
# Visualize progress using `tqdm` instead of `range`
for _ in trange(epochs, leave=False):
# Training phase
for x_train, y_train in train_loader:
optimizer.zero_grad()
pred_y = model(x_train)
loss = loss_func(pred_y, y_train)
loss.backward()
optimizer.step()
# Validation phase
with torch.no_grad():
loss = 0
for x_valid, y_valid in valid_loader:
out_valid = model(x_valid)
loss += loss_func(out_valid, y_valid)
if loss < min_loss:
# Save the model parameter is the loss function value is updated
best_state = copy.deepcopy(model.state_dict())
min_loss = loss
# Update model based on the trainded parameters
model.load_state_dict(best_state)
```

## Solving the model with Amplify¶

Next, we implement the `anneal`

function to minimize the inferred machine learning model,
corresponding to the lower left part of the FMQA cycle diagram (reproduced below), where the input is
the
model `TorchFM`

class after training, and the output is a vector $x$ that minimizes the
model.

Solve the following optimization problem for the class `TorchFM`

model learned earlier to
find
the input $x$ such that the inferred model is minimized.

$$ \underset{x}{\mathrm{argmin}} \quad \underset{\color{red}{\mathtt{out\_linear}}}{\underline{ w_0 + \sum_{i=1}^d w_i x_i} } + \underset{\color{red}{\mathtt{out\_quadratic}}}{\underline{\frac{1}{2} \left[\underset{\color{red}{\mathtt{out\_1}}}{\underline{ \sum_{f=1}^k\left(\sum_{i=1}^d v_{i f} x_i\right)^2 }} - \underset{\color{red}{\mathtt{out\_2}}}{\underline{ \sum_{f=1}^k\sum_{i=1}^d v_{i f}^2 x_i^2 }} \right] }} $$

The decision variable in this optimization problem is $x$. It is a one-dimensional binary variable vector of length $d$, just like the input vector to the black-box function. Also, $v$, $w$, and $w_0$, parameters in the learning phase, are constants here.

We define the `anneal`

function that performs the optimization using Amplify AE for the
given
model as follows. The `anneal`

function uses the `VariableGenerator`

to create a
1D
binary vector `x`

of length $d$, and then uses the binary variable array $x$ and $v$, $w$,
and
$w_0$ obtained from the `TorchFM`

class to create a Factorization Machine The objective
function to be optimized is created according to the formula of `Factorization Machine`

.

After building the optimization model, we will perform minimization of the objective function using
the
created QUBO model and the solver client `FixstarsClient`

defined earlier.

```
from amplify import VariableGenerator, Model, solve, Poly
def anneal(torch_model: TorchFM) -> np.ndarray:
"""Take the parameters of the FM model and find x that gives the minimum value of the FM model as the QUBO model described by those parameters."""
# Generate a binary variable array of size d
gen = VariableGenerator()
x = gen.array("Binary", torch_model.d)
# Obtain FM parameters v, w, w0 from TorchFM
v, w, w0 = torch_model.get_parameters()
# Create the objective function
out_linear = w0 + (x * w).sum()
out_1 = ((x[:, np.newaxis] * v).sum(axis=0) ** 2).sum() # type: ignore
out_2 = ((x[:, np.newaxis] * v) ** 2).sum()
objective: Poly = out_linear + (out_1 - out_2) / 2
# Construct combinatorial optimization model
amplify_model = Model(objective)
# Execute the optimization (pass the constructed QUBO model and the solver client created in the beginning)
result = solve(amplify_model, client)
if len(result.solutions) == 0:
raise RuntimeError("No solution was found.")
# Return the input that minimizes the model
return x.evaluate(result.best.values).astype(int)
```

Now that we have defined the `train`

function for machine learning and the
`anneal`

function for optimization, which are the core of FMQA, we can execute FMQA using these functions.

The black-box function (corresponding to experiment or simulation) to be used for th epresent
black-box
optimization is the `blackbox_func()`

already defined above. This function takes a NumPy
one-dimensional binary vector of length $d = 100$ consisting of $0$ or $1$ and returns the negative
value
of the critical temperature.

### Creating initial training data¶

Next, the `black-box`

function $y = f(\boldsymbol{x})$ is evaluated on the input vector
$\boldsymbol{x}$ to create $N_0$ initial training data. The `black-box`

function
corresponds to
the results of a simulation or experiment, so you can create these data based on previous experiments
or
simulations.

We define the `init_training_data`

function to create the initial training data using a
random
$N_0$ input vector $x$ as follows. This function takes the `black-box`

function and the
number
of initial training data $N_0$ and returns $N_0$ input vectors $\boldsymbol{x}$ and corresponding
output
$y$ as initial training data.

```
def init_training_data(d: int, n0: int):
"""Craete initial training data of n0 samples"""
assert n0 < 2**d
# Generate n0 random inputs of size d
x = rng.choice(np.array([0, 1]), size=(n0, d))
# Check if these random inputs are unique. If not, modify corresponding inputs
x = np.unique(x, axis=0)
while x.shape[0] != n0:
x = np.vstack((x, np.random.randint(0, 2, size=(n0 - x.shape[0], d))))
x = np.unique(x, axis=0)
# Evaluate black-bok function to obtain n0 output values
y = np.zeros(n0)
for i in range(n0):
y[i] = blackbox_func(x[i])
return x, y
N0 = 10 # Number of samples in the initlal training data
x_init, y_init = init_training_data(num_materials, N0)
```

Now, we have $N_0$-pair of initial training data (10 pairs of the input vector of size 100 and the output value).

## Running an FMQA cycle¶

Using the $N_0$ pairs of initial training data created above, we will run the FMQA cycles, by performing the following operations per FMQA cycle.

- Train the model
- Construct a model of class
`TorchFM`

and train it by calling the`train`

function on the model with initial training data`x = x_init`

,`y = y_init`

.

- Construct a model of class
- Minimize the model
- Minimize the model by passing a trained model of class
`TorchFM`

to the`anneal`

function to get $\hat{x}$ that minimizes the model. - If $\hat{x}$ is already included in the training data
`x`

, change part of $\hat{x}$ to avoid duplicate training data.

- Minimize the model by passing a trained model of class
- Evaluation of black-box function
- Input $\hat{x}$ to the black-box function and get output $\hat{y}$.
- Add $\hat{x}$ and $\hat{y}$ to the traiing data
`x`

and`y`

, respectively.

An example implementation of the above is as follows: Since the execution takes approximately 10 minutes of computation time, an example output is shown in Execution example.

```
# Number of FMQA cycles
N = 40
# Initialize training data
x, y = x_init, y_init
# Perform N FMQA cycles
# Use `tqdm` instead of `range` to visualize the progress
for i in trange(N):
# Create machine-learning model
model = TorchFM(num_materials, k=10)
# Train the machine-learning model
train(x, y, model)
# Obtain the input that yields the minimum value of the trained model
x_hat = anneal(model)
# Ensure the uniqueness of x_hat
while (x_hat == x).all(axis=1).any():
flip_idx = rng.choice(np.arange(num_materials))
x_hat[flip_idx] = 1 - x_hat[flip_idx]
# Evaluate the black-box function with the estimated x_hat
y_hat = blackbox_func(x_hat)
# Add the evaluation results to the training data set
x = np.vstack((x, x_hat))
y = np.append(y, y_hat)
tqdm.write(f"FMQA cycle {i}: found y = {y_hat}; current best = {np.min(y)}")
```

After running the above cell, `x`

and `y`

have input-output pairs of black-box
function evaluations. The number of pairs is $N_0 + N = 50$.

```
print(x.shape, y.shape)
```

The following gives the input that gives the minimum value of the black box function evaluation and its value.

```
min_idx = np.argmin(y)
print(f"best x = {x[min_idx]}")
print(f"best y = {y[min_idx]}")
```

### Plotting the evolution of the evaluated values¶

Below are $N_0$ initial training data and the evolution of the optimized evaluation values of the black-box function over $N$ FMQA cycles. The initial training data is shown in blue, and the evaluation values obtained by the FMQA cycles are shown in red.

```
import matplotlib.pyplot as plt
fig = plt.figure(figsize=(6, 4))
ax = fig.add_subplot()
# Evaluation results obtained during initial training data generation
ax.plot(
range(N0),
y[:N0],
marker="o",
linestyle="-",
color="b",
)
# Evaluation results during the FMQA cycles
ax.plot(
range(N0, N0 + N),
y[N0:],
marker="o",
linestyle="-",
color="r",
)
ax.set_xlabel("number of iterations", fontsize=18)
ax.set_ylabel("f(x)", fontsize=18)
ax.tick_params(labelsize=18)
plt.show()
```

### Execution example¶

The output of a typical plot obtained by running this sample code is shown below. The optimization uses a QUBO solver, so the results will vary from run to run. However, you can see that the evaluated value of the black-box function decreases with the number of FMQA cycles.