変数配列を用いた分割問題の解法¶

分割問題

\(w=(w_0, w_1, \ldots, w_{n-1})\) を \(n\) 個の正の数とします。 分割問題は、これらの数を2つの集合 \(P\) と \(Q\) (\(=\overline{P}\)) に分割し、2つの集合の要素の和ができるだけ近くなるようにする問題です。より具体的には、次を最小化する部分集合 \(L \subseteq \lbrace 0,1,\ldots, n-1\rbrace\) を見つける問題です：

\[\begin{split} \begin{aligned} P(L) &= \sum_{i\in L}w_i \\ Q(L) &= \sum_{i\not\in L}w_i \\ f(L) &= \left| P(L)-Q(L) \right| \end{aligned} \end{split}\]

この問題はQUBO問題として定式化できます。 \(x=(x_0, x_1, \ldots, x_{n-1})\) を集合 \(L\) を表すバイナリ変数とし、\(i\in L\) は \(x_i=1\) のときかつそのときに限ります。 \(x\) を用いて \(P(L)\)、\(Q(L)\)、\(f(L)\) を次のように書き換えることができます：

\[\begin{split} \begin{aligned} P(x) &= \sum_{i=0}^{n-1} w_ix_i \\ Q(x) &= \sum_{i=0}^{n-1} w_i \overline{x_i} \\ f(x) &= \left( P(x)-Q(x) \right)^2 \end{aligned} \end{split}\]

ここで \(\overline{x_i}\) は \(x_i\) の否定リテラルであり、\(x_i=0\) のとき \(1\)、\(x_i=1\) のとき \(0\) を取ります。数学的には \(\overline{x_i} = 1 - x_i\) ですが、Hi-QUBOは否定リテラルを ~ 演算子（例: ~x[i]）でネイティブに扱えるため、\(1 - x_i\) への展開を避けて効率的に処理できます。詳細は**否定リテラル**を参照してください。

明らかに \(f(x)=f(L)^2\) が成り立ちます。関数 \(f(x)\) は \(x\) の2次式であり、\(f(x)\) を最小化する最適解は元の分割問題の最適解も与えます。

分割問題のHi-QUBOプログラム

以下のHi-QUBOプログラムは、固定された8個の数の分割問題のQUBO定式化を作成し、Exhaustive Solverを用いて解を求めます。

integer-partition-program1.cpp¶

#include <qbpp/qbpp.hpp>
#include <qbpp/exhaustive_solver.hpp>

int main() {
  std::vector<uint32_t> w = {64, 27, 47, 74, 12, 83, 63, 40};

  auto x = qbpp::var("x", w.size());
  auto p = qbpp::toExpr(0);
  auto q = qbpp::toExpr(0);
  for (size_t i = 0; i < w.size(); ++i) {
    p += w[i] * x[i];
    q += w[i] * ~x[i];
  }
  auto f = qbpp::sqr(p - q);
  std::cout << "f = " << f.simplify_as_binary() << std::endl;

  auto solver = qbpp::ExhaustiveSolver(f);
  auto sol = solver.search();

  std::cout << "Solution: " << sol << std::endl;
  std::cout << "f(sol) = " << f(sol) << std::endl;
  std::cout << "p(sol) = " << p(sol) << std::endl;
  std::cout << "q(sol) = " << q(sol) << std::endl;

  std::cout << "P :";
  for (size_t i = 0; i < w.size(); ++i) {
    if (sol(x[i]) == 1) {
      std::cout << " " << w[i];
    }
  }
  std::cout << std::endl;

  std::cout << "Q :";
  for (size_t i = 0; i < w.size(); ++i) {
    if (sol(x[i]) == 0) {
      std::cout << " " << w[i];
    }
  }
  std::cout << std::endl;
}

このプログラムでは、w は8個の数を持つ std::vector オブジェクトとして定義されています。 w.size()=8 個のバイナリ変数の配列 x が定義されます。 2つの qbpp::Expr オブジェクト p と q が定義され、\(P(x)\) と \(Q(x)\) の式がforループ内で構築されます。ここで ~x[i] は x[i] の否定リテラル \(\overline{x_i}\) を表します。 qbpp::Expr オブジェクト f は \(f(x)\) の式を格納します。

f に対するExhaustive Solverオブジェクト solver が作成され、search() メンバ関数の呼び出しにより解 sol（qbpp::Sol オブジェクト）が得られます。

\(f(x)\)、\(P(x)\)、\(Q(x)\) の値はそれぞれ f(sol)、p(sol)、q(sol) の呼び出しにより評価されます。集合 \(L\) と \(\overline{L}\) の数はforループを使って表示されます。これらのループでは、sol(x[i]) が sol における x[i] の値を返します。

