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

整数変数

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}\]

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

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{aligned} x + y = 10\\ 2x+4y = 28 \end{aligned}\]

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

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

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

\[\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}\]

したがって、結合式

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

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

PyQBPP プログラム

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

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

まず、整数変数 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 についても同様です。

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

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

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

注意 PyQBPPの qbpp.constrain(expr, equal=n) は、expr が式で n が整数の場合にのみサポートされています。 式 == 式 の形式の等価制約は直接サポートされていません。qbpp.constrain(expr1 - expr2, equal=0) と書き換えてください。 詳細は比較制約で説明しています。