数分割問題

数分割問題#

JijModelingとJijZeptSolverを用いて、数分割問題を解く方法をご紹介します。 この問題はLucas, 2014, "Ising formulations of many NP problems"の 2.1. Number Partitioningでも触れられています。

数分割問題とは#

数の集合を2つに分割したときに、それらの集合の数の和が等しくなるように分割することを考えるのが、数分割問題です。 単純な例として、\(A = \{1, 2, 3, 4\}\)のような数の集合を分割することを考えましょう。 これはすぐに\(\{1, 4\}, \{2, 3\}\)のように分割され、これらの集合の和は両方とも5であることがわかります。 このように、数の集合のサイズが小さなときは、簡単に解を発見することができます。 しかし、数の集合サイズが大きくなると、その解の発見は困難となります。 このような理由から、このチュートリアルでは任意の数の集合サイズを解けるような実装を考えていきましょう。

数理モデル#

最初に、この問題のハミルトニアンをモデル化しましょう。 \(A\)を分割したい数の集合とし、この\(A\)は\(a_i \ (i = \{0, 1, \dots, N-1\})\)の要素を持つとします。 ここで\(N\)は\(A\)の要素数です。 そして\(A\)を\(A_0, A_1\)の2つの集合に分けるとします。 \(x_i = 0\)のとき\(a_i\)を\(A_0\)に入れ、\(x_i = 1\)のとき\(a_i\)を\(A_1\)に含めるバイナリ変数\(x_i\)を定義します。 これを用いると、\(A_0\)の数の総和は\(\sum_i a_i (1-x_i)\)となり、\(A_1\)の数の総和は\(\sum_i a_i x_i\)となります。 これら2つの総和が等しいという制約条件を満たす、解を探す必要があります。

\[ \sum_{i=0}^{N-1} a_i (1-x_i) = \sum_{i=0}^{N-1} a_i x_i \ \Longrightarrow \ \sum_{i=0}^{N-1}a_i (2 x_i - 1) = 0 \tag{1} \]

式(1)にペナルティ法を適用すれば、数分割問題のハミルトニアンを得ることができます。

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

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で定義します。 最後に、(2)式で用いられている添字をiで定義します。
次に、数分割問題のハミルトニアンを実装しましょう。

problem += (jm.sum(N, lambda i: a[i]*(2*x[i] - 1))) ** 2

Jupyter Notebook上であれば、実装された数理モデルを確認することができます。

problem
\[\begin{split}\begin{array}{rl} \text{Problem}\colon &\text{number partition}\\\displaystyle \min &\displaystyle {\left(\sum _{i=0}^{N-1}{{a}_{i}\cdot \left(2\cdot {x}_{i}-1\right)}\right)}^{2}\\&\\\text{where}&\\&\text{Decision Variables:}\\&\qquad \begin{alignedat}{2}x&\in \mathop{\mathrm{Array}}\left[N;\left\{0, 1\right\}\right]&\quad &1\text{-dim binary variable}\\\end{alignedat}\\&\\&\text{Placeholders:}\\&\qquad \begin{alignedat}{2}a&\in \mathop{\mathrm{Array}}\left[(-);\mathbb{R}\right]&\quad &1\text{-dimensional array of placeholders with elements in }\mathbb{R}\\\end{alignedat}\\&\\&\text{Dependent Variables:}\\&\qquad \begin{alignedat}{2}N&=\mathop{\mathtt{len\_{}at}}\left(a,0\right)&\quad &\in \mathbb{N}\\\end{alignedat}\end{array} \end{split}\]

インスタンスの準備#

数の集合\(A\)を準備します。 ここでは\(1 \sim 40\)の整数を分割することにしましょう。 \(N_i\)から\(N_f\)までの連続した整数かつ項数が偶数の場合、2つに分割された各集合の総和は

\[ (\mathrm{total \ value}) = \frac{(N_f + N_i) (N_f - N_i + 1)}{4} \]

のように計算することができます。 \(1 \sim 40\)では、この値は410となります。

import numpy as np

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

JijZeptSolverで解く#

jijzept_solverを用いて、この数分割問題を解いてみましょう。

import jijzept_solver

instance = problem.eval(instance_data)
solution = jijzept_solver.solve(instance, solve_limit_sec=1.0)

解の可視化#

それでは、得られた結果を確認してみましょう。 得られたバイナリ変数から、\(A_0, A_1\)に数を分割します。 最後に、それらの総和を計算しましょう。

df = solution.decision_variables_df
x_1_indices = np.ravel(df[df["value"]==1.0]["subscripts"].to_list())
x_0_indices = np.ravel(df[df["value"]==0.0]["subscripts"].to_list())
A_1 = instance_data["a"][x_1_indices]
A_0 = instance_data["a"][x_0_indices]
print(f"class 1 : {A_1} , total value = {np.sum(A_1)}")
print(f"class 0 : {A_0} , total value = {np.sum(A_0)}")

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

予想通り、これら2つの集合の各総和が410となりました。

目次