このプログラムの出力は以下のとおりです：

f = 168100 -88576*x[0] -41364*x[1] -68244*x[2] -99456*x[3] -19104*x[4] -108564*x[5] -87444*x[6] -59200*x[7] +13824*x[0]*x[1] +24064*x[0]*x[2] +37888*x[0]*x[3] +6144*x[0]*x[4] +42496*x[0]*x[5] +32256*x[0]*x[6] +20480*x[0]*x[7] +10152*x[1]*x[2] +15984*x[1]*x[3] +2592*x[1]*x[4] +17928*x[1]*x[5] +13608*x[1]*x[6] +8640*x[1]*x[7] +27824*x[2]*x[3] +4512*x[2]*x[4] +31208*x[2]*x[5] +23688*x[2]*x[6] +15040*x[2]*x[7] +7104*x[3]*x[4] +49136*x[3]*x[5] +37296*x[3]*x[6] +23680*x[3]*x[7] +7968*x[4]*x[5] +6048*x[4]*x[6] +3840*x[4]*x[7] +41832*x[5]*x[6] +26560*x[5]*x[7] +20160*x[6]*x[7]
Solution: 0:{{x[0],0},{x[1],0},{x[2],1},{x[3],0},{x[4],1},{x[5],1},{x[6],1},{x[7],0}}
f(sol) = 0
p(sol) = 205
q(sol) = 205
P : 47 12 83 63
Q : 64 27 74 40

注意 qbpp::Expr オブジェクト f と qbpp::Sol オブジェクト sol に対して、f(sol) と sol(f) はどちらも sol 上で f を評価した値を返します。同様に、qbpp::Var オブジェクト a に対して、a(sol) と sol(a) はどちらも解 sol における a の値を返します。 f(sol) の形式は、関数をある点で評価することに対応するため、数学的な観点からは自然です。一方、sol(f) は解オブジェクトが式を評価するという、オブジェクト指向プログラミングの観点からは自然です。好みに応じてどちらの形式を使用しても構いません。

配列演算を用いたHi-QUBOプログラム

Hi-QUBOにはコードを簡潔にするための豊富な配列演算があります。以下のコードでは、w は qbpp::array() で定義し、x は qbpp::var() でバイナリ変数の配列として定義しています。配列に対してオーバーロードされた演算子 * は要素ごとの積を返すため、qbpp::sum(w * x) は \(P(L)\) を表す qbpp::Expr オブジェクトを返します。変数の配列に対する ~ 演算子は否定リテラルの配列を返します。したがって、qbpp::sum(w * ~x) は \(Q(L)\) を格納する qbpp::Expr オブジェクトを返します。

  auto w = qbpp::array({64, 27, 47, 74, 12, 83, 63, 40});
  auto x = qbpp::var("x", w.size());
  auto p = qbpp::sum(w * x);
  auto q = qbpp::sum(w * ~x);
  auto f = qbpp::sqr(p - q);

これらの配列演算を使用することで、Hi-QUBOプログラムを簡潔にできます。さらに、大きな配列に対する配列演算はマルチスレッドにより並列化されるため、QUBOモデルの作成プロセスを高速化できます。

注意演算子 +、-、* は、2つの配列オブジェクト間、およびスカラーと配列オブジェクト間の両方でオーバーロードされています。 2つの配列オブジェクトの場合、オーバーロードされた演算子は要素ごとの演算を行います。スカラーと配列オブジェクトの場合、オーバーロードされた演算子は配列の各要素にスカラー演算を適用します。

この問題を非線形関数 abs やネイティブ制約 cons() でも定式化し、3 通りの定式化を比較するケーススタディは、Python 版（ページ上部の言語切替で Python を選択）を参照してください。

PyQBPP では、同じ問題を従来の 2 乗ペナルティ（sqr）、 非線形関数（abs）、ネイティブ制約（cons()）のいずれでも定式化できます。このケーススタディでは分割問題を題材に 3 通りの定式化を並べて書き、エネルギーの意味・モデルの大きさ・得られる情報の違いを比較します。非線形関数とネイティブ制約そのものの説明は非線形関数とネイティブ制約を参照してください。

分割問題は、\(n\) 個の正の数 \(w_0, w_1, \ldots, w_{n-1}\) を 2 つの集合に分割し、両集合の和 \(P\) と \(Q\) ができるだけ近くなるようにする問題です（詳細は分割問題）。バイナリ変数 \(x_i\) を用いると、2 つの和は

