整数変数と連立方程式の求解

このページでは，バイナリエンコーディングによる整数変数を説明します． 整数値をそのまま保持するネイティブ整数変数も利用できます．

整数変数

Hi-QUBOは整数変数をサポートしており,内部的には複数のバイナリ変数を用いて実装されています. 整数値の表現には従来のバイナリエンコーディングが使用されます. \(n\)個のバイナリ変数\(x_0, x_1, \ldots, x_{n-1}\)があるとします. これらの変数は,以下の線形式を用いて\(0\)から\(2^n-1\)までのすべての整数を表現できます:

\[ \begin{aligned} 2^0x_0+2^1x_1+\cdots 2^{n-1}x_{n-1} \end{aligned} \]

定数オフセット\(l\)を導入し,\(x_{n-1}\)の係数を任意の値\(d\)に置き換えると,次のようになります:

\[ \begin{aligned} l+2^0x_0+2^1x_1+\cdots +2^{n-2}x_{n-2}+dx_{n-1} \end{aligned} \]

この式は\(l\)から\(l+2^{n-1}+d-1\)までのすべての整数を表現できます. このエンコーディングに基づき,整数範囲が\([l,u]\)の変数は,以下を満たす適切な\(n\)と\(d\)（\(1\leq d\leq 2^{n-1}\)）を選ぶことで構成できます:

\[ \begin{aligned} u &= l+2^{n-1}+d-1 \end{aligned} \]

以下のHi-QUBOプログラムは整数変数の定義方法を示しています:

integer-variables-and-linear-systems-program1.cpp
#include <qbpp/qbpp.hpp>

int main() {
  auto x = 1 <= qbpp::var_int("x") <= 8;
  auto y = -10 <= qbpp::var_int("y") <= 10;
  std::cout << "x = " << x << " uses " << x.var_count() << " variables.\n";
  std::cout << "y = " << y << " uses " << y.var_count() << " variables.\n";
}

整数変数は範囲演算子 **<= <=**を用いて定義され,変数が取りうる整数の範囲を指定します. **qbpp::var_int("name")を範囲演算子と組み合わせた式全体が,指定されたnameを持つqbpp::Expr**オブジェクトを作成し,バイナリ変数でエンコードされた線形式を表現します. プログラムの出力は次のとおりです:

x = 1 +x[0] +2*x[1] +4*x[2] uses 3 variables.
y = -10 +y[0] +2*y[1] +4*y[2] +8*y[3] +5*y[4] uses 5 variables.

WARNING 整数変数に必要なバイナリ変数の数は,その範囲に対して対数的に増加します. \(u−l\)が大きい場合,QUBOのサイズが増大するため,広い整数範囲はできる限り避けるべきです.

連立方程式を解くためのQUBO定式化

Hi-QUBOは,変数を整数変数として表現することで連立方程式を解くことができます. 例として,解が\(x=4\),\(y=6\)である以下の方程式に対するQUBO定式化を構築します:

\[\begin{split} \begin{aligned} x + y = 10\\ 2x+4y = 28 \end{aligned} \end{split}\]

これらの方程式を解くために,範囲\([0,10]\)の整数変数\(x\)と\(y\)を定義し,それぞれ4つのバイナリ変数でエンコードします:

\[\begin{split} \begin{aligned} x = x_0 +2x_1 +4x_2 +3x_3\\ y = y_0 +2y_1 +4y_2 +3y_3 \end{aligned} \end{split}\]

以下の各ペナルティ式は,対応する方程式が満たされるとき,かつそのときに限り最小値0をとります:

\[\begin{split} \begin{aligned} f(x,y) &= (x+y-10)^2\\ &=(x_0 +2x_1 +4x_2 +3x_3+y_0 +2y_1 +4y_2 +3y_3-10)^2\\ g(x,y) &= (2x+4y -28)^2\\ &= (2\cdot(x_0 +2x_1 +4x_2 +3x_3)+4\cdot( y_0 +2y_1 +4y_2 +3y_3)-28)^2 \end{aligned} \end{split}\]

したがって,結合式

\[ \begin{aligned} h(x,y) &= f(x,y) +g(x,y) \end{aligned} \]

は,両方の方程式が同時に満たされるとき,正確にその最小値0を達成します.

Hi-QUBOプログラム

