FMQA によるブラックボックス最適化
FMQA（Factorization Machine with Quantum Annealing）は、評価に時間や費用のかかる関数を、少ない評価回数で最小化する方法です。 評価した結果から関数の近似モデルを学習し、そのモデルを最小にする点を QUBO ソルバーで求めて次に評価する、という手順を繰り返します。 名前の QA は量子アニーリングのことですが、QUBO を解く部分にはどの QUBO ソルバーも使えます。このページでは EasySolver を使います。
前提: 評価が高価なブラックボックス関数
$n$ 個のバイナリ変数 $x = (x_1, \ldots, x_n)$ の関数 $f(x)$ を最小化したいとします。 ただし $f$ の式は分からず、$x$ を与えると値が返ってくるだけです（ブラックボックス関数）。 しかも 1 回の評価に実験や時間のかかるシミュレーションが必要で、評価できる回数は限られています。 たとえば、20 種類の材料それぞれを使うか使わないかを決めて試作品を作り、その性能を測る、といった状況です。
$n = 20$ でも $x$ は $2^{20}$（約 100 万）通りあり、すべてを試すことはできません。 このページでは、評価できる回数を $300$ 回として、その中でできるだけ小さい $f(x)$ を見つけます。
FMQA の手順
- ランダムに選んだ点をいくつか評価する。
- それまでに評価したすべての点 $(x, f(x))$ に、Factorization Machine（FM）という 2 次のモデルを当てはめる（学習）。
- 学習した FM の予測式は $x$ の 2 次式、つまり QUBO なので、QUBO ソルバーで予測値を最小にする $x$ を求める。
- その $x$ を評価してデータに加え、2. に戻る。
モデルが小さいと予測した点を実際に評価し、予測が外れていれば、その結果が次の学習に入ってモデルが直ります。 この「予測、確認、修正」を繰り返して、少ない評価回数で最小値に近づきます。 QUBO ソルバーは $2^n$ 通り全体の中からモデルの最小点を選ぶので、何ビットも同時に変えた点へ一度に移れます。
この例のブラックボックス関数
このページでは、高価な評価の代わりに次の式を使います（$n = 20$）:
\[f(x) = 100 \sum_{k=1}^{3} \sin\Bigl(\sum_{i=1}^{20} a_{ki}\, x_i\Bigr)\]係数 $a_{ki}$ は $-1$ から $1$ までの数で、プログラムの中に表として書いてあります。 この式は評価の代わりと答え合わせにだけ使い、FMQA の部分は式の中身を使いません。
- $f(x)$ は $-300$ 以上 $300$ 以下で、3 つの和 $\sum_i a_{ki} x_i$ がそろって $\sin$ の谷（$-\pi/2 + 2\pi m$）の近くにあるとき、$-300$ に近くなります。 全 $2^{20}$ 通りを調べると、最小値は $-299.92$ です。
- $\sin$ の和なので、$x$ の 2 次式（QUBO）ではありません。FMQA はこれを 2 次式で近似しながら探します。
- 良い解はまれです。最小値の 99% 以内（$-296.92$ 以下）に入る $x$ は $51$ 個（全体の約 $0.005\%$）しかなく、 ランダムに 300 個選んでも、その中に入っている確率は $1.4\%$ です。
Factorization Machine
FM は、$f(x)$ を次の 2 次式で予測するモデルです:
\[\hat{y}(x) = w_0 + \sum_{i} w_i x_i + \sum_{i<j} \langle v_i, v_j\rangle\, x_i x_j\]$w_0$ と $w_i$ は実数、$v_i$ は変数 $x_i$ ごとの $D$ 次元の実数ベクトルで、$\langle v_i, v_j\rangle$ はその内積です。
- 2 次の係数を内積で表す: 2 次の係数を 1 つずつ推定すると、$n = 20$ では $190$ 個あります。 FM ではベクトル $v_i$ の成分だけを推定するので、$D = 3$ なら $60$ 個です。 $w_0$ と $w_i$ を合わせても $81$ 個なので、評価した点が少ないうちから学習できます。
- $D = 3$ にした理由: この例の $f$ は、3 つの 1 次式 $\sum_i a_{ki} x_i$ を通してだけ $x$ に依存します。 それぞれの $\sin$ を谷の近くで放物線 $\gamma_k \bigl(\sum_i a_{ki} x_i - c_k\bigr)^2$（$\gamma_k > 0$）で近似すると、 $x_i x_j$ の係数は $2 \sum_k \gamma_k a_{ki} a_{kj}$ となり、3 次元のベクトルの内積でちょうど表せます。 実際の問題ではこのような構造は分からないので、$D$ はいくつか試して決めます。 また、$x_i^2 = x_i$ を使うと、2 次の項は次のように書き換えられます（$v_{id}$ は $v_i$ の第 $d$ 成分）:
学習でも QUBO の式を作るときも、この形を使います。
学習
評価した点の $f$ の値を平均 $0$、標準偏差 $1$ にそろえたものを $y$ とします。 予測値と $y$ の二乗誤差の平均に、正則化項 $\lambda \bigl(\sum_i w_i^2 + \sum_{i,d} v_{id}^2\bigr)$（$\lambda = 10^{-4}$）を加えたものを、 Adam という勾配法で最小化します（$1000$ 回の更新）。 勾配は予測式を微分したもので、たとえば $\partial \hat{y} / \partial v_{id} = x_i \bigl(\sum_j v_{jd}\, x_j - v_{id}\bigr)$ です。 FM は評価のたびに、$v_{id}$ を乱数で初期化し直して最初から学習します。
QUBO++ プログラム
以下のプログラムは、評価回数 300 回で FMQA を実行し、最後に答え合わせとして全 $2^{20}$ 通りの中の最小値を表示します:
import numpy as np
import pyqbpp.d as qbpp
N = 20 # ビット数
K = 3 # 黒箱の中の sin の数
D = 3 # FM のベクトルの次元
INIT = 20 # 最初にランダムに選んで評価する点の数
EVALS = 300 # 黒箱を評価できる回数
# 黒箱の係数 a_ki
A = np.array([
[0.02, 0.90, -0.71, 0.90, -0.38, -0.15, 0.66, -0.18, 0.10, -0.94,
0.51, 0.08, -0.34, 0.58, -0.39, -0.09, -0.73, -0.19, -0.59, -0.48],
[0.50, -0.44, -0.03, 0.96, 0.92, 0.45, 0.08, -0.45, -0.68, 0.94,
0.03, -0.77, 0.25, 0.55, 0.23, 0.83, -0.92, 0.06, -0.08, -0.88],
[0.28, 0.71, 0.19, -0.48, 0.68, 0.02, 0.02, 0.51, -0.70, 0.64,
0.37, 0.57, -0.62, 0.60, -0.62, -0.84, 0.71, 0.72, 0.75, -0.06]])
# 黒箱 f(x) = 100 * Σ_k sin(Σ_i a_ki x_i)
# 本当は高価な測定やシミュレーションで、ここでは式で代用する
def blackbox(b):
return 100 * np.sin(A @ b).sum()
# 予測値と y の二乗誤差の平均 + 正則化項 を Adam で最小化して FM を学習する
def train(X, y, rng):
p = [np.zeros(1), np.zeros(N), rng.normal(0, 0.1, (N, D))] # w0, w, v
m = [np.zeros_like(a) for a in p]
u = [np.zeros_like(a) for a in p]
lr, b1, b2, l2 = 0.05, 0.9, 0.999, 1e-4
for t in range(1, 1001):
w0, w, v = p
# 予測値 w0 + Σ_i w_i x_i + (1/2) Σ_d [(Σ_i v_id x_i)^2 - Σ_i v_id^2 x_i]
s = X @ v # s[n, d] = Σ_i v_id x_i
pred = w0[0] + X @ w + 0.5 * ((s**2).sum(1) - X @ (v**2).sum(1))
# 勾配
r = 2 * (pred - y) / len(y)
g = [np.array([r.sum()]),
X.T @ r + 2 * l2 * w,
X.T @ (r[:, None] * s) - v * (X.T @ r)[:, None] + 2 * l2 * v]
# Adam による更新
for j in range(3):
m[j] = b1 * m[j] + (1 - b1) * g[j]
u[j] = b2 * u[j] + (1 - b2) * g[j]**2
mh = m[j] / (1 - b1**t)
uh = u[j] / (1 - b2**t)
p[j] = p[j] - lr * mh / (np.sqrt(uh) + 1e-8)
return p
rng = np.random.default_rng(7)
X = [] # 評価した点
y = [] # その点の f の値
seen = set() # 評価済みの点
def evaluate(b):
X.append(b)
y.append(blackbox(b))
seen.add(tuple(b))
# 1. ランダムに選んだ点を評価する
while len(X) < INIT:
b = rng.integers(0, 2, N)
if tuple(b) not in seen:
evaluate(b)
best = min(y)
print(f"eval {len(X)}: f = {best:.2f}")
x = qbpp.var("x", N)
while len(X) < EVALS:
# 2. f の値を平均 0、標準偏差 1 にそろえて FM を学習する
z = (np.array(y) - np.mean(y)) / np.std(y)
w0, w, v = train(np.array(X), z, rng)
# 3. FM の予測式を QUBO++ の式にして、最小にする x を求める
f = w0[0]
for i in range(N):
f += w[i] * x[i]
for d in range(D):
s = 0
q = 0
for i in range(N):
s += v[i, d] * x[i]
q += v[i, d] ** 2 * x[i]
f += 0.5 * (qbpp.sqr(s) - q)
f.simplify_as_binary()
sol = qbpp.EasySolver(f).search(time_limit=0.1)
b = np.array(sol(x), dtype=int)
# 4. 評価済みの点なら、ランダムなビットを反転してから評価する
while tuple(b) in seen:
i = rng.integers(N)
b[i] = 1 - b[i]
evaluate(b)
if y[-1] < best:
best = y[-1]
print(f"eval {len(X)}: f = {best:.2f}")
# 答え合わせ: 全 2^N 通りの中の最小値
B = (np.arange(2**N)[:, None] >> np.arange(N)) & 1
fmin = (100 * np.sin(B @ A.T).sum(axis=1)).min()
print(f"minimum over all 2^{N} points: {fmin:.2f}")
import pyqbpp.dで、式の係数を実数（float）にしています（実数（double）係数）。 学習した FM のパラメータは実数なので、そのまま QUBO++ の式の係数にできます。学習には numpy を使います（pip install numpy）。blackboxが評価の代わりの式です。FMQA の部分はblackboxが返す値だけを使います。- FM のパラメータは
[w0, w, v]（w0は長さ 1、wは長さN、vはN×Dの numpy 配列）として持ちます。trainは、評価したすべての点の予測値と勾配を行列の計算でまとめて求め、Adam で更新します。 - 手順 3 では、学習した予測式をそのまま QUBO++ の式
fに書き、simplify_as_binary()で $x_i^2 = x_i$ として整理してからEasySolverで解きます。 変数は 20 個なので、求解は 0.1 秒で十分です。 - 学習が進むと、すでに評価した点が提案されることがよくあります。 同じ点を評価し直しても情報は増えないので、手順 4 で、まだ評価していない点になるまでランダムなビットを反転します。
- 実行には 1〜2 分ほどかかります。
出力結果
eval 20: f = -200.02
eval 28: f = -236.44
eval 32: f = -255.80
eval 37: f = -261.30
eval 41: f = -292.94
eval 89: f = -293.16
eval 113: f = -297.93
minimum over all 2^20 points: -299.92
最初にランダムに選んだ 20 点の最良値は $-200.02$ で、113 回目の評価で $-297.93$ に達しています（全 $2^{20}$ 通りの最小値は $-299.92$）。 学習の初期値などに乱数を使っているので、乱数の種を変えると結果も変わります。
ランダム探索・山登り法との比較
乱数の種を変えて上のプログラムを 20 回実行し、評価回数ごとの最良値（中央値）を、次の 2 つの方法と比べました:
- ランダム探索: 各ビットを確率 $1/2$ で決めた $x$ を評価することを繰り返します（値は全 $2^{20}$ 通りの分布から計算）。
- 山登り法: ランダムな $x$ から始め、ランダムな順に 1 ビットずつ反転して評価し、$f$ が小さくなればそこへ移ります。 どのビットを反転しても小さくならない点（局所最適解）に着いたら、新しいランダムな $x$ からやり直します。 評価済みの点は評価し直さず、500 回実行した中央値です。
| 評価回数 | FMQA | ランダム探索 | 山登り法 |
|---|---|---|---|
| 50 | $-285.32$ | $-209.88$ | $-290.96$ |
| 100 | $-293.61$ | $-235.38$ | $-294.11$ |
| 200 | $-297.93$ | $-255.48$ | $-296.06$ |
| 300 | $-298.55$ | $-264.01$ | $-296.67$ |
| 300 回以内に最小値の 99% 以内に入った割合 | $90\%$ | $1.4\%$ | $49\%$ |
- FMQA はランダム探索よりずっと小さい値に届きます。
- 山登り法は、評価回数が少ないうちは FMQA より良い値を出しますが、$-296$ 付近で止まりがちです。 この $f$ には局所最適解が $357$ 個あり、その多くは $-290$ 前後の値です。 そこから最小値の近くへ行くには、3 つの和を同時に谷底へ合わせる必要があり、何ビットも同時に変えなければなりません。
- FMQA は谷の形を FM で学習し、その最小点として何ビットも離れた点へ一度に移れるので、評価 100〜200 回で山登り法を追い越します。
この例の QUBO は変数が 20 個と小さく、ソルバーにとっては簡単な問題です。 結果を左右しているのは主に学習の部分です。 変数が数百以上に増えたり、$x$ に制約があったりすると、QUBO を正しく解く部分の役割が大きくなります。