数分割問題#

こちらでは、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になるということがわかります。 このように、集合のサイズが小さい場合には、比較的簡単に答えがもとまりますが、これが大きくなるとすぐには解けません。 そこで、このチュートリアルでは、この問題をアニーリングを使って解いてみましょう。
まず初めに、この問題のハミルトニアンを考えます。 分割する集合を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(12xi)=0\sum_i a_i (1-2x_i)=0と書くことができ、さらに、Penalty法を用いて、この制約条件を2乗したものをハミルトニアンとすると、結局、数分割問題のハミルトニアンは、

H=(i=0N1ai(12xi))2(1) H=\left( \sum_{i=0}^{N-1} a_i (1-2x_i)\right)^2 \tag{1}

となります。

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.sum(i, a[i]*(1-2*x[i])) ** 2

Jupyter Notebookで実装の確認を行いましょう。

problem
Problem:Number Partitionmin((i=0N1ai(2xi+1))2)wherex1-dim binary variable\begin{array}{cccc}\text{Problem:} & \text{Number Partition} & & \\& & \min \quad \displaystyle \left(\left(\sum_{i = 0}^{N - 1} a_{i} \cdot \left(- 2 \cdot x_{i} + 1\right)\right)^{2}\right) & \\\text{{where}} & & & \\& x & 1\text{-dim binary variable}\\\end{array}

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

ここでは、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となります。 実際にそれを確かめてみましょう。

import numpy as np

inst_N = 40
instance_data = {"a": np.arange(1, inst_N+1)}

JijModeling transpilerによるPyQUBOへの変換#

ここまで行われてきた実装は、全てJijModelingによるものでした。 これをPyQUBOに変換することで、OpenJijはもちろん、他のソルバーを用いた組合せ最適化計算を行うことが可能になります。

import jijmodeling_transpiler as jmt

# compile
compiled_model = jmt.core.compile_model(problem, instance_data, {})
# get qubo model
pubo_builder = jmt.core.pubo.transpile_to_pubo(compiled_model=compiled_model, relax_method=jmt.core.pubo.RelaxationMethod.AugmentedLagrangian)
qubo, const = pubo_builder.get_qubo_dict(multipliers={})

OpenJijによる最適化計算の実行#

OpenJijを用いて計算してみます。

import openjij as oj

sampler = oj.SASampler()
response = sampler.sample_qubo(qubo, num_reads=1)

デコードと解の表示#

計算結果をデコードします。 ここでは、AAの中でA1A_1に分類されたindexとA0A_0に分類されたindexを分けて、それらについて和をとっています。

# decode a result to JijModeling sampleset
sampleset = jmt.core.pubo.decode_from_openjij(response, pubo_builder, compiled_model)
# get the indices of x == 1
class_1_indices = sampleset.record.solution["x"][0][0][0]
class_0_indices = [i for i in range(0, inst_N) if i not in class_1_indices]

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 : [ 1  3  5  6  7  9 10 11 12 16 18 19 20 22 23 25 27 29 34 36 37 40] , total value = 410
class 0 : [ 2  4  8 13 14 15 17 21 24 26 28 30 31 32 33 35 38 39] , total value = 410

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