以下のHi-QUBOプログラムはQUBO式\(h(x,y)\)を構築し,それを解き,結果の\(x\)と\(y\)の値をデコードします:

integer-variables-and-linear-systems-program2.cpp
#include <qbpp/qbpp.hpp>
#include <qbpp/easy_solver.hpp>

int main() {
  auto x = 0 <= qbpp::var_int("x") <= 10;
  auto y = 0 <= qbpp::var_int("y") <= 10;

  auto f = x + y == 10;
  auto g = 2 * x + 4 * y == 28;
  auto h = f + g;
  h.simplify_as_binary();

  auto solver = qbpp::EasySolver(h);
  auto sol = solver.search({{"target_energy", 0}});

  std::cout << "sol = " << sol << std::endl;
  std::cout << "x = " << x << " = " << sol(x) << std::endl;
  std::cout << "y = " << y << " = " << sol(y) << std::endl;
  std::cout << "f = " << f << " = " << sol(f) << std::endl;
  std::cout << "g = " << g << " = " << sol(g) << std::endl;
  std::cout << "f.body() = " << f.body() << " = " << f.body(sol) << std::endl;
  std::cout << "g.body() = " << g.body() << " = " << g.body(sol) << std::endl;
}
```

{% endraw %}
まず、整数変数を表す`qbpp::Expr`オブジェクト**`x`**と**`y`**が範囲$[0,10]$で定義されます。
`qbpp::Expr`オブジェクト**`f`**は制約**`x + y == 10`**を表すために作成されます。
内部的には、これはQUBO式`qbpp::sqr(x + y -10)`と等価です。
同様に、**`g`**は制約**`2 * x + 4 * y == 28`**を表します。
結合式**`h = f + g`**は両方の方程式をエンコードします。
Easy Solverのインスタンスが`h`で作成され、最適解がすべての制約を満たすため、目標エネルギーが`0`に設定されます。
`search()`を呼び出すと、すべてのバイナリ変数の最適な割り当てを格納する`qbpp::Sol`オブジェクトsolが返されます。
最後に、プログラムは`sol`、`sol(x)`、`sol(y)`、`sol(f)`、`sol(g)`、`f.body(sol)`、`g.body(sol)`の値を出力します。
ここで、
- **`f`**: `x + y = 10`を強制するペナルティ式。したがって、方程式が満たされるとき、かつそのときに限り`sol(f) = 0`となります。
- **`f.body()`**: 線形式`x + y`。したがって`f.body(sol)`は`x + y`の実際の評価値を返します。

**`g`**と**`g.body()`**についても同様です。

プログラムの出力結果は次のとおりです:

{% raw %}
```
sol = 0:{{x[0],0},{x[1],1},{x[2],1},{x[3],0},{y[0],0},{y[1],0},{y[2],1},{y[3],0}}
x = x[0] +2*x[1] +4*x[2] +3*x[3] = 6
y = y[0] +2*y[1] +4*y[2] +3*y[3] = 4
f = 100 -19*x[0] -36*x[1] -64*x[2] -51*x[3] -19*y[0] -36*y[1] -64*y[2] -51*y[3] +4*x[0]*x[1] +8*x[0]*x[2] +6*x[0]*x[3] +2*x[0]*y[0] +4*x[0]*y[1] +8*x[0]*y[2] +6*x[0]*y[3] +16*x[1]*x[2] +12*x[1]*x[3] +4*x[1]*y[0] +8*x[1]*y[1] +16*x[1]*y[2] +12*x[1]*y[3] +24*x[2]*x[3] +8*x[2]*y[0] +16*x[2]*y[1] +32*x[2]*y[2] +24*x[2]*y[3] +6*x[3]*y[0] +12*x[3]*y[1] +24*x[3]*y[2] +18*x[3]*y[3] +4*y[0]*y[1] +8*y[0]*y[2] +6*y[0]*y[3] +16*y[1]*y[2] +12*y[1]*y[3] +24*y[2]*y[3] = 0
g = 784 -108*x[0] -208*x[1] -384*x[2] -300*x[3] -208*y[0] -384*y[1] -640*y[2] -528*y[3] +16*x[0]*x[1] +32*x[0]*x[2] +24*x[0]*x[3] +16*x[0]*y[0] +32*x[0]*y[1] +64*x[0]*y[2] +48*x[0]*y[3] +64*x[1]*x[2] +48*x[1]*x[3] +32*x[1]*y[0] +64*x[1]*y[1] +128*x[1]*y[2] +96*x[1]*y[3] +96*x[2]*x[3] +64*x[2]*y[0] +128*x[2]*y[1] +256*x[2]*y[2] +192*x[2]*y[3] +48*x[3]*y[0] +96*x[3]*y[1] +192*x[3]*y[2] +144*x[3]*y[3] +64*y[0]*y[1] +128*y[0]*y[2] +96*y[0]*y[3] +256*y[1]*y[2] +192*y[1]*y[3] +384*y[2]*y[3] = 0
f.body() = x[0] +2*x[1] +4*x[2] +3*x[3] +y[0] +2*y[1] +4*y[2] +3*y[3] = 10
g.body() = 2*x[0] +4*x[1] +8*x[2] +6*x[3] +4*y[0] +8*y[1] +16*y[2] +12*y[3] = 28

