量子ゲートの問題設定の前に、量子アニーリングのそれを出来るようになっておきたい…
ということで、逆評定問題（勝手に命名）をMDR社のQUBOシミュレータであるWildqatを用いて解いてみた。
https://github.com/Blueqat/Wildqat
問題設定
大学生活における「最小努力で必要単位数を取得する」ケース※$\tiny{1}$について考える。
或る履修可能科目群があるとして、その履修によって得られる単位を$u_i$、それらを取得するために必要な工数※$\tiny{2}$を$w_i$とし、科目$i$を履修するか否かを$x_i（x_i \in \{0, 1\}）$で表すとする。
最低でも$U$単位は取得しなければならないという前提条件にたつと、これは 取得単位の制約条件を満たしつつ総工数を最小化する最適化問題（0-1整数計画問題） と考えられる。命名、 逆評定問題 。
※$\tiny{1}$
ちなみに私は社会人なのでもうそんなケースはない。
※$\tiny{2}$
「各科目ごとの単位の取りやすさをまとめた何か」を想定。
自分の思いつく事例でいう逆評定。の評価を数値化したもの。
適宜、適切な用語に置き換えて頂けると。。
関係式の整理
定式化
最小努力で必要単位を取得したいので、今回の最小化対象は
\displaystyle \sum_{ i = 1 }^{ N } w_i x_i, x_i \in \{0, 1\}
またそれに際して、必要単位以上は取得しなければならないので制約条件は
\displaystyle \sum_{ i = 1 }^{ N } u_i x_i \geq U
これに調整パラメータ$\alpha$を加えて目的関数（＝最小化対象）とコスト関数（＝制約条件）を1つのハミルトニアンとしてまとめ、通例の1つ？に従って最小化対象を
\begin{align}
H &= H_{obj} + H_{cost}\\
&= \displaystyle \sum_{ i = 1 }^{ N } w_i x_i + \alpha (\sum_{ i = 1 }^{ N } u_i x_i - U)^2
\end{align}
と置いてみると、このままでは
- 取得単位に1満たなかった場合
- 取得単位を1超えた場合
を同様の重みで評価してしまうため、もう少し式変形を加える。
\begin{align}
H &= \displaystyle \sum_{ i = 1 }^{ N } w_i x_i + \alpha \{\sum_{ i = 1 }^{ N } u_i x_i - (U + \sum_{ j = 1 }^{ 3 } \lambda_j)\}^2\\
&= \displaystyle \mathbf{x}^T \mathbf{w} + \alpha\{\mathbf{x}^T \mathbf{u} - (5 + \mathbf{\lambda})\}^2, \lambda_j \in \{0, 1\}
\end{align}
補足
「取得単位の組み合わせの有益な候補として挙がるものは、高々最低取得単位数よりも3多い程度に収まるだろう」という想定のもと、チューニングパラメータ$\lambda$を3個に制限し、$\lambda$自体も最適化対象に加えた。
QUBO形式に変形
上記でまとめた最適化対象のハミルトニアンに対して、
\mathbf{P} = \mathbf{x^{'}}^T C \mathbf{x^{'}}
の形になるよう式変形を行う。ただし $\mathbf{x^{'}} = [x_1, x_2, ... ,x_N, \lambda_1, \lambda_2, \lambda_3] $
色々頑張ると、下記のようになった。
\mathbf{C} =
\left(\begin{array}{ccc}
w_{1} + \alpha u_{1}^2 - 10u_{1}& 2\alpha u_{1}u_{2} & 2\alpha u_{1}u_{3} & \cdots & 2\alpha u_{1}u_{N} & -2\alpha u_{1} & \cdots & & -2\alpha u_{1} \\
& \ddots & 2\alpha u_{2}u_{3} & \ddots & \vdots & \vdots & & & \vdots \\
& & \ddots & & 2\alpha u_{N-2}u_{N-1} \\
& & & \ddots & 2\alpha u_{N-1}u_{N}\\
& & & & w_{N}+\alpha u_{N}^2 - 10u_{N} & -2\alpha u_{N} & \cdots & & -2\alpha u_{N} \\
& & & & & 11\alpha & 2\alpha & \cdots & 2\alpha \\
& & & & & & 11\alpha & & \vdots \\
& & & & & & & \ddots & 2\alpha\\
& & & & & & & & 11\alpha \\
\end{array}\right)
あんまりきれいじゃない…
計算の実行
単位と工数の想定
Wildqatに読み込ませるところまで整理が出来たので、実際の数値を決めて関数を定義していく。
# u_i : 科目iの履修による取得単位数
u = [1, 3, 2, 2, 1]
# w_i : 科目iの単位取得難易度＝単位取得所要時間（工数）
w = [1, 5, 3, 1, 2]
# 行列サイズ
N = len(u)
# 最低取得単位数+3まで許容する
U = 3
# 最適化パラメータα
alpha = 1.5
つまり、下記のようなテーブルが与えられていると考える。
（「想定工数」箇所は、仏が1.0とか鬼が5.0、とかそういうイメージ）
この程度なら暗算でも確認可能であり、科目[A, C, D]の組み合わせが最小工数を取りそうである。
つまりシミュレーション実行結果としては $\mathbf{x}=[1,0,1,1,0]$ が出てくれるのが望ましい。
| 科目 | 取得単位[単位] | 想定工数[人月] |
|---|---|---|
| A | 1.0 | 1.0 |
| B | 3.0 | 5.0 |
| C | 2.0 | 3.0 |
| D | 2.0 | 1.0 |
| E | 1.0 | 2.0 |
関数の定義
状況が整理できたので、QUBO形式である$C$を計算するため関数を定義する。
import numpy as np
def get_qubo(u: list,w: list, N: int, U: int):
# 行列の合計サイズ
size = N + U
qubo = np.zeros((size, size))
# 1. x vector part
# 対角成分
for i in range(N):
for j in range(N):
if j == i:
qubo[i][j] = w[i] + alpha * (u[i] ** 2) - 10 * alpha * u[i]
# 上三角行列部分
for i in range(N):
for j in range(N):
if j >= i + 1:
qubo[i][j] = 2 * alpha * u[i] * u[j]
for i in range(N):
for j in range(N, N + U):
qubo[i][j] = -2 * alpha * u[i]
# 2. lambda vector part
# 対角成分
for i in range(N, N + U):
for j in range(N, N + U):
if j == i:
qubo[i][j] = 11 * alpha
# 上三角行列部分
for i in range(N, N + U):
for j in range(N, N + U):
if j >= i + 1:
qubo[i][j] = 2 * alpha
return qubo
計算
QUBOの定義も完了したので、あとはWlidqatに渡すだけ。
100回試行してみて、$\mathbf{x}（\mathbf{x^{'}}）$の分布を見てみる。
from blueqat import opt
a = opt.opt()
a.qubo = qubo
result = a.sa(shots = 100, sampler = "fast")
opt.counter(result)
今回の実行結果は下記の通りとなった。
# 実行結果 ： 実行のたびに解の分布は変わることに注意
# α=1.0の場合
Counter({'10110000': 25,'00110000': 31,'10011000': 24,'01010000': 19,'11000000': 1})
# α=1.5の場合
Counter({'01010000': 29,'11100010': 1,'10110000': 55,'10011000': 10,'11010001': 1,'11010010': 1,'01100000': 3})
# α=2.0の場合
Counter({'01010000': 23,'10110000': 47,'10011000': 6,'01110110': 2,'11100100': 2,'11010100': 6,'11010010': 6,'01100000': 3,'11010001': 2,'11000000': 2,'01110101': 1})
どうやら$\alpha$=1.5 辺りで$[1,0,1,1,0]$ の確率が最大となりそうである。
解釈
$\alpha$ = 1.0の場合に多く導出されている$[0,0,1,1,0],[1,0,0,1,1]$は単位数自体が足りていないため、コスト項がうまく影響しきれていない。また$\alpha$ = 1.5, 2.0の場合の次点である$[0,1,0,1,0]$は単位数は満たしているものの、最適解より工数が1多いので局所最適解である。
科目数や取得難易度のバラつきによってはペナルティが二乗で効いてくるコスト項の影響が大きくなりすぎるので、$\alpha$の値調整が難しい。。
解釈違い、誤りなどあればご指摘頂けると非常に有難いです。
それと、数式下のスクロールバーどうやって消すんだろう・・・