D-WaveのLeapを使ってナップサック問題を量子アニーリングで解いてみた

概要

前回記事(以下)では、D-Wave社が提供する量子アニーリングの開発環境であるLeapのAPIをPythonを使って動作させるところまで確認しました。本記事ではさらに、NP困難な組合せ最適化問題の1つであるナップサック問題を解かせてみて、どのような結果が出るかを検証してみました。

QUBO形式への定式化

定式化の方針

現状、最適化問題をLeapの量子アニーリングで解くためには、問題を以下のように評価関数Eの最小化問題、かつ最適化変数xに対する2次系式として定式化する必要があります。(下記サイト内ではQUBO形式という呼ばれ方をしています。)
https://docs.dwavesys.com/docs/latest/c_gs_3.html

$$Minimize(E) \\
E=x^{T}Cx \tag{eq1}$$

ここで、xは最適化変数(n行1列で要素は0か1)、Cはコストや制約等を諸々加味した行列(n行n列の3角行列で要素は実数)です。この形になるように問題を定式化することが必要になります。

$$x^{T}=[x_1, x_2 ,,, x_n] \\
x_{i} \in [0, 1] \tag{eq2}$$

$$C = \left[ \begin{array}{cccc} c_{11} & c_{12} & \ldots & c_{1n} \\ c_{21} & c_{22} & \ldots & c_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ c_{n1} & c_{n2} & \ldots & c_{nn} \end{array} \right]\\
c_{ij} \in R \tag{eq3}$$

ナップサック問題の場合、おおよその理解としては、どの荷物を詰めるかをxに対応させ、また詰め込んだ荷物の価値と荷物の重量制限の満足を考慮した式となるようにCを設計することになります。

ナップサック問題の定式化

各品目に対する入れる/入れないをx(n行1列)で表し、各品目の価値をv(n行1列)、各品目の重量をw(n行1列)、ナップサックの重量制限をl(スカラー)と置くと、量子アニーリングで解くことを考慮しなかった場合の、一般的なナップサック問題の定式化は以下になります。

$$Maximize(P) \\
P = x^{T} v \\
\text{subject to } x^{T} w < l \tag {eq4}$$

(eq1)と比較すると、最大化問題になっていることや、制約条件があること、2次系式となっていないことに違いがあります。以下順を追って変えていきます。

最大化問題⇒最小化問題

これは簡単で、評価関数にマイナスをかければOKです。

$$Minimize(P’) \\
P’ = -x^{T} v \\
\text{subject to } x^{T} w < l \tag {eq5}$$

不等式制約⇒ペナルティ関数化

不等式制約をそのまま組み込むことは出来ないので、制約を守っていれば0、制約を逸脱するとその逸脱幅に応じた大きな値になる、というペナルティ関数へ変換して組み込みます。

$$Minimize(P’) \\
P’ = -x^{T} v + \alpha (x^{T} w – l)^2\tag {eq6}$$

α以降がペナルティ部分になります。αは0以上の値とし、大きいほどペネルティの影響度が高くなります。ペナルティ部分では、重みがナップサックの制限であるl以上になると正の値となります。ただし、このままだと制限以下でも正の値となってしまうため、スラック変数と呼ばれる調整用の変数λを追加します。
https://docs.dwavesys.com/docs/latest/c_handbook_9.html?highlight=slack#csp-conversion-by-penalty-functions

$$Minimize(P’) \\
P’ = -x^{T} v + \alpha (x^{T} w – \lambda_1 – \lambda_2 – … – \lambda_m)^2\\
\lambda_{i} \in [0, 1] \tag {eq7}$$

コーディングを楽にするため、上記URLとは少し異なる形状に変形していますが、やりたいことは同じです。λは0か1を取る変数で、全てのλが1のとき、l-λの和が0になる、つまりλを最適化変数として加えることで、重量がl以下であればペナルティが0となるようにλを調整するように最適化することを意味します。

2次系式へ式変形

最後にeq7を展開して整理し、eq1に相当する形状へ式変形します。

$$Minimize(P’) \\
P’ = x’^{T} C x’ \\
x’ = [x_1, x_2 ,,, x_n, \lambda_1, \lambda_2 ,,, \lambda_m] \\
C = \left[\begin{array}{ccccccc} -v_{1}+\alpha w_{1}^{2} & 2\alpha w_{1}w_{2} & \ldots & 2\alpha w_{1}w_{n} & 2\alpha w_{1}c_{1} & \ldots & 2\alpha w_{1}c_{m}\\
0 & -v_{2}+\alpha w_{2}^{2} & \ldots & 2\alpha w_{2}w_{n} & 2\alpha w_{2}c_{1} & \ldots & 2\alpha w_{2}c_{m} \\
\vdots & 0 & \ddots & \vdots & \vdots & \ldots & \vdots \\
\vdots & \vdots & \ldots & -v_{2}+\alpha w_{n}^{2} & \vdots & \ldots & \vdots \\
\vdots & \vdots & \ldots & 0 & \alpha c_{1}^{2} & \ldots & 2\alpha c_{1}c_{m}\\
\vdots & \vdots & \ldots & \vdots & 0 & \ddots & \vdots \\
0 & 0 & \ldots & 0 & 0 & \ldots & 2\alpha c_{m}^{2}\\
\end{array} \right] \tag{eq8}$$

ここまでで定式化が出来たので、引き続き実装に移ります。

実装と検証

引き続き、Pythonを用いて実装を行います。ソースコードはGithubでも公開しています。
https://github.com/Rosyuku/Dwave_leap_trial/blob/master/knapsack%20problem/quantum.py

先ずパラメータを定義します。品目の種類は3種類とし、それぞれの価値(value)と重さ(weight)は乱数で与えます。また、重量制限はn*2としておきます。alphaもとりあえず1としておきます。

