数分割問題#
こちらでは、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による定式化#
次に、JijModelingを使って上記の数理モデルを定式化してみましょう。まずは、数理モデルに現れる変数およびパラメーターを定義します。
import jijmodeling as jm
problem = jm.Problem("Number Partition")
a = problem.Float("a", ndim=1)
N = problem.DependentVar("N", a.len_at(0))
x = problem.BinaryVar("x", shape=(N, ))
aは、\(A\)にある要素を表す1次元配列です。
aの長さから、\(A\)に含まれる要素数Nを得ています。
最適化に用いるバイナリ変数をxで定義します。
制約#
制約式は以下のように定式化しましょう。
problem += problem.Constraint("equal", (jm.sum(N, lambda i: a[i]*(1-x[i])) == jm.sum(N, lambda 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)}
OpenJijによる最適化計算の実行#
OpenJijのシミュレーテッド・アニーリングを用いて、最適化問題を解いてみましょう。
from ommx_openjij_adapter import OMMXOpenJijSAAdapter
instance = problem.eval(instance_data)
adapter = OMMXOpenJijSAAdapter(instance)
best_sample = adapter.sample(instance, num_reads=10, uniform_penalty_weight=0.8).best_feasible_unrelaxed
解の可視化#
ここでは、\(A\)の中で\(A_0\)に分類されたindexと\(A_1\)に分類されたindexを分けて、それらについて和をとっています。
# decode a result to JijModeling sampleset
# get the indices of x == 1
df = best_sample.decision_variables_df
class_0_indices = [row['subscripts'][0] for _, row in df.iterrows() if row['name'] == 'x' and row['value'] <= 0.5]
class_1_indices = [row['subscripts'][0] for _, row in df.iterrows() if row['name'] == 'x' and row['value'] > 0.5]
class_0 = instance_data['a'][class_0_indices]
class_1 = instance_data['a'][class_1_indices]
print(f"class 0 : {class_0} , total value = {np.sum(class_0)}")
print(f"class 1 : {class_1} , total value = {np.sum(class_1)}")
class 0 : [ 1 2 5 9 18 21 22 25 26 29 30 32 36 37 38 39 40] , total value = 410
class 1 : [ 3 4 6 7 8 10 11 12 13 14 15 16 17 19 20 23 24 27 28 31 33 34 35] , total value = 410
合計値410が得られていることがわかりました。