\[ P = \sum_{i=0}^{n-1} w_i x_i, \qquad Q = \sum_{i=0}^{n-1} w_i \overline{x_i} \]

と書けます。ここで \(\overline{x_i}\) は \(x_i\) の否定リテラルです。以下の 3 つのプログラムは、式 p・q の構築までは共通で、目的関数 f を作る 1 行だけが異なります。

定式化 1: 2 乗ペナルティ sqr

差の 2 乗 \((P - Q)^2\) を最小化する、従来の QUBO 定式化です:

partition-formulations-program1.py¶

import pyqbpp as qbpp

w = qbpp.array([64, 27, 47, 74, 12, 83, 63, 40])
x = qbpp.var("x", shape=len(w))
p = qbpp.sum(w * x)
q = qbpp.sum(w * ~x)
f = qbpp.sqr(p - q)
f.simplify_as_binary()
solver = qbpp.ExhaustiveSolver(f)
sol = solver.search()
print(f"(P - Q)^2 = {sol.energy}")
print(f"P = {sol(p)}, Q = {sol(q)}")

プログラムの出力は以下の通りです:

(P - Q)^2 = 0
P = 205, Q = 205

simplify_as_binary() により、f は 2 乗が展開された 2 次の QUBO 多項式になります（展開後の式は分割問題のページに示されています）。展開後は通常の 2 次式なので、Hi-QUBO のすべてのソルバーはもちろん、ネイティブ制約に対応しない外部の QUBO ツールでもそのまま扱える、最も汎用的な形です。一方、エネルギーは差の 2 乗であり、差そのものを知るには p・q を評価する必要があります。

定式化 2: 絶対値 abs

非線形関数 qbpp.abs を使うと、分割問題の目的関数 \(|P - Q|\) をそのまま書けます:

partition-formulations-program2.py¶

import pyqbpp as qbpp

w = qbpp.array([64, 27, 47, 74, 12, 83, 63, 40])
x = qbpp.var("x", shape=len(w))
p = qbpp.sum(w * x)
q = qbpp.sum(w * ~x)
f = qbpp.abs(p - q)
f.simplify_as_binary()
solver = qbpp.ExhaustiveSolver(f)
sol = solver.search()
print(f"|P - Q| = {sol.energy}")
print(f"P = {sol(p)}, Q = {sol(q)}")

プログラムの出力は以下の通りです:

|P - Q| = 0
P = 205, Q = 205

エネルギーが差 \(|P - Q|\) そのものになり、解の質を直接読み取れます。2 乗の展開は起きず、線形式 \(P - Q\) は展開されないまま関数の本体として扱われます。 abs を含む式は、Hi-QUBO にバンドルされているソルバーが関数値を直接扱って探索します（対応ソルバーの一覧は非線形関数とネイティブ制約を参照）。

定式化 3: ネイティブ制約 cons

「\(P - Q = 0\)（等分割）」を制約として宣言する書き方です。目的関数はなく、モデルは制約だけからなります:

partition-formulations-program3.py¶

import pyqbpp as qbpp

w = qbpp.array([64, 27, 47, 74, 12, 83, 63, 40])
x = qbpp.var("x", shape=len(w))
p = qbpp.sum(w * x)
q = qbpp.sum(w * ~x)
f = qbpp.cons(p - q == 0)
f.simplify_as_binary()
solver = qbpp.ExhaustiveSolver(f)
sol = solver.search()
print(f"Energy = {sol.energy}")
print(f"P = {sol(p)}, Q = {sol(q)}")
print(f"violated constraints = {f.cons(sol)}")

プログラムの出力は以下の通りです:

Energy = 0
P = 205, Q = 205
violated constraints = 0

cons() で宣言した等式制約の値は違反量の 2 乗 \((P - Q)^2\) なので、エネルギーの値は定式化 1 と同じです。違いは意味論です — 式は制約として宣言されているため、違反本数を返す f.cons(sol) や制約ごとの違反量を報告する violations() が使え、 target_energy は「エネルギーが target 以下かつ全制約充足」のときに探索を停止します。大きなインスタンスを EasySolver で解くときに target_energy に 0 を指定すると、完全分割が見つかった時点で探索を打ち切れます。

モデルの大きさの比較

3 通りの定式化でソルバーに渡されるモデルの大きさを比べてみます。 sol.info の var_count・term_count が、それぞれモデルの変数数と目的関数の多項式の項数です:

partition-formulations-program4.py¶

import pyqbpp as qbpp

