数分割問題#

こちらでは、Lucas, 2014, “Ising formulations of many NP problems”の 2.1. Number Partitioning を OpenJij と JijModeling、そしてJijModeling transpiler を用いて解く方法をご紹介します。

概要: 数分割問題とは#

数分割問題は、与えられた数字の集合を足した合計値が等しくなるように2つの集合に分割する問題です。 ここで、簡単な例を考えてみましょう。

例えば、A={1,2,3,4}A=\{1,2,3,4\}という数字の集合AAがあるとします。 この集合を合計値が等しくなるように分割するのは簡単で、{1,4},{2,3}\{1,4\},\{2,3\}とすれば、それぞれの集合の合計値が5になるということがわかります。 このように、集合のサイズが小さい場合には、比較的簡単に答えがもとまりますが、これが大きくなるとすぐには解けません。 そこで、このチュートリアルでは、この問題をアニーリングを使って解いてみましょう。

from jijmodeling.transpiler.pyqubo import to_pyqubo
import openjij as oj
import jijmodeling as jm
import numpy as np

まず初めに、この問題のハミルトニアンを考えます。

分割する集合をAAとし、その要素をai(i={0,1,,N1})a_i (i = \{0,1,\dots,N-1\})とします。 ここでNNはこの集合の要素数です。 この集合AAを二つの集合をA0A_0A1A_1に分割するとします。 この時、xix_iAAii番目の要素が、集合A0A_0に含まれる時0、A1A_1に含まれる時1となる変数とします。 この変数xix_iを用いると、A0A_0に入っている数の合計値はiai(1xi)\sum_i a_i (1-x_i)とかけ、A1A_1iaixi\sum_i a_i x_iとなることがわかります。 この問題は、A0A_0A1A_1に含まれている数の合計値が等しくなるという制約を満たす解を求める問題ですので、これを式にすると、

iai(1xi)=iaixi\sum_i a_i (1-x_i)=\sum_i a_i x_i

という制約条件を満たすxix_iを求めよという問題であることがわかります。 これを式変形すると、iai(2xi)=0\sum_i a_i (2-x_i)=0と書くことができ、さらに、Penalty法を用いて、この制約条件を2乗したものをハミルトニアンとすると、結局、数分割問題のハミルトニアンは、

H=(i=0N1ai(2xi))2H=\left( \sum_{i=0}^{N-1} a_i (2-x_i)\right)^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
Problem: number partitionmin(i=0ashape(0)1ai(2xi1))2xi0{0,1}\begin{alignat*}{4}\text{Problem} & \text{: number partition} \\\min & \quad \left( \sum_{ i = 0 }^{ a_{\mathrm{shape}(0)} - 1 } a_{i} \cdot \left( 2 \cdot x_{i} - 1 \right) \right) ^ { 2 } \\& x_{i_{0}} \in \{0, 1\}\end{alignat*}

インスタンスデータの作成#

数理モデルができたのでここでは、1から40までの数字を分割するという問題を考えてみます。 NiN_{i}からNfN_{f}まで連続する数を分割する問題(連続する数の合計数が偶数の時)では、分割の仕方はいろんなパターンがありますが分割された集合の合計値は、

total value=(Ni+Nf)(NfNi+1)4\mathrm{total\ value} = \frac{(N_{i} + N_{f})(N_{f} - N_{i} + 1)}{4}

と計算することができます。 今考えている場合には、合計値は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)

ここでは、AAの中でA1A_1に分類されたindexとA0A_0に分類された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 : [ 6  8 10 11 13 15 17 18 19 20 24 28 31 36 37 38 39 40] , total value = 410
class 0 : [ 1  2  3  4  5  7  9 12 14 16 21 22 23 25 26 27 29 30 32 33 34 35] , total value = 410

我々の予想通り、それぞれの合計値410が得られていることがわかりました。