数分割問題#
こちらでは、Lucas, 2014, “Ising formulations of many NP problems”の 2.1. Number Partitioning を OpenJij と JijModeling、そしてJijModeling transpiler を用いて解く方法をご紹介します。
概要: 数分割問題とは#
数分割問題は、与えられた数字の集合を足した合計値が等しくなるように2つの集合に分割する問題です。 ここで、簡単な例を考えてみましょう。
例えば、という数字の集合があるとします。 この集合を合計値が等しくなるように分割するのは簡単で、とすれば、それぞれの集合の合計値が5になるということがわかります。 このように、集合のサイズが小さい場合には、比較的簡単に答えがもとまりますが、これが大きくなるとすぐには解けません。 そこで、このチュートリアルでは、この問題をアニーリングを使って解いてみましょう。
from jijmodeling.transpiler.pyqubo import to_pyqubo
import openjij as oj
import jijmodeling as jm
import numpy as np
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[1], line 1
----> 1 from jijmodeling.transpiler.pyqubo import to_pyqubo
2 import openjij as oj
3 import jijmodeling as jm
ModuleNotFoundError: No module named 'jijmodeling.transpiler'
まず初めに、この問題のハミルトニアンを考えます。
分割する集合をとし、その要素をとします。 ここではこの集合の要素数です。 この集合を二つの集合をとに分割するとします。 この時、をの番目の要素が、集合に含まれる時0、に含まれる時1となる変数とします。 この変数を用いると、に入っている数の合計値はとかけ、のとなることがわかります。 この問題は、とに含まれている数の合計値が等しくなるという制約を満たす解を求める問題ですので、これを式にすると、
という制約条件を満たすを求めよという問題であることがわかります。 これを式変形すると、と書くことができ、さらに、Penalty法を用いて、この制約条件を2乗したものをハミルトニアンとすると、結局、数分割問題のハミルトニアンは、
となります。
JijModelingによるモデル構築#
上記のハミルトニアンをJijModelingを用いて定式化していきます。
problem = jm.Problem("number partition")
a = jm.Placeholder("a",dim = 1)
N = a.shape[0]
x = jm.Binary("x",shape=(N,))
i = jm.Element("i",(0,N))
s_i = 2*x[i] - 1
problem += (jm.Sum(i,a[i] * s_i))**2
problem
インスタンスデータの作成#
数理モデルができたのでここでは、1から40までの数字を分割するという問題を考えてみます。 からまで連続する数を分割する問題(連続する数の合計数が偶数の時)では、分割の仕方はいろんなパターンがありますが分割された集合の合計値は、
と計算することができます。 今考えている場合には、合計値は410になることが予想されます。 実際に確かめてみましょう。
N = 40
instance_data = {"a":np.arange(1,N+1)}
JijModeling transpilerによるPyQUBOへの変換#
問題の入力データができたので、これをQUBOに変換します。
model,cache = to_pyqubo(problem,instance_data,{})
Q,offset = model.compile().to_qubo()
OpenJijによる最適化計算の実行#
OpenJijを用いて計算してみます。
sampler = oj.SASampler(num_reads=1)
res = sampler.sample_qubo(Q=Q)
デコードと解の表示#
得られた結果を見てみましょう。返された計算結果をデコードし、解析を行いやすくします。
decoded = cache.decode(res)
ここでは、の中でに分類されたindexとに分類されたindexを分けて、それらについて和をとっています。
class_1_index = decoded.record.solution['x'][0][0][0]
class_0_index = [i for i in range(0,N) if i not in class_1_index]
class_1 = instance_data['a'][class_1_index]
class_0 = instance_data['a'][class_0_index]
print(f"class 1 : {class_1} , total value = {np.sum(class_1)}")
print(f"class 0 : {class_0} , total value = {np.sum(class_0)}")
class 1 : [ 1 4 5 8 10 11 12 14 16 21 22 23 26 30 31 32 33 35 37 39] , total value = 410
class 0 : [ 2 3 6 7 9 13 15 17 18 19 20 24 25 27 28 29 34 36 38 40] , total value = 410
我々の予想通り、それぞれの合計値410が得られていることがわかりました。