数分割問題

数分割問題#

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

概要: 数分割問題とは#

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

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

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

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

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

\sum_i a_i (1-x_i)=\sum_i a_i x_i

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

H=\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

\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までの数字を分割するという問題を考えてみます。 $N_{i}$ から $N_{f}$ まで連続する数を分割する問題(連続する数の合計数が偶数の時)では、分割の仕方はいろんなパターンがありますが分割された集合の合計値は、

\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)

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

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