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

整数変数について

QUBO++ は整数変数をサポートしており、これらは内部的には複数のバイナリ変数を用いて実装されています。

PyQBPP は整数変数をサポートしており、これらは内部的には複数のバイナリ変数を用いて実装されています。

整数値の表現には、従来のバイナリエンコーディングが使用されます。

\(n\) 個のバイナリ変数 \(x_0,\cdots,x_{n-1}\) があるとします。
これらの変数を用いることで、次の線形表現により \(0\) から \(2^n-1\) までのすべての整数を表すことができます:

\[ 2^0x_0+\cdots+2^{n-2}x_{n-2}+2^{n-1}x_{n-1} \]

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

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

この式は、\(l\) から \(l+2^{n-1}+d-1\) までのすべての整数を表すことができます。 このエンコーディングに基づき、整数範囲 \([l, u]\) を持つ変数は、\(n\) および \(d\)\(1 \leq d \leq 2^{n-1}\) )を適切に選ぶことで、次を満たすように構築できます:

\[ u=l+2^{n-1}+d-1 \]

次の 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") 関数は、指定した名前を持つ qbpp::VarInt オブジェクトを生成します。
このオブジェクトは、バイナリ変数によってエンコードされた線形式を表します。

このプログラムは、次の式を出力します:

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.

次の PyQBPP プログラムは、整数変数の定義方法を示しています。

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

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

整数変数はbetween() 関数を使って定義し、変数がとりうる整数範囲を指定します。
この演算子により、その変数が取り得る整数範囲を指定します。

var_int("name") 関数は、指定した名前を持つ VarIntCore オブジェクトを生成し、between(var_int("name"), min, max) はバイナリ変数でエンコードされた線形式を表す VarInt オブジェクトを作成します。 このプログラムは、次の式を出力します:

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を用いて連立方程式を解く

QUBO++ は、変数を整数変数として表現することで、連立方程式を解くことができます。

PyQBPP は、変数を整数変数として表現することで、連立方程式を解くことができます。

例として、解が \(x=4\) および \(y=6\) である次の方程式を考えます:

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

これらの方程式を解くために、整数範囲 \([0, 10]\) を持つ整数変数 \(x\) および \(y\) を定義します。
それぞれは 4 個のバイナリ変数によってエンコードされます:

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

次の各ペナルティ式は、対応する方程式が満たされる場合に限り、最小値 0 を取ります:

\[\begin{split} 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(x_0+2x_1+4x_2+3x_3)+4(y_0+2y_1+4y_2+3y_3)-28)^2. \end{split}\]

したがって、これらを組み合わせた式

\[ h(x, y)=f(x, y)+g(x, y) \]

は、両方の方程式が同時に満たされる場合に正確に最小値 0 を達成します。

QUBO++ program

次の 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;
}

まず、範囲 \([0, 10]\) を持つ qbpp::VarInt オブジェクト x および y を定義します。 qbpp::Expr オブジェクト f は制約 x + y == 10 を表します。
内部的には、これは QUBO 式 qbpp::sqr(x + y - 10) と等価です。

同様に、g は制約 2 * x + 4 * y == 28 を表します。
結合式 h = f + g は、両方の方程式を同時にエンコードします。

h を用いて Easy Solver のインスタンスを生成し、最適解はすべての制約を満たすため、目標エネルギーを 0 に設定します。
search() を呼び出すと、すべてのバイナリ変数の最適な割り当てを格納した qbpp::Sol オブジェクト sol が返されます。

最後に、solsol(x)sol(y)sol(f)sol(g)sol(*f)sol(*g) の値を出力します。

ここで:

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

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

g および *g についても同様です。

このプログラムは、次の結果を出力します:

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 = x[0] +2*x[1] +4*x[2] +3*x[3] +y[0] +2*y[1] +4*y[2] +3*y[3] = 10
*g = 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

これにより、xy、および制約式 fg*f*g の値が、解と整合していることを確認できます。

PyQBPP program

次の 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))

まず、範囲 \([0, 10]\) を持つ VarInt オブジェクト x および y を定義します。 Expr オブジェクト f は制約 x + y == 10 を表します。
内部的には、これは QUBO 式 sqr(x + y - 10) と等価です。

同様に、g は制約 2 * x + 4 * y == 28 を表します。
結合式 h = f + g は、両方の方程式を同時にエンコードします。

h を用いて Easy Solver のインスタンスを生成し、最適解はすべての制約を満たすため、目標エネルギーを 0 に設定します。
search() を呼び出すと、すべてのバイナリ変数の最適な割り当てを格納した Sol オブジェクト sol が返されます。

最後に、solsol(x)sol(y)sol(f)sol(g)sol(*f)sol(*g) の値を出力します。

ここで:

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

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

g および g.body についても同様です。

このプログラムは、次の結果を出力します:

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 = x[0] +2*x[1] +4*x[2] +3*x[3] +y[0] +2*y[1] +4*y[2] +3*y[3] = 10
*g = 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

これにより、xy、および制約式 fg*f*g の値が、解と整合していることを確認できます。


WARNING

QUBO++ では、== 演算子は「左辺が式、右辺が整数」の場合にのみサポートされています。

PyQBPP では、== 演算子は「左辺が式、右辺が整数」の場合にのみサポートされています。

integer == expressionexpression == expression の形式の比較はサポートされていません。

詳細は 比較演算子 を参照してください。