n = 3
value = np.random.randint(1, n, n)
weight = np.random.randint(1, n, n)
weight_limit = n * 2
alpha = 1
c = np.array([1]*weight_limit)

なお、valueとweightの値は以下になっています。

value:array([1, 2, 2])
weight:array([1, 1, 2])

上記のパラメータを用いて、eq8と等価になるようにQUBOモデルを構築します。

Q = dict()
for i in range(n+weight_limit):
    for j in range(n+weight_limit):
        if i==j:
            if i &amp;amp;amp;lt; n:
                Q.update({("x"+str(i), "x"+str(i)):-value[i]+alpha*weight[i]**2})
            else:
                Q.update({("x"+str(i), "x"+str(i)):alpha*c[i-n]**2})
        elif i &amp;amp;amp;gt; j:
            if (i &amp;amp;amp;lt; n) &amp;amp;amp;amp; (j &amp;amp;amp;lt; n):
                Q.update({("x"+str(i), "x"+str(j)):2*alpha*weight[i]*weight[j]})
            elif (i &amp;amp;amp;gt;= n) &amp;amp;amp;amp; (j &amp;amp;amp;lt; n):
                Q.update({("x"+str(i), "x"+str(j)):-2*alpha*c[i-n]*weight[j]})
            else:
                Q.update({("x"+str(i), "x"+str(j)):2*alpha*c[i-n]*c[j-n]})
        else:
            continue

最後にAPIを呼び出して結果を受け取ります。

response = EmbeddingComposite(DWaveSampler(token=setting.tokencode)).sample_qubo(Q, num_reads=1000)
print(list(response.data())[0])
print("Total_real_time ", response.info["timing"]["total_real_time"], "us")

n=3の結果

n=3を指定して計算した結果、以下の解が得られました。

Sample(sample={'x0': 1, 'x1': 1, 'x2': 1, 'x3': 0, 'x4': 1, 'x5': 1, 'x6': 1, 'x7': 1, 'x8': 0}, energy=-5.0, num_occurrences=1, chain_break_fraction=0.6666666666666666)

1つ目の品目を入れる[x0=1]、2つ目の品目を入れる[x1=1]、3つ目の品目を入れる[x2=1]が良く、スラック変数[x3以降]は4つが1になっており、重さの合計は4である。またその時の価値(energy)は-5.0になっており、1000回中1回(num_occurrences)がこの解にたどり着いた、になります。これは、厳密解と一致しており、正しく解けた結果と言えます。一方で、1件しか本解にたどり着いていませんが、スラック変数の組み合わせは冗長なので、これらを集約した結果も求めました。

x0 x1 x2 Energy             
 1  1  1   5.0             26
 0  1  1   4.0              6
 1  1  1   4.0            181
 0  1  1   3.0              7
 1  1  0   3.0              1
       1   1.0            487
 0  1  1  -0.0              1
 1  1  1  -4.0            291

x0、x1、x2を全て1にする解が全体の(985件/1000件)となっています。一方で、スラック変数の最適化はうまく出来ておらず、Energyが5.0になっているものは26件に限られました。このことから、量子アニーリングに投げ込みさえすれば最適解が得られる、というほど簡単なものではなく、適切なチューニングを実施しないと、周辺の局所解に陥るケースが多そうである、ということが想像できます。

n=4の結果

続いて、n=4の場合を検討します。value、weightは以下のようになっています。

value:array([1, 2, 1, 3])
weight:array([3, 1, 3, 2])

得られた結果は以下になります。

Sample(sample={'x0': 1, 'x1': 1, 'x2': 0, 'x3': 1, 'x4': 1, 'x5': 0, 'x6': 1, 'x7': 0, 'x8': 1, 'x9': 1, 'x10': 1, 'x11': 1}, energy=-6.0, num_occurrences=2, chain_break_fraction=0.8333333333333334)
Total_real_time  171733 us

n=4でも厳密解に一致する結果が得られました。

n=5の結果

続いて、n=5の場合を検討します。value、weightは以下のようになっています。

value:array([1, 4, 2, 1, 3])
weight:array([4, 3, 4, 1, 4])

得られた結果は以下になります。

Sample(sample={'x0': 1, 'x1': 0, 'x2': 1, 'x3': 1, 'x4': 1, 'x5': 1, 'x6': 1, 'x7': 1, 'x8': 0, 'x9': 1, 'x10': 1, 'x11': 1, 'x12': 1, 'x13': 0, 'x14': 1}, energy=18.0, num_occurrences=1, chain_break_fraction=0.8666666666666667)
Total_real_time  171762 us

energyの値が18.0と正の値になっており、制約条件を満足できていない(つまり、解けていない)ことが分かります。実際この問題の厳密解は[x0=0, x1=1, x2=0, x3=1, x4=1]の時に、value合計=8、weight合計=8であり、正しい答えが求まっていません。なんと、純粋にナップサック問題を解くと、(ライブラリの理解等が進めば変わるのかもしれませんが、単純にはn=5の時点で答えが求まらなくなってしまいました。

まとめ

ナップサック問題を定式化し、量子アニーリングを用いて最適解を算出しました。すると、n=3及び4では解が求まりましたが、n=5というまだ簡単な条件の時点で解が求まらなくなるというやや驚きの結果となりました。また解が求まったケースでも不等式制約の部分は1000件のうちの数十件と、数%の割合でしか解けておらず、デフォルト設定のまま解かせるとかなり局所解に陥りやすい傾向を確認しました。次回は、最適化変数を減らすために、不等式制約を厳しくした条件で求解させてみて、どのような傾向になるかを確認してみようと思います。

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です

CAPTCHA