まず、整数変数を表すqbpp::Exprオブジェクト**xとyが範囲\([0,10]\)で定義されます。 qbpp::Exprオブジェクトfは制約x + y == 10を表すために作成されます。 内部的には、これはQUBO式qbpp::sqr(x + y -10)と等価です。 同様に、gは制約2 * x + 4 * y == 28を表します。 結合式h = f + g**は両方の方程式をエンコードします。 Easy Solverのインスタンスがhで作成され、最適解がすべての制約を満たすため、目標エネルギーが0に設定されます。 search()を呼び出すと、すべてのバイナリ変数の最適な割り当てを格納するqbpp::Solオブジェクトsolが返されます。 最後に、プログラムはsol、sol(x)、sol(y)、sol(f)、sol(g)、f.body(sol)、g.body(sol)の値を出力します。 ここで、

  • f: x + y = 10を強制するペナルティ式。したがって、方程式が満たされるとき、かつそのときに限りsol(f) = 0となります。

  • f.body(): 線形式x + y。したがってf.body(sol)はx + yの実際の評価値を返します。

**gとg.body()**についても同様です。

プログラムの出力結果は次のとおりです:

これにより,x,y,および制約式f,g,f.body(),g.body()の値が解と整合していることが確認できます.

WARNING Hi-QUBOは,左辺が式で右辺が整数の場合にのみ==演算子をサポートしています. 整数==式や式==式の形式の比較はサポートされていません. 詳細は比較演算子で説明しています.

このページでは，バイナリエンコーディングによる整数変数を説明します． 整数値をそのまま保持するネイティブ整数変数も利用できます．

整数変数

PyQBPPは整数変数をサポートしており,内部的には複数のバイナリ変数を用いて実装されています. 整数値の表現には従来のバイナリエンコーディングが使用されます. \(n\)個のバイナリ変数\(x_0, x_1, \ldots, x_{n-1}\)があるとします. これらの変数は,以下の線形式を用いて\(0\)から\(2^n-1\)までのすべての整数を表現できます:

\[ \begin{aligned} 2^0x_0+2^1x_1+\cdots 2^{n-1}x_{n-1} \end{aligned} \]

定数オフセット\(l\)を導入し,\(x_{n-1}\)の係数を任意の値\(d\)に置き換えると,次のようになります:

\[ \begin{aligned} l+2^0x_0+2^1x_1+\cdots +2^{n-2}x_{n-2}+dx_{n-1} \end{aligned} \]

この式は\(l\)から\(l+2^{n-1}+d-1\)までのすべての整数を表現できます. このエンコーディングに基づき,整数範囲が\([l,u]\)の変数は,以下を満たす適切な\(n\)と\(d\)（\(1\leq d\leq 2^{n-1}\)）を選ぶことで構成できます:

\[ \begin{aligned} u &= l+2^{n-1}+d-1 \end{aligned} \]

以下のプログラムは整数変数の定義方法を示しています.

integer-variables-and-linear-systems-program1.py
import pyqbpp as qbpp

x = qbpp.var("x", between=(1, 8))
y = qbpp.var("y", between=(-10, 10))
print(f"x = {x} uses {x.var_count} variables.")
print(f"y = {y} uses {y.var_count} variables.")

整数変数は between= キーワード引数を使って定義し,変数がとりうる整数範囲を指定します. qbpp.var("name", between=(min, max)) は,指定された name を持つ整数変数 Expr オブジェクトを作成し,バイナリ変数でエンコードされた線形式を表現します. プログラムの出力は以下の通りです.

