Jupyter Notebookで実行する
サンプルコード本サンプルコードでは、ブラックボックス最適化の一手法である FMQA を紹介するために、ある代数式を未知のブラックボックス関数と見なし、その出力値を最小化するような入力値を推定します。より現実的なモデルケースにおける FMQA の実施例やサンプルコードは、以下のリンクをご覧ください。
また、最適化対象が既知関数で2次式である場合、二次制約なし二値最適化(QUBO: Quadratic Unconstrained Binary Optimization)形式での最適化が可能です。QUBO問題の解説やサンプルコード、Amplify の使い方は、以下のリンクをご覧ください(一部抜粋)。
本ノートブックは、以下の章立ての構成となっています。
FMQA は、ベイズ最適化に似たブラックボックス最適化法の一つです。通常、数理最適化では、何らかの目的関数 $f(\boldsymbol{x})$ を最小化(あるいは最大化)するような決定変数 $\boldsymbol{x}$ を推定することを目的とします。ここで、目的関数 $f(\boldsymbol{x})$ に関する情報(関数形、勾配、劣モジュラ性、凸性等)が与えられている場合、効率的な最適化が可能です。
$$ \begin{aligned} \mathrm{Minimize}&\,\,f(\boldsymbol{x}) \\ \mathrm{subject\,\,to\,\,}&\boldsymbol{x} \in [0,1]^D \end{aligned} $$例えば、Amplify のデモ・チュートリアルで紹介しているいくつかの最適化問題のように、$f(\boldsymbol{x})$ の関数が既知(かつ $\boldsymbol{x}$ の2次式)の場合、$f(\boldsymbol{x})$ を目的関数とすることで、直接、二次制約なし二値最適化(QUBO: Quadratic Unconstrained Binary Optimization)としての最適化実施が可能です。
※ここで、$\boldsymbol{x}$ としてバイナリ変数ベクトルを仮定しますが、非バイナリ変数を one-hot エンコーディング等を使い、バイナリ変数に変換することができます。そのようなサンプルは、ブラックボックス最適化と流体シミュレーションによる翼形状の最適化でご覧になれます。
一方、物理現象や社会現象に対するシミュレーションや実験によって得られる値を最小化(または最大化)する最適化の場合、目的関数 $f(\boldsymbol{x})$ はシミュレーションあるいは実験ということになり、目的関数を具体的な式で記述することはできません。このような未知の目的関数 $f(\boldsymbol{x})$ に対して行う数理最適化のことをブラックボックス最適化と呼びます。また、そのような目的関数の評価(シミュレーションや実験の実施)には、一般的に比較的大きなコストが必要なため、決定変数の集合が有限であっても、全検索による最適化は困難な場合が多く、できるだけ少ない目的関数の評価回数での最適化が要求されます。
ベイズ最適化では、以下の最適化サイクルを繰り返すことでブラックボックス最適化を実施します。
このサイクルの繰り返しに伴い、最適化点近傍における獲得関数 $g(\boldsymbol{x})$ の予測精度が向上し、その結果、得られる $\hat{\boldsymbol{x}}$ は、目的関数 $f(\boldsymbol{x})$ を最小化する真の決定変数に近い値を取ることが期待されます。一方で、このベイズ最適化のサイクルにおいては、次の2つの課題があります。
ベイズ最適化で重要となるこれら2つの課題を解決し、ブラックボックス最適化を実現する汎用的な手法として、次に説明する FMQA があります。
ベイズ最適化で必要な獲得関数 $g(\boldsymbol{x})$ として、次式のような機械学習モデルの一種である Factorization Machine (FM) を用いる場合を考えます。
$$ \begin{aligned} g(\boldsymbol{x} | \boldsymbol{w}, \boldsymbol{v}) &= w_0 + \langle \boldsymbol{w}, \boldsymbol{x}\rangle + \sum_{i=1}^D \sum_{j=i+1}^D \langle \boldsymbol{v}_i, \boldsymbol{v}_j \rangle x_i x_j \\ &=w_0 + \sum_{i=1}^D w_i x_i + \sum_{i=1}^D \sum_{j=i+1}^D \sum_{f=1}^k v_{if}v_{jf}x_ix_j \\ &=w_0 + \sum_{i=1}^D w_i x_i + \frac{1}{2}\sum_{f=1}^k\left(\left(\sum_{i=1}^D v_{i f} x_i\right)^2 - \sum_{i=1}^D v_{i f}^2 x_i^2\right) \end{aligned} $$ここで、FM は $\boldsymbol{x}$ の2次式であるため、上式は QUBO による最適化が可能な関数形となります。式中の $\boldsymbol{w}$ や $\boldsymbol{v}$($v_{ij}$, $w_i$)は、上式のモデルを機械学習した後に得られる FM パラメータ(機械学習における重みやバイアス)であり、$k$ はハイパーパラメータです。
FM パラメータ数は、ハイパーパラメータである $k$ に依存し、$k=D$ のとき、FM は QUBO の相互作用項と同じ自由度がある一方、$k$ を小さくすることで FM パラメータ数を減らし過学習を抑制する効果があります。
獲得関数 $g(\boldsymbol{x})$ に FM を採用し、その最適化を量子アニーリング(QA)やイジングマシンを用いて実施することで、ベイズ最適化における前述の課題を解決し、一般的な問題に適用することができます。このように、量子アニーリング・イジングマシンと機械学習を融合して行うブラックボックス最適化手法を FMQA と呼ぶ場合があります(4. 参考文献 参照)。
FMQA の実施手順は、上述ベイズ最適化のサイクルと同様ですが、次のように進めます。
まず、目的関数評価に掛かる時間・金銭コストから、最適化において実施可能な目的関数評価の回数 $N$ を見積ります。例えば、一度の目的関数評価(実験又はシミュレーション)に1時間かかり、FMQA による最適化を1日で終える必要がある場合、最大でも $N=24$ であると考えられます。そして、$N_0<N$ であるような初期教師データのサンプル数 $N_0$ を決定し、以下の通り、初期教師データを準備、そして、FMQA サイクルを $N-N_0$ 回実施します。
初期教師データの事前準備 ($N_0$ セット)
FMQA による最適化サイクルの実施 ($N-N_0$ 回)
上記1~4を $N-N_0$ 回繰り返す。
上記 FMQA サイクルの繰り返しと共に、最適化点近傍における FM の予測精度が向上し、量子アニーリング・イジングマシンによるより良い $\hat{\boldsymbol{x}}$ の推定が期待されます。
import os
import torch
import numpy as np
def seed_everything(seed=0):
os.environ["PYTHONHASHSEED"] = str(seed)
np.random.seed(seed)
torch.manual_seed(seed)
torch.cuda.manual_seed(seed)
torch.backends.cudnn.deterministic = True
from amplify.client import FixstarsClient
client = FixstarsClient()
client.parameters.timeout = 1000 # タイムアウト1秒
# client.token = "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx" # ローカル環境等で使用する場合は、Amplify AEのアクセストークンを入力してください。
FM の学習と推論を PyTorch で行います。TorchFM
クラスでは、機械学習モデルとしての $g(\boldsymbol{x})$
を定義します。下式の通り、$g(\boldsymbol{x})$ 内の各項は、TorchFM
クラス内の
out_lin
、out_1
、out_2
、out_inter
に直接対応します。
import torch.nn as nn
class TorchFM(nn.Module):
def __init__(self, d: int, k: int):
super().__init__()
self.V = nn.Parameter(torch.randn(d, k), requires_grad=True)
self.lin = nn.Linear(d, 1) # 右辺第1項及び2項は全結合ネットワーク
def forward(self, x):
out_1 = torch.matmul(x, self.V).pow(2).sum(1, keepdim=True)
out_2 = torch.matmul(x.pow(2), self.V.pow(2)).sum(1, keepdim=True)
out_inter = 0.5 * (out_1 - out_2)
out_lin = self.lin(x)
out = out_inter + out_lin
return out
次に、入出力データから FM を機械学習する関数 train()
を定義します。一般的な機械学習と同様に、教師データを学習データと検証データに分割し、学習データを用いてパラメータの最適化、検証データを用いて学習中のモデル検証を行います。train()
関数は、検証データに対して最も予測精度の高かったモデルを返します。
from torch.utils.data import TensorDataset, DataLoader
from sklearn.model_selection import train_test_split
import copy
def train(
X,
y,
model_class=None,
model_params=None,
batch_size=1024,
epochs=3000,
criterion=None,
optimizer_class=None,
opt_params=None,
lr_sche_class=None,
lr_sche_params=None,
):
X_tensor, y_tensor = (
torch.from_numpy(X).float(),
torch.from_numpy(y).float(),
)
indices = np.array(range(X.shape[0]))
indices_train, indices_valid = train_test_split(
indices, test_size=0.2, random_state=42
)
train_set = TensorDataset(X_tensor[indices_train], y_tensor[indices_train])
valid_set = TensorDataset(X_tensor[indices_valid], y_tensor[indices_valid])
loaders = {
"train": DataLoader(train_set, batch_size=batch_size, shuffle=True),
"valid": DataLoader(valid_set, batch_size=batch_size, shuffle=False),
}
model = model_class(**model_params)
best_model_wts = copy.deepcopy(model.state_dict())
optimizer = optimizer_class(model.parameters(), **opt_params)
if lr_sche_class is not None:
scheduler = lr_sche_class(optimizer, **lr_sche_params)
best_score = 1e18
for epoch in range(epochs):
losses = {"train": 0.0, "valid": 0.0}
for phase in ["train", "valid"]:
if phase == "train":
model.train()
else:
model.eval()
for batch_x, batch_y in loaders[phase]:
optimizer.zero_grad()
out = model(batch_x).T[0]
loss = criterion(out, batch_y)
losses[phase] += loss.item() * batch_x.size(0)
with torch.set_grad_enabled(phase == "train"):
if phase == "train":
loss.backward()
optimizer.step()
losses[phase] /= len(loaders[phase].dataset)
with torch.no_grad():
model.eval()
if best_score > losses["valid"]:
best_model_wts = copy.deepcopy(model.state_dict())
best_score = losses["valid"]
if lr_sche_class is not None:
scheduler.step()
with torch.no_grad():
model.load_state_dict(best_model_wts)
model.eval()
return model
入力値 $\boldsymbol{x}$ に対して目的関数 $f(\boldsymbol{x})$ を評価し、$N_0$ 個の入出力ペア(初期教師データ)を作成します。ここでの入力値 $\boldsymbol{x}$ の決め方は様々ですが、乱数を用いたり、現象に対する知見に基づき機械学習に適した値を用いたりします。過去に実施した実験やシミュレーションの結果から、教師データを構築しても構いません。
def gen_training_data(D: int, N0: int, true_func):
assert N0 < 2**D
# N0個の入力値を乱数を用いて取得
X = np.random.randint(0, 2, size=(N0, D))
# 取得した入力値のうち重複しているものを除外し、除外した分の入力値を乱数を用いて追加
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)
y = np.zeros(N0)
# N0個の入力値に対応する出力値を目的関数を評価して取得
for i in range(N0):
if i % 10 == 0:
print(f"Generating {i}-th training data set.")
y[i] = true_func(X[i])
return X, y
FMQA.cycle()
では、事前に準備した初期教師データを用い、FMQA サイクルを $N-N_0$ 回実施します。FMQA.step()
は、FMQA
を1サイクルのみ行う関数で、FMQA.cycle()
から $N-N_0$ 回呼び出されます。
from amplify import (
Solver,
BinarySymbolGenerator,
sum_poly,
BinaryMatrix,
BinaryQuadraticModel,
)
import matplotlib.pyplot as plt
import sys
class FMQA:
def __init__(self, D: int, N: int, N0: int, k: int, true_func, solver) -> None:
assert N0 < N
self.D = D
self.N = N
self.N0 = N0
self.k = k
self.true_func = true_func
self.solver = solver
self.y = None
# 教師データに基づいて N-N0 回のFMQAを教師データを追加しながら繰り返し実施するメンバー関数
def cycle(self, X, y, log=False) -> np.ndarray:
print(f"Starting FMQA cycles...")
pred_x = X[0]
pred_y = 1e18
for i in range(self.N - self.N0):
print(f"FMQA Cycle #{i} ", end="")
try:
x_hat = self.step(X, y)
except RuntimeError:
sys.exit(f"Unknown error, i = {i}")
# x_hat として既に全く同じ入力が教師データ内に存在する場合、その周辺の値を x_hat とする。
is_identical = True
while is_identical:
is_identical = False
for j in range(i + self.N0):
if np.all(x_hat == X[j, :]):
change_id = np.random.randint(0, self.D, 1)
x_hat[change_id.item()] = 1 - x_hat[change_id.item()]
if log:
print(f"{i=}, Identical x is found, {x_hat=}")
is_identical = True
break
# hat{x} で目的関数 f() を評価
y_hat = self.true_func(x_hat)
# 最適点近傍における入出力ペア [x_hat, y_hat] を教師データに追加
X = np.vstack((X, x_hat))
y = np.append(y, y_hat)
# 目的関数の評価値が最小値を更新したら、その入出力ペアを [pred_x, pred_y] へコピー
if pred_y > y_hat:
pred_y = y_hat
pred_x = x_hat
print(f"variable updated, {pred_y=}")
else:
print("")
# 全ての入力を全探索済みの場合は、for文を抜ける
if len(y) >= 2**self.D:
print(f"Fully searched at {i=}. Terminating FMQA cycles.")
break
self.y = y
return pred_x
# 1回のFMQAを実施するメンバー関数
def step(self, X, y) -> np.ndarray:
# FM を機械学習
model = train(
X,
y,
model_class=TorchFM,
model_params={"d": self.D, "k": self.k},
batch_size=8,
epochs=2000,
criterion=nn.MSELoss(),
optimizer_class=torch.optim.AdamW,
opt_params={"lr": 1},
)
# 学習済みモデルから、FM パラメータの抽出
v, w, w0 = list(model.parameters())
v = v.detach().numpy()
w = w.detach().numpy()[0]
w0 = w0.detach().numpy()[0]
# ここから量子アニーリング・イジングマシンによる求解を実施
gen = BinarySymbolGenerator() # BinaryPoly の変数ジェネレータを宣言
q = gen.array(self.D) # BinaryPoly から決定変数の作成
cost = self.__FM_as_QUBO(q, w0, w, v) # FM パラメータから QUBO として FM を定義
result = self.solver.solve(cost) # 目的関数を Amplify のソルバーに受け渡し
if len(result.solutions) == 0:
raise RuntimeError("No solution was found.")
values = result.solutions[0].values
q_values = q.decode(values)
return q_values
# FM パラメータから QUBO として FM を定義する関数。前定義の TorchFM クラスと同様に、g(x) の関数形通りに数式を記述。
def __FM_as_QUBO(self, x, w0, w, v):
lin = w0 + (x.T @ w)
D = w.shape[0]
out_1 = sum_poly(self.k, lambda i: sum_poly(D, lambda j: x[j] * v[j, i]) ** 2)
# 次式において、x[j] はバイナリ変数なので、x[j] = x[j]^2 であることに注意。
out_2 = sum_poly(
self.k, lambda i: sum_poly(D, lambda j: x[j] * v[j, i] * v[j, i])
)
return lin + (out_1 - out_2) / 2
"""上記の __FM_as_QUBO で用いられている sum_poly は、計算速度やメモリの観点から非効率。一般的に決定変数の相互作用項が非ゼロである FM の場合、BinaryMatrix を使う次の書き方が効率的。ここで、BinaryMatrixでの2次項は、上三角行列で表される非対角項に対応するため、FM式の2次の項に対する x(1/2) は不要。また、上の __FM_as_QUBO(sum_poly を使う実装)と関数のシグネチャを合わせるために、x を引数に取っているが、BinaryMatrix を使う本実装では本来は不要。
def __FM_as_QUBO(self, x, w0, w, v):
out_1_matrix = v @ v.T
out_2_matrix = np.diag((v * v).sum(axis=1))
matrix = BinaryMatrix(out_1_matrix - out_2_matrix + np.diag(w))
# 定数項 w0 を忘れずに BinaryQuadraticModel の2つ目の引数に入れる。
model = BinaryQuadraticModel(matrix, w0)
return model
"""
# 初期教師データ及び各 FMQA サイクル内で実施した i 回の目的関数評価値の履歴をプロットする関数
def plot_history(self):
assert self.y is not None
fig = plt.figure(figsize=(6, 4))
plt.plot(
[i for i in range(self.N0)],
self.y[: self.N0],
marker="o",
linestyle="-",
color="b",
) # 初期教師データ生成時の目的関数評価値(ランダム過程)
plt.plot(
[i for i in range(self.N0, self.N)],
self.y[self.N0 :],
marker="o",
linestyle="-",
color="r",
) # FMQA サイクル時の目的関数評価値(FMQA サイクル過程)
plt.xlabel("i-th evaluation of f(x)", fontsize=18)
plt.ylabel("f(x)", fontsize=18)
plt.tick_params(labelsize=18)
return fig
本サンプルコードで最適化対象とする代数式は、次の2次式です。$Q$ は乱数で生成される成分の平均値がゼロな $d$ 次元対称行列であり、make_Q
で作成されます。
本来 FMQA は、未知関数に対して実施されるものであることに注意してください。今回は、簡単な理解のため、既知関数である上記 $f(\boldsymbol{x})$ を、未知関数(=ブラックボックス)として取り扱っています。
以下の条件($D=100$、$N=100$、$N_0=70$)では、FMQA のサイクル完了までにおよそ数分の計算時間を要しますので、ご注意下さい。出力例を、『3.3. FMQA サンプルコード実行例』に紹介しています。
# d次元対称行列であって成分の平均が0のものを出力
def make_Q(d) -> np.ndarray:
Q_true = np.random.rand(d, d)
Q_true = (Q_true + Q_true.T) / 2
Q_true = Q_true - np.mean(Q_true)
return Q_true
# 乱数シード値を初期化
seed_everything(0)
# 入力値の次元(問題サイズ)
D = 100
# 真の関数で使われる行列 Q
Q = make_Q(D)
# 目的関数(xQx)の定義。
def true_func(x):
# 本来は、未知関数(シミュレーションや実験)の結果値やその後処理結果値を cost とする。
cost = x @ Q @ x
return cost
N = 70 # 関数を評価できる回数
N0 = 60 # 初期教師データのサンプル数
k = 10 # FMにおけるベクトルの次元(ハイパーパラメータ)
# client:先に作成した Amplify クライアント
solver = Solver(client)
# 初期教師データの生成
X, y = gen_training_data(D, N0, true_func)
# FMQA のインスタンス化
fmqa_solver = FMQA(D, N, N0, k, true_func, solver)
# FMQA サイクルの実行
pred_x = fmqa_solver.cycle(X, y)
# 最適化結果の出力
print("pred x:", pred_x)
print("pred value:", true_func(pred_x))
初期教師データ作成時にランダムに生成した入力値に対して得られた $N_0$ 個の目標関数値及び $N-N_0$ サイクルの FMQA 最適化過程における目標関数値の推移を以下にプロットします。それぞれ、青色及び赤色で示されています。
FMQA 最適化サイクルにより得られた入力値 $\hat{x}$ により、目的関数値の最小値が次々と更新される様子が示されています(出力例を、『3.3. FMQA サンプルコード実行例』に紹介しています)。
fig = fmqa_solver.plot_history()
一般的に、FixstarsClient
で採用されているヒューリスティクスというアルゴリズムの原理上、得られる解に完全な再現性はありませんが、本サンプルコードを実行した際に得られる、典型的な標準出力及び画像出力を以下に紹介します。※得られる値が多少異なる場合がございます。
『3.1. $\boldsymbol{x}$の2次式に対する最適化』に記載の FMQA コードを与えられた条件のまま実行すると、次のような標準出力が逐次出力されます。
Generating 0-th training data set.
Generating 10-th training data set.
Generating 20-th training data set.
Generating 30-th training data set.
Generating 40-th training data set.
Generating 50-th training data set.
Starting FMQA cycles...
FMQA Cycle #0 variable updated, pred_y=-59.15752919611154
FMQA Cycle #1
FMQA Cycle #2 variable updated, pred_y=-72.66802296872575
FMQA Cycle #3
FMQA Cycle #4
FMQA Cycle #5
FMQA Cycle #6
FMQA Cycle #7
FMQA Cycle #8 variable updated, pred_y=-76.81540215271143
FMQA Cycle #9
pred x: [0. 0. 0. 0. 0. 0. 1. 1. 1. 1. 1. 0. 1. 0. 1. 1. 1. 1. 1. 1. 1. 0. 1. 0.
1. 0. 0. 0. 0. 0. 1. 1. 1. 0. 1. 1. 1. 1. 0. 1. 1. 1. 0. 1. 1. 0. 1. 1.
2. 1. 1. 1. 0. 1. 0. 1. 0. 1. 0. 0. 0. 1. 1. 1. 1. 1. 1. 0. 0. 0. 0. 0.
3. 0. 1. 1. 0. 1. 0. 0. 1. 0. 1. 1. 0. 1. 0. 1. 0. 1. 0. 1. 1. 0. 1. 0.
4. 0. 0. 1.]
pred value: -76.81540215271143
『3.2.FMQA 最適化過程における目的関数値の推移』に記載の fmqa_solver.plot_history()
の出力画像は次のようになります。
今回は、比較的単純な既知関数を対象として、FMQA による最適化を実施しました。しかし、本来、FMQA が力を発揮するのは、強い非線形性を有する物理現象や社会現象など、評価にシミュレーションや実験計測が必要な未知関数に対する最適化です。Amplify では、そのようなより現実的なモデルケースにおける FMQA の実施例やサンプルコードも紹介しています。
今回、最適化の対象であった $f(x) = x^{\top}Qx$ は2次式であり、数式が既知であるため、FMQA を使わずに直接組み合わせ最適化による最適入力値の探索が可能です。以下のコードでは、本関数に対して直接、量子アニーリング・イジングマシンによる最適化を行います。
# BinaryPoly の変数ジェネレータを宣言
gen = BinarySymbolGenerator()
# サイズ D を有する決定変数の1次元配列を作成
q = gen.array(D)
# xQx を QUBO の目的関数として定式化
cost = sum_poly(D, lambda i: sum_poly(D, lambda j: Q[i, j] * q[i] * q[j]))
# 目的関数を Amplify に渡し、求解を行う。
result = solver.solve(cost)
if len(result.solutions) == 0:
raise RuntimeError("No solution was found.")
# 推定された最適解を抽出し表示する。
values = result.solutions[0].values
true_x = q.decode(values)
print("true x:", true_x)
print("true value:", true_func(true_x))
本サンプルコードでご紹介したブラックボックス最適化手法は、次の研究で Factorization Mechine with Quantum Anealing (FMQA) として提案されたものです。
この論文では、FMQA により『メタマテリアル』の探索が実施されており、論文内では従来のブラックボックス最適化手法であるベイズ最適化に比べて優位な性能を示すことも報告されています。また、次の研究では、フォトニック結晶の設計において、本ブラックボックス最適化手法が活用されています。
これらの研究報告は、様々な分野におけるブラックボックス最適化問題において、FM と組合せ最適化に基づく本最適化手法 (FMQA) が汎用的に適用できる可能性を示唆しています。
また、機械学習とアニーリングに基づくブラックボックス最適化に関して、以下の解説記事をご紹介します。この解説記事において、本ブラックボックス最適化手法は、量子アニーリング (Quantum Annealing) から一般的なアニーリング (Annealing) に拡張した名称である Factorization Machine with Annealing (FMA) として紹介されています。
Fixstars Amplify では、材料探索だけでなく、化学反応、流体力学分野に関するブラックボックス最適化例も以下のとおり公開しています。本サンプルコードで紹介したブラックボックス最適化の応用例としてご覧ください。