数分割問題#
こちらでは、Lucas, 2014, “Ising formulations of many NP problems”の 2.1. Number Partitioning を OpenJij と JijModeling、そしてommx-openjij-adapter を用いて解く方法をご紹介します。
概要: 数分割問題とは#
数分割問題は、与えられた数字の集合を足した合計値が等しくなるように2つの集合に分割する問題です。 ここで、簡単な例を考えてみましょう。
例えば、\(A=\{1,2,3,4\}\)という数字の集合\(A\)があるとします。
この集合を合計値が等しくなるように分割するのは簡単で、\(\{1,4\},\{2,3\}\)とすれば、それぞれの集合の合計値が5になるということがわかります。
このように、集合のサイズが小さい場合には、比較的簡単に答えがもとまりますが、これが大きくなるとすぐには解けません。
そこで、このチュートリアルでは、この問題をアニーリングを使って解いてみましょう。
まず初めに、この問題の数理モデルを考えます。
分割する集合を\(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\)に含まれている数の合計値が等しくなるという制約を満たす解を求める問題ですので、これを式にすると、
という制約条件を満たす\(x_i\)を求めよという問題になります。
JijModelingによるモデル構築#
式(1)をJijModelingを用いて定式化していきます。
import jijmodeling as jm
problem = jm.Problem("Number Partition")
a = jm.Placeholder("a", ndim=1)
N = a.len_at(0, latex="N")
x = jm.BinaryVar("x", shape=(N, ))
i = jm.Element("i", (0, N))
problem += jm.Constraint("equal", jm.sum(i, a[i]*(1-x[i])) == jm.sum(i, a[i]*x[i]))
Jupyter Notebookで実装の確認を行いましょう。
problem
インスタンスデータの作成#
ここでは、1から40までの数字を分割する問題を考えましょう。 \(N_{i}\)から\(N_{f}\)まで連続する数を分割する問題(連続する数の合計数が偶数の時)では、分割の仕方はいろんなパターンがありますが分割された集合の合計値は
と計算することができます。 今の場合には、合計値は410となります。 実際にそれを確かめてみましょう。
import numpy as np
inst_N = 40
instance_data = {"a": np.arange(1, inst_N+1)}
JijModeling Interpreterを用いてOMMX Instanceを作成#
instance = jm.Interpreter(instance_data).eval_problem(problem)
OpenJijによる最適化計算の実行#
OpenJijを用いて計算してみます。
import ommx_openjij_adapter as oj_ad
import openjij as oj
adapter = oj_ad.OMMXOpenJijSAAdapter(instance)
best_sample = adapter.sample(instance, num_reads=10).best_feasible_unrelaxed()
デコードと解の表示#
計算結果をデコードします。 ここでは、\(A\)の中で\(A_1\)に分類されたindexと\(A_0\)に分類されたindexを分けて、それらについて和をとっています。
# decode a result to JijModeling sampleset
# get the indices of x == 1x
x_value = best_sample.extract_decision_variables("x")
class_1_indices = [k[0] for k, v in x_value.items() if v == 1]
class_0_indices = [k[0] for k, v in x_value.items() if v == 0]
class_1 = instance_data['a'][class_1_indices]
class_0 = instance_data['a'][class_0_indices]
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 : [ 3 4 6 7 8 10 11 12 13 14 16 17 18 19 21 28 30 31 32 33 37 40] , total value = 410
class 0 : [ 1 2 5 9 15 20 22 23 24 25 26 27 29 34 35 36 38 39] , total value = 410
合計値410が得られていることがわかりました。