x = 1 +x[0] +2*x[1] +4*x[2] uses 3 variables.
y = -10 +y[0] +2*y[1] +4*y[2] +8*y[3] +5*y[4] uses 5 variables.

注意 整数変数に必要なバイナリ変数の数は,その範囲に対して対数的に増加します. max - min が大きい場合,QUBOのサイズが増大するため,広い整数範囲はできる限り避けるべきです.

連立方程式を解くためのQUBO定式化

PyQBPPは,変数を整数変数として表現することで連立方程式を解くことができます. 例として,解が \(x=4\),\(y=6\) である以下の方程式に対するQUBO定式化を構築します.

\[\begin{split} \begin{aligned} x + y = 10\\ 2x+4y = 28 \end{aligned} \end{split}\]

これらの方程式を解くために,範囲 \([0,10]\) の整数変数 \(x\) と \(y\) を定義し,それぞれ4つのバイナリ変数でエンコードします:

\[\begin{split} \begin{aligned} x = x_0 +2x_1 +4x_2 +3x_3\\ y = y_0 +2y_1 +4y_2 +3y_3 \end{aligned} \end{split}\]

以下の各ペナルティ式は,対応する方程式が満たされるとき,かつそのときに限り最小値0をとります:

\[\begin{split} \begin{aligned} f(x,y) &= (x+y-10)^2\\ &=(x_0 +2x_1 +4x_2 +3x_3+y_0 +2y_1 +4y_2 +3y_3-10)^2\\ g(x,y) &= (2x+4y -28)^2\\ &= (2\cdot(x_0 +2x_1 +4x_2 +3x_3)+4\cdot( y_0 +2y_1 +4y_2 +3y_3)-28)^2 \end{aligned} \end{split}\]

したがって,結合式

\[ \begin{aligned} h(x,y) &= f(x,y) +g(x,y) \end{aligned} \]

は,両方の方程式が同時に満たされるとき,正確にその最小値0を達成します.

PyQBPP プログラム

以下のプログラムはQUBO式 \(h(x,y)\) を構築し,それを解いて結果の \(x\) と \(y\) の値をデコードします.

integer-variables-and-linear-systems-program2.py
import pyqbpp as qbpp

x = qbpp.var("x", between=(0, 10))
y = qbpp.var("y", between=(0, 10))

f = qbpp.constrain(x + y, equal=10)
g = qbpp.constrain(2 * x + 4 * y, equal=28)
h = f + g
h.simplify_as_binary()

solver = qbpp.EasySolver(h)
sol = solver.search(target_energy=0)