w = qbpp.array([64, 27, 47, 74, 12, 83, 63, 40])
x = qbpp.var("x", shape=len(w))
p = qbpp.sum(w * x)
q = qbpp.sum(w * ~x)

f1 = qbpp.sqr(p - q)
f1.simplify_as_binary()
s1 = qbpp.ExhaustiveSolver(f1).search()

f2 = qbpp.abs(p - q)
f2.simplify_as_binary()
s2 = qbpp.ExhaustiveSolver(f2).search()

f3 = qbpp.cons(p - q == 0)
f3.simplify_as_binary()
s3 = qbpp.ExhaustiveSolver(f3).search()

print(f"sqr : var_count = {s1.info['var_count']}, term_count = {s1.info['term_count']}")
print(f"abs : var_count = {s2.info['var_count']}, term_count = {s2.info['term_count']}")
print(f"cons: var_count = {s3.info['var_count']}, term_count = {s3.info['term_count']}")

プログラムの出力は以下の通りです:

sqr : var_count = 8, term_count = 36
abs : var_count = 8, term_count = 0
cons: var_count = 8, term_count = 0

sqr では 2 乗の展開により、8 個の 1 次項と \(\binom{8}{2} = 28\) 個の 2 次項、計 36 項の QUBO 多項式になります。 abs・cons では展開は起きません — 目的関数の多項式は空（term_count = 0）で、8 項の線形式 \(P - Q\) が関数・制約の本体としてそのまま保持されます。

一般に \(n\) 個の数では、sqr の項数は \(n(n + 1)/2\) に増えます — \(n = 1000\) なら約 50 万項です。abs・cons の本体は \(n\) 項の線形式のままなので、モデルの構築も解 1 つの評価も軽く、インスタンスが大きいほど差が広がります。

完全分割が存在しない場合

3 通りのエネルギーの意味の違いは、完全分割（\(P = Q\) となる分割）が存在しないインスタンスではっきり現れます。次のプログラムは、最良の分割でも差が 5 になる 8 個の数で 3 通りの定式化を解きます:

partition-formulations-program5.py¶

import pyqbpp as qbpp

w = qbpp.array([17, 18, 39, 61, 84, 85, 88, 89])
x = qbpp.var("x", shape=len(w))
p = qbpp.sum(w * x)
q = qbpp.sum(w * ~x)

f1 = qbpp.sqr(p - q)
f1.simplify_as_binary()
print(f"sqr : {qbpp.ExhaustiveSolver(f1).search().energy}")

f2 = qbpp.abs(p - q)
f2.simplify_as_binary()
print(f"abs : {qbpp.ExhaustiveSolver(f2).search().energy}")

f3 = qbpp.cons(p - q == 0)
f3.simplify_as_binary()
sol = qbpp.ExhaustiveSolver(f3).search()
print(f"cons: {sol.energy}, violated constraints = {f3.cons(sol)}")
print(f"P = {sol(p)}, Q = {sol(q)}")

プログラムの出力は以下の通りです:

sqr : 25
abs : 5
cons: 25, violated constraints = 1
P = 243, Q = 238

最適な分割は \(P = 243\)、\(Q = 238\) で、差は 5 です。 sqr と cons のエネルギーは \(5^2 = 25\)、abs のエネルギーは 5 と、それぞれの定義どおりの値になります。cons では制約 \(P - Q = 0\) を満たす分割が存在しないため違反本数が 1 と報告されますが、ソルバーはエラーにはならず、違反量が最小の解（= 最良の分割）を返します。制約が満たせないときは最も制約に近い解が得られる — これがネイティブ制約のソフトな意味論です。

まとめ

定式化	エネルギー	モデル	特徴
`qbpp.sqr(p - q)`	\((P - Q)^2\)	2 次式に展開（\(n(n+1)/2\) 項）	最も汎用的。外部の QUBO ツールでも扱える
`qbpp.abs(p - q)`	\(\lvert P - Q \rvert\)	展開なし（\(n\) 項の本体）	エネルギーが差そのもの
`qbpp.cons(p - q == 0)`	\((P - Q)^2\)	展開なし（\(n\) 項の本体）	制約として宣言。違反の集計・`target_energy` に参加

どの定式化でも、得られる最適な分割そのものは同じです。目的の量（差）を直接エネルギーとして読みたいなら abs、「等分割」という制約を宣言して充足を追跡したいなら cons()、ネイティブ制約に対応しない外部の QUBO ツールに渡す必要があるなら sqr が適しています。対応ソルバーの一覧と詳細な演算規則は非線形関数とネイティブ制約を参照してください。不等式制約についての同様の比較は、ナップサック問題を参照してください。