print("sol =", sol)
print("x =", x, "=", sol(x))
print("y =", y, "=", sol(y))
print("f =", f, "=", sol(f))
print("g =", g, "=", sol(g))
print("f.body =", f.body, "=", sol(f.body))
print("g.body =", g.body, "=", sol(g.body))
```

{% endraw %}
まず、整数変数 **`x`** と **`y`** を範囲 $[0,10]$ で定義します。
式 **`f`** は制約 **`qbpp.constrain(x + y, equal=10)`** を表すために作成されます。
内部的には、これはQUBO式 `qbpp.sqr(x + y - 10)` と等価です。
同様に、**`g`** は制約 **`qbpp.constrain(2 * x + 4 * y, equal=28)`** を表します。
結合式 **`h = f + g`** は両方の方程式をエンコードします。
Easy Solver のインスタンスが `h` で作成され、最適解がすべての制約を満たすため、目標エネルギーが `0` に設定されます。
`search()` を呼び出すと、すべてのバイナリ変数の最適な割り当てを格納する解オブジェクト `sol` が返されます。
最後に、プログラムは `sol`、`sol(x)`、`sol(y)`、`sol(f)`、`sol(g)`、`sol(f.body)`、`sol(g.body)` の値を出力します。
ここで、
- **`f`**: `x + y = 10` を強制するペナルティ式。したがって、方程式が満たされるとき、かつそのときに限り `sol(f) = 0` となります。
- **`f.body`**: 線形式 `x + y`。したがって `sol(f.body)` は `x + y` の実際の評価値を返します。

**`g`** と **`g.body`** についても同様です。

プログラムの出力結果は次のとおりです。

{% raw %}
```
sol = Sol(energy=0, {x[0]: 0, x[1]: 1, x[2]: 1, x[3]: 0, y[0]: 0, y[1]: 0, y[2]: 1, y[3]: 0})
x = x[0] +2*x[1] +4*x[2] +3*x[3] = 6
y = y[0] +2*y[1] +4*y[2] +3*y[3] = 4
f = 100 -19*x[0] -36*x[1] -64*x[2] -51*x[3] -19*y[0] -36*y[1] -64*y[2] -51*y[3] +4*x[0]*x[1] +8*x[0]*x[2] +6*x[0]*x[3] +2*x[0]*y[0] +4*x[0]*y[1] +8*x[0]*y[2] +6*x[0]*y[3] +16*x[1]*x[2] +12*x[1]*x[3] +4*x[1]*y[0] +8*x[1]*y[1] +16*x[1]*y[2] +12*x[1]*y[3] +24*x[2]*x[3] +8*x[2]*y[0] +16*x[2]*y[1] +32*x[2]*y[2] +24*x[2]*y[3] +6*x[3]*y[0] +12*x[3]*y[1] +24*x[3]*y[2] +18*x[3]*y[3] +4*y[0]*y[1] +8*y[0]*y[2] +6*y[0]*y[3] +16*y[1]*y[2] +12*y[1]*y[3] +24*y[2]*y[3] = 0
g = 784 -108*x[0] -208*x[1] -384*x[2] -300*x[3] -208*y[0] -384*y[1] -640*y[2] -528*y[3] +16*x[0]*x[1] +32*x[0]*x[2] +24*x[0]*x[3] +16*x[0]*y[0] +32*x[0]*y[1] +64*x[0]*y[2] +48*x[0]*y[3] +64*x[1]*x[2] +48*x[1]*x[3] +32*x[1]*y[0] +64*x[1]*y[1] +128*x[1]*y[2] +96*x[1]*y[3] +96*x[2]*x[3] +64*x[2]*y[0] +128*x[2]*y[1] +256*x[2]*y[2] +192*x[2]*y[3] +48*x[3]*y[0] +96*x[3]*y[1] +192*x[3]*y[2] +144*x[3]*y[3] +64*y[0]*y[1] +128*y[0]*y[2] +96*y[0]*y[3] +256*y[1]*y[2] +192*y[1]*y[3] +384*y[2]*y[3] = 0
f.body = x[0] +2*x[1] +4*x[2] +3*x[3] +y[0] +2*y[1] +4*y[2] +3*y[3] = 10
g.body = 2*x[0] +4*x[1] +8*x[2] +6*x[3] +4*y[0] +8*y[1] +16*y[2] +12*y[3] = 28

まず、整数変数 x と y を範囲 \([0,10]\) で定義します。 式 f は制約 qbpp.constrain(x + y, equal=10) を表すために作成されます。 内部的には、これはQUBO式 qbpp.sqr(x + y - 10) と等価です。 同様に、g は制約 qbpp.constrain(2 * x + 4 * y, equal=28) を表します。 結合式 h = f + g は両方の方程式をエンコードします。 Easy Solver のインスタンスが h で作成され、最適解がすべての制約を満たすため、目標エネルギーが 0 に設定されます。 search() を呼び出すと、すべてのバイナリ変数の最適な割り当てを格納する解オブジェクト sol が返されます。 最後に、プログラムは sol、sol(x)、sol(y)、sol(f)、sol(g)、sol(f.body)、sol(g.body) の値を出力します。 ここで、

  • f: x + y = 10 を強制するペナルティ式。したがって、方程式が満たされるとき、かつそのときに限り sol(f) = 0 となります。

  • f.body: 線形式 x + y。したがって sol(f.body) は x + y の実際の評価値を返します。

g と g.body についても同様です。

プログラムの出力結果は次のとおりです。

これにより,x,y,および制約式 f,g,f.body,g.body の値が解と整合していることが確認できます.

注意 PyQBPPの qbpp.constrain(expr, equal=n) の n には,整数のほか式（Var/Term/Expr）も指定できます. penalty は sqr(expr - n),body は expr になります（例: qbpp.constrain(x + y, equal=2*z)）. 演算子形式の 式 == 式 は直接サポートされていないため,式同士の等価制約には equal= か qbpp.constrain(expr1 - expr2, equal=0) を使ってください. 詳細は比較制約で説明しています.