Constraint

このセクションでは制約条件を表すクラスについて説明します。

参考

このセクションを読む前に Tutorial にて Amplify SDK の基本的な使用方法と Polynomial にて多項式オブジェクトの使用方法を確認してください。

制約条件オブジェクト

量子アニーリングマシン・イジングマシンで問題を解く際に、使用する変数の間に制約条件を課す必要が生じることがあります。

例えば、\(x_0, x_1, x_2\) を変数とする目的関数 \(f\left(x_0, x_1, x_2 \right)\) の最適解を探索する際に、変数の間に \(x_0 + x_1 + x_2 = 1\) などの等式制約や、\(x_0 + x_1 + x_2 \leq 2\) などの不等式制約が課されていると、制約条件を満たす「実行可能解」の中から最適解を見つける必要があります。

ところが、量子アニーリングマシン・イジングマシンの扱える QUBO やイジングモデルは、そのような制約条件を直接は扱うことが出来ません。そのため、基本的なアプローチとして、制約条件を満たす場合に最小値をとるようなペナルティ関数 \(g\) に重み \(\lambda\) を与えて元の目的関数 \(f\) に追加する方法が行われます。

つまり、\(f\) の代わりに、\(f' = f + \lambda g \quad (\lambda > 0)\) の最適解を求めることで、ペナルティ関数 \(g\) が最小、すなわち制約条件を満たす実行可能解の取得が可能になります。実際には、マシンの出力が必ずしも制約条件を満たす解であるとは限らないことや、重み \(\lambda\) が適切であるかは解らないので、\(f'\) の解を \(g\) で評価したときに最小値であるかを確認することで、実行可能解かどうかを識別します。

Amplify SDK では、このようなペナルティ関数の定式化、ペナルティ関数の重み、解を評価する際の充足チェック等、構築から管理までを行う制約条件オブジェクトという形でこれを抽象化します。制約条件オブジェクトの管理する構成要素は次の通りです。

  • 制約条件式: 制約条件を表現する式

  • ペナルティ関数: 制約条件に相当するペナルティ関数

  • 重み: ペナルティ関数に与える重み (上記における \(\lambda\))

制約条件オブジェクトは、BinaryQuadraticModel 等に与えることで論理模型に変換できるほか、多項式クラスや行列クラスと組み合わせることで、与えられた多項式・行列オブジェクトに対する付加的な制約条件として扱われます。最終的に、論理模型から物理模型に渡されるタイミングにおいて、「重み」 \(\times\) 「ペナルティ関数」が目的関数に足され、イジングマシンへの入力となります。

多項式クラス・行列クラスとの対応関係は次の通りです。

多項式クラス

行列クラス

制約条件クラス

BinaryPoly

BinaryMatrix

BinaryConstraint

IsingPoly

IsingMatrix

IsingConstraint

BinaryIntPoly

BinaryIntMatrix

BinaryIntConstraint

IsingIntPoly

IsingIntPoly

IsingIntConstraint

Amplify SDK では、典型的な制約条件に関して、制約条件オブジェクトの作成を行うヘルパー関数が用意されています。

  • equal_to() は与えられた多項式 \(f\) が定数値 \(c\) になるように制約します。one_hot() は \(f = 1\) の場合におけるエイリアスとして与えられます。

  • less_equal() は与えられた多項式 \(f\) が定数値 \(c\) 以下になるように制約します。ただし、制約範囲と \(f\) の値域が重なっている必要があります。

  • greater_equal() は与えられた多項式 \(f\) が定数値 \(c\) 以上になるように制約します。ただし、制約範囲と \(f\) の値域が重なっている必要があります。

  • clamp() は与えられた多項式 \(f\) が定数値 \(c_1\) 以上 \(c_2\) 以下になるように制約します。ただし、制約範囲と \(f\) の値域が重なっている必要があります。

制約

生成関数

使用可能条件

\(f \left( \mathbf{q} \right) = c\)

equal_to()

\(\exists \mathbf{q}, f\left( \mathbf{q} \right) = c\)

\(f \left( \mathbf{q} \right) \le c\)

less_equal()

\(\exists \mathbf{q}, f\left( \mathbf{q} \right) \le c\)

\(f \left( \mathbf{q} \right) \ge c\)

greater_equal()

\(\exists \mathbf{q}, f\left( \mathbf{q} \right) \ge c\)

\(c_1 \le f \left( \mathbf{q} \right) \le c_2\)

clamp()

\(\exists \mathbf{q}, f\left( \mathbf{q} \right) \le c_2\) かつ \(\exists \mathbf{q}, f\left( \mathbf{q} \right) \ge c_1\)

これらのヘルパー関数は、与えられた制約条件式そのものと制約条件に対応するペナルティ関数を持ち、重み \(1\) で初期化された制約条件オブジェクトを返却します。それぞれのヘルパー関数の使い方の詳細については、以下で解説します。

注釈

ヘルパー関数で与えられるペナルティ関数は、制約条件を破った場合のペナルティの最低値として \(1\) が課せられるように規格化されています。目的関数の取り得る値に応じて、制約条件オブジェクトに適切な重みを与える必要があることに注意してください。

等式制約

等式制約を実現するために equal_to() 関数は、与えられた任意の多項式 \(f\) が定数値 \(c\) になるように制約します。次の例では \(3 q_0 - q_1 - 2 q_2 + q_3 = 2\) を表す制約条件を作成し、これを満たす解を求めます。

from amplify import BinarySymbolGenerator
from amplify.constraint import equal_to

gen = BinarySymbolGenerator()  # BinaryPoly の変数ジェネレータを宣言
q = gen.array(4)  # 長さ 4 の変数配列を生成
f = 3 * q[0] - q[1] - 2 * q[2] + q[3]  # 制約条件の左辺 f
g = equal_to(f, 2)  # f = 2 を表す制約条件オブジェクトを作成

>>> g
3 q_0 - q_1 - 2 q_2 + q_3 == 2

equal_to() 関数には第一引数に制約対象を与え、第二引数に定数値を設定します。次のようにクライアントを設定し、作成した制約条件オブジェクトから論理模型を構築することで、実行結果を確認できます。

from amplify import BinaryQuadraticModel, Solver
from amplify.client import FixstarsClient

client = FixstarsClient()
client.token = "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx"
client.parameters.timeout = 1000
client.parameters.outputs.duplicate = True  # エネルギー値が同一の解を重複して出力する
client.parameters.outputs.num_outputs = 0

model = BinaryQuadraticModel(g)
solver = Solver(client)
result = solver.solve(model)

>>> for solution in result:
>>>     values = solution.values
>>>     print(f"q = {q.decode(values)}, f = {f.replace_all(values)}")
q = [1. 1. 0. 0.], f = 2.0
q = [1. 0. 1. 1.], f = 2.0

確かに \(f = 2\) を満たす解が得られていることがわかります。

制約対象として多項式配列を与えることも可能です。この場合は多項式配列に対して、暗黙的に sum() メソッドが呼ばれ、配列の要素和が制約対象と解釈されます。

from amplify import BinarySymbolGenerator
from amplify.constraint import equal_to

gen = BinarySymbolGenerator()  # BinaryPoly の変数ジェネレータを宣言
q = gen.array(4)  # 長さ 4 の変数配列を生成

>>> equal_to(q, 2)
q_0 + q_1 + q_2 + q_3 == 2
>>> equal_to(q.sum(), 2)
q_0 + q_1 + q_2 + q_3 == 2

注釈

等式制約に対応する equal_to() 関数は、与えられた多項式 \(f\) と定数 \(c\) に対して \(\min f\) を見積もり、次のようにペナルティ関数 \(g\) を生成します。

  1. \(\min f < c\) の場合:

\[g = \left( f - c \right)^2\]

  1. \(\min f = c\) の場合:

\[g = f - c\]

ペナルティ関数 \(g\) は \(f = c\) の時に \(g = 0\) となり、それ以外の時に \(g > 0\) となります。

巡回セールスマン問題の定式化

等式制約で表現される代表的な制約条件として One-hot 制約があげられます。複数の変数に対し、一つだけ \(1\) になり他は \(0\) とする制約です。複数の変数を用いて多値を表現する際に現れます。

ここでは、巡回セールスマン問題の定式化 に現れる One-hot 制約を例に制約条件の定式化を説明します。

巡回セールスマン問題の目的関数は次の通りです。

\[f = \sum_{n=0}^{N-1}{\sum_{i=0}^{N-1}{\sum_{j=0}^{N-1}{ d_{ij} q_{n, i} q_{n+1, j} }}} \quad \left( q_{N,j} = q_{0,j} \right)\]

ここで、\(N\) は問題サイズ (都市数) を表します。\(q\) は目的関数の決定変数でバイナリ変数とします。また、\(d\) は都市間の距離を表す距離行列とします。 目的関数の定式化は 多項式配列を用いた目的関数の構築 にて説明しました。最も洗練された書き方は次の通りです。

import numpy as np
from amplify import BinarySymbolGenerator, sum_poly

# 問題サイズ (都市数)
ncity = 10

# 距離行列を乱数から作成
d = np.triu(np.random.rand(ncity, ncity), k=1)
d += np.triu(d).T  # 対称行列に変換

# 変数配列の作成
gen = BinarySymbolGenerator()
q = gen.array(ncity, ncity)  # ncity x ncity の変数配列を生成

# 目的関数の定式化
qr = q.roll(-1, axis=0)  # q を列方向 (axis=0) に -1 ずらす
f = sum_poly("ij,ni,nj->", d, q, qr)  # アインシュタイン縮約規則

巡回セールスマン問題の定式化 は、目的関数だけでは不十分です。変数配列に対して次のような制約条件が課せられます。これは変数の配列をどの都市を何番目に訪問するという順列として解釈するために必要です。

\[ \begin{align}\begin{aligned}\sum_{i=0}^{N-1}{q_{n, i}} = 1 \quad n \in \left\{0, 1, \cdots, N - 1 \right\}\\\sum_{n=0}^{N-1}{q_{n, i}} = 1 \quad i \in \left\{0, 1, \cdots, N - 1 \right\}\end{aligned}\end{align} \]

上記の制約条件は二次元変数配列のそれぞれの行・列において、\(1\) となる変数が必ず一つだけ存在することを表します。これは次のように等式制約関数を用いて表現されます。

from amplify.constraint import equal_to

g = sum(equal_to(q[i, :].sum(), 1) + equal_to(q[:, i].sum(), 1) for i in range(ncity))

各行と各列に対して和を sum() メソッドで計算し、equal_to() 関数を呼びます。複数の制約条件オブジェクトは足し合わせることが出来るので、sum() 関数にジェネレータ式を渡します。

equal_to() 関数などの制約条件オブジェクトの生成関数は、多項式配列を直接与えることも出来ます。この場合、多項式配列に対して sum() メソッドを呼んだ場合と等価です。また、equal_to() 関数の特別な場合として、One-hot 制約に表す one_hot() 関数が用意されています。これらを用いると、次のように簡略化されます。

from amplify.constraint import one_hot

g = sum(one_hot(q[i]) + one_hot(q[:, i]) for i in range(ncity))

次に制約条件オブジェクトの重みを決めます。デフォルトで生成される制約条件オブジェクトは制約違反の場合に最小のペナルティ値として \(1\) が課せられるように規格化されています。これは目的関数の大きさによっては、過小あるいは過大なことがあります。重みが小さすぎる場合には、最適化した結果が制約条件を満たしにくくなり、大きすぎる場合には、目的関数の最小化が不十分な解になってしまうことがあります。

適切な制約条件の重みを事前に求めることは一般に難しいですが、次のような指針が参考になります。

制約を違反した場合の目的関数の最大利得 \(<\) 制約を違反した場合の最小ペナルティ値

マシンの出力する (準) 最適解が事前にわかっていないため、上記の左辺を事前に正確に見積もることは一般に難しいです。そこで大雑把に左辺の最大値を見積ることで、ペナルティの重みを決めることにします。

巡回セールスマン問題の場合において制約を違反した場合の利得とは次のように見積もられます。まず、全ての制約を満たす実行可能解が得られたとします。その解を改変することで目的関数が最大でどれだけ得をするかを考えます。目的関数の値は距離を加算しなければ小さく出来るので、例えば都市 \(i \rightarrow j \rightarrow k\) という訪問順において、辺 \(i \rightarrow j\) と \(j \rightarrow k\) を削除する、あるいは \(j\) をスキップして \(i \rightarrow k\) とする変更が考えられます。このような変更は少なくとも一つの制約条件が違反となるためペナルティが課せられます。一方で、全ての実行可能解に対する辺の削除や変更による目的関数の最大利得は、辺の最大長 \(\max_{(i,j)} d_{ij}\) で上限を見積もることが可能です。詳細な議論は省きますが、結果としてペナルティ重みは次のようにすれば十分であることがわかります。

g *= np.max(d)

制約条件オブジェクトには定数をかけることで重みを変更できます。この機能の詳細は後述する 制約条件の重みの設定 を参照してください。目的関数の定式化と合わせて、巡回セールスマン問題の定式化 をまとめると次の通りです。

import numpy as np
from amplify import BinarySymbolGenerator, sum_poly
from amplify.constraint import one_hot

# 問題サイズ (都市数)
ncity = 10

# 距離行列を作成
d = np.triu(np.random.rand(ncity, ncity), k=1)
d += np.triu(d).T  # 対称行列に変換

# 変数配列の作成
gen = BinarySymbolGenerator()
q = gen.array(ncity, ncity)  # ncity x ncity の変数配列を生成

# 目的関数の定式化
qr = q.roll(-1, axis=0)  # q を列方向 (axis=0) に -1 ずらす
f = sum_poly("ij,ni,nj->", d, q, qr)  # sum_poly の代わりに einsum を用いても同様

# 制約条件の定式化
g = np.max(d) * sum(one_hot(q[i]) + one_hot(q[:, i]) for i in range(ncity))

# 論理模型の構築
m = f + g

参考

制約条件オブジェクトの重みの設定は大規模な問題ほど重要になります。制約を満たす結果を得るため、またはより良い解を得るためには、重みの調整に対する試行錯誤が必要になることがあります。上記で述べた制約条件の重みの見積もりや調整による解精度の向上については下記を参照してください。

K. Takehara, D. Oku, Y. Matsuda, S. Tanaka and N. Togawa, "A Multiple Coefficients Trial Method to Solve Combinatorial Optimization Problems for Simulated-annealing-based Ising Machines," 2019 IEEE 9th International Conference on Consumer Electronics (ICCE-Berlin), Berlin, Germany, 2019, pp. 64-69, doi: 10.1109/ICCE-Berlin47944.2019.8966167.

不等式制約

次の生成関数は、整数値をとる多項式に対して、不等式(以上・以下・範囲)で制約します。

  • less_equal() は与えられた多項式 \(f\) が値 \(c\) 以下 (\(f \le c\)) になるように制約します。

  • greater_equal() は与えられた多項式 \(f\) が値 \(c\) 以上 (\(f \ge c\)) になるように制約します。

  • clamp() は与えられた多項式 \(f\) が値 \(c_1\) 以上 \(c_2\) 以下 (\(c_1 \le f \le c_2\)) になるように制約します。

less_equal() 関数には第一引数に制約対象を与え、第二引数に定数値を設定します。greater_equal() 関数も同様です。clamp() 関数には、第二引数と第三引数それぞれに最小値と最大値を与えます。

次の例では、整数係数の多項式に対する制約として \(4 q_0 + 3 q_1 + 2 q_2 + q_3 \le 3\) を表す制約条件を作成し、これを満たす解を求めます。

from amplify import BinarySymbolGenerator
from amplify.constraint import less_equal

gen = BinarySymbolGenerator()  # BinaryPoly の変数ジェネレータを宣言
q = gen.array(4)  # 長さ 4 の変数配列を生成
f = 4 * q[0] + 3 * q[1] + 2 * q[2] + q[3]
g = less_equal(f, 3) # f <= 3 を表す制約条件オブジェクトを作成

>>> g
4 q_0 + 3 q_1 + 2 q_2 + q_3 <= 3

次のようにクライアントを設定し、作成した制約条件オブジェクトから論理模型を構築することで、実行結果を確認できます。

from amplify import BinaryQuadraticModel, Solver
from amplify.client import FixstarsClient

client = FixstarsClient()
client.token = "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx"
client.parameters.timeout = 1000
client.parameters.outputs.duplicate = True  # エネルギー値が同一の解を重複して出力する
client.parameters.outputs.num_outputs = 0

model = BinaryQuadraticModel(g)
solver = Solver(client)
result = solver.solve(model)

>>> for solution in result:
>>>     values = solution.values
>>>     print(f"q = {q.decode(values)}, f = {f.replace_all(values)}")
q = [0. 0. 0. 1.], f = 1.0
q = [0. 0. 1. 0.], f = 2.0
q = [0. 0. 0. 0.], f = 0.0
q = [0. 0. 1. 1.], f = 3.0
q = [0. 1. 0. 0.], f = 3.0

多項式 \(f\) について、\(f \le 3\) を満たす解が得られました。

下記は greater_equal() 関数と clamp() 関数の使用例です。

from amplify.constraint import greater_equal

g = greater_equal(f, 3) # f >= 3 を表す制約条件オブジェクトを作成
result = solver.solve(g)

>>> g
4 q_0 + 3 q_1 + 2 q_2 + q_3 >= 3

>>> for solution in result:
>>>     values = solution.values
>>>     print(f"q = {q.decode(values)}, f = {f.replace_all(values)}")
q = [1. 1. 0. 1.], f = 8.0
q = [1. 0. 0. 1.], f = 5.0
q = [0. 1. 1. 1.], f = 6.0
q = [1. 1. 0. 0.], f = 7.0
q = [1. 0. 1. 0.], f = 6.0
q = [1. 1. 1. 0.], f = 9.0
q = [1. 0. 1. 1.], f = 7.0
q = [1. 1. 1. 1.], f = 10.0
q = [1. 0. 0. 0.], f = 4.0
q = [0. 1. 1. 0.], f = 5.0
q = [0. 1. 0. 1.], f = 4.0
q = [0. 1. 0. 0.], f = 3.0
q = [0. 0. 1. 1.], f = 3.0

from amplify.constraint import clamp

g = clamp(f, 2, 3)  # 2 <= f <= 3 を表す制約条件オブジェクトを作成
result = solver.solve(g)

>>> g
2 <= 4 q_0 + 3 q_1 + 2 q_2 + q_3 <= 3

>>> for solution in result:
>>>     values = solution.values
>>>     print(f"q = {q.decode(values)}, f = {f.replace_all(values)}")
q = [0. 0. 1. 0.], f = 2.0
q = [0. 0. 1. 1.], f = 3.0
q = [0. 1. 0. 0.], f = 3.0

制約対象として多項式配列を与えると、equal_to() 関数と同様に配列の要素の和が制約対象の多項式として解釈されます。

from amplify import BinarySymbolGenerator
from amplify.constraint import less_equal

gen = BinarySymbolGenerator()  # BinaryPoly の変数ジェネレータを宣言
q = gen.array(4)  # 長さ 4 の変数配列を生成

>>> less_equal(q, 2)
q_0 + q_1 + q_2 + q_3 <= 2
>>> less_equal(q.sum(), 2)
q_0 + q_1 + q_2 + q_3 <= 2

注釈

制約式が整数係数多項式でない場合、不等式制約を正確に表現できないことがあります。 また、制約式が整数係数の場合であっても、係数が大きすぎる場合、イジングマシンでは良い解が得られない傾向にあります。 制約式の両辺に適切な値を掛けてから制約条件オブジェクトを作成することで、これらが改善される可能性があります。

不等式制約の実現方法

Amplify SDK には、不等式制約のペナルティ関数を生成する方法が大きく分けて 2 種類実装されています。 1 つは補助変数を発行して不等式制約を等式制約に変換する方法 (整数エンコーディング) で、 もう 1 つは制約を満たしているときに比較的小さな値を取る関数を見つけて、それをペナルティ関数とする方法 (緩和法) です。

これらの方法の詳細に関しては次節以降で述べますが、それぞれ以下のような特徴があります。

アルゴリズム

補助変数の発行

ペナルティ関数の次数

厳密性

整数エンコーディング (制約式が整数係数)

有

制約式の 2 乗

o

整数エンコーディング (制約式が実数係数)

有

制約式の 2 乗

x

緩和法 (制約式が整数係数)

無

制約式と同じまたは 2 乗

x

緩和法 (制約式が実数係数)

無

制約式と同じまたは 2 乗

x

補助変数をただちに発行しない方法であっても、ペナルティ関数が 3 次以上となる場合は次数下げのために追加で補助変数が発行されます。 また、上表で厳密性が x となっているものについては、制約条件の重みをどれだけ調整しても論理模型の最適解と問題の最適解を一致させることができない場合があります。

デフォルトでは、制約式の次数や係数の値に関わらず、常に整数エンコーディングが選択されます。

整数エンコーディング

greater_equal(), less_equal(), clamp() 関数は、デフォルトでは等式制約と同様のペナルティ関数 \(g\) で実現されます。

まず、制約条件を満たす \(f\) がとりうる整数範囲 \(\left[c_1, c_2 \right]\) を、制約対象の多項式 \(f\) の値域と与えられた制約範囲から決定します。 ただし、\(c_1, c_2\) が整数でない場合は、四捨五入されて整数に丸められます。 \(c_2 - c_1\) の値によって次の二通りの方法でペナルティ関数 \(g\) が構成されます。

  1. \(c_2 - c_1 > 1\) の場合

    制約対象の整数範囲 \(\left[c_1, c_2 \right]\) が値域となる多項式 \(h\) を用いて、次のようにペナルティ関数 \(g\) を作成します。

    \[g = \left( f - h \right)^2\]

    ペナルティ関数 \(g\) は \(f = h\) つまり \(c_1 \le f \le c_2\) の時に \(g = 0\) となり、それ以外の時に \(g > 0\) となります。

  2. \(c_2 - c_1 = 1\) の場合

\(c_1\) と \(c_2\) が隣接した整数である場合、次のようにペナルティ関数 \(g\) を作成します。

\[g = \left( f - c_1 \right)\left( f - c_2 \right)\]

ペナルティ関数 \(g\) は \(f = c_1\) または \(f = c_2\) の場合に \(g = 0\) となり、それ以外の時に \(g > 0\) となります。

前者の \(c_2 - c_1 > 1\) の場合は、整数範囲 \(\left[c_1, c_2 \right]\) が値域となる多項式 \(h\) を実現するために、整数値をバイナリ多項式 (またはイジング多項式) で表すための整数エンコーディングが必要であることを意味します。

Amplify SDK では多項式 \(h\) を作成するための整数エンコーディングとして次の方法が指定できます。

方法

多項式

補助変数の数

係数の絶対値の最大

列挙型名

備考

Unary法

\(c_1 + \sum_{i=0}^{c_2 - c_1 - 1}{q_i}\)

\(c_2 - c_1\)

1

Unary

Binary法

\(c_1 + \sum_{i=0}^{\lfloor \log_2 \left( c_2 - c_1 + 1 \right) \rfloor - 1}{2^i q_i} + \cdots\)

\(O \left( \log_2 \left( c_2 - c_1 + 1 \right) \right)\)

\(O \left( c_2 - c_1 \right)\)

Binary

イジング多項式では未実装

Linear法

\(c_1 + \sum_{i=1}^{d := \lfloor \sqrt{c_2 - c_1} \rfloor}{ i q_i } + \sum_{i=1}^{d - 1}{ i q_{i+d} } + \cdots\)

\(O \left( \sqrt{c_2 - c_1} \right)\)

\(O \left( \sqrt{c_2 - c_1} \right)\)

Linear

イジング多項式では未実装

Unary 法以外の方法では、上記の定式化において端数が生じることがあります。Binary法では端数 \(r\) について \(r q_x\) を追加します。Linear法では端数 \(r\) について \(r\) を2等分した値を2つ追加します。ただし \(r\) が奇数の場合は \(\left( r + 1 \right)/2\) および \(\left( r - 1 \right)/2\) が追加されます。

不等式制約を生成する greater_equal(), less_equal(), clamp() 関数でキーワード引数 method に InequalityFormulation を与えることで整数エンコーディング法を指定できます。整数エンコーディングには表現できる整数範囲に応じた補助変数が必要になります。補助変数の管理は Amplify SDK が自動で行いますが、どの程度の補助変数が必要になり係数値の絶対値の大きさがどの程度なのかは留意する必要があるかもしれません。補助変数の数が多い、または係数値の絶対値が (目的関数と比較して) 大きいと求解精度に影響を与える可能性があるためです。

次の実行例は Unary 法での不等式制約オブジェクトの構築例です。

from amplify import BinarySymbolGenerator, InequalityFormulation
from amplify.constraint import less_equal

gen = BinarySymbolGenerator()  # BinaryPoly の変数ジェネレータを宣言
q = gen.array(5)  # 長さ 5 の変数配列を生成
f = q[0] + 2 * q[1] + 3 * q[3] + 4 * q[4] + 5 * q[5]

>>> less_equal(f, 10, method=InequalityFormulation.Unary)
q_0 + 2 q_1 + 3 q_2 + 4 q_3 + 5 q_4 <= 10

エンコーディング手法が未指定または Default の場合は、次の指針に基づいて決定されます。

  1. Unary 法、Linear 法、Binary 法のうち必要な補助変数が最も少ない方法を選ぶ

  2. 必要な補助変数の数が同値の場合は、優先順位を Unary 法、Linear 法、Binary 法の順に選ぶ

制約式が整数値以外を取りうる場合や制約範囲が整数でない場合は、生成されるペナルティ関数 \(g\) は制約が満たされていても必ずしも最小値を取らないことに注意してください。

緩和法

制約式が高次であったり取りうる値の範囲が大きい場合や、制約式が整数係数でない場合に、上記の整数エンコーディング法では定式化が難しいことがあります。そのため、Amplify SDK では制約式の違反量に応じたペナルティを与える緩和法によって不等式制約を実現する方法が用意されています。

等式制約や整数エンコーディング法による不等式制約も、制約を満たすときのみペナルティ関数が0となるような緩和法に当てはまります。一方、ここでの緩和法は制約を満たしていても必ずしもペナルティが0とならないことから、一種の近似的手法になります。

まず、制約対象の多項式 \(f\) の値域と与えられた制約範囲から \(f\) の制約範囲 \(\left[c_1, c_2 \right]\) を決定します。その上で、制約 \(f \le c_2\) について次の二通りの方法でペナルティ関数 \(g\) が作成可能です。

  1. 制約式に対する1次のペナルティ関数

    \[g = \left(f - c_1 \right) / (c_2 - c_1)\]

  2. 制約式に対する2次のペナルティ関数

\[g = \left(f - m \right)^2 / (c_2 - m)^2, \quad m = \left(c_1 + c_2 \right) / 2\]

どちらの方法でも上限値 \(f = c_2\) の場合にペナルティが \(g = 1\) になるように規格化し、制約式の違反量に応じてペナルティが大きくなるように定式化されます。なお、制約 \(f \ge c_1\) については、内部的に \(-f \le -c_1\) に変換してペナルティ関数を作成します。

上記の緩和法を使用するには greater_equal(), less_equal() 関数の method 引数に RelaxationLinear, RelaxationQuadra または Relaxation を指定します。RelaxationLinear は制約式に対する1次のペナルティ関数、RelaxationQuadra は2次のペナルティ関数に対応します。イジングマシンでは3次以上の多項式は次数下げが必要になることを考慮して、Relaxation は制約対象の多項式が1次なら後者を、2次以上なら前者を自動選択します。

注釈

緩和法においては、ペナルティ関数の係数 (ラグランジュ乗数) の調整がこれまでの方法以上に重要となります。重みを変更しながら何度も求解し、実行可能解の中から最適解を探索するという方針が考えられます。

参考

ペナルティ関数の重みの指定方法については 制約条件の重みの設定 を参照してください。

ナップサック問題の定式化

不等式制約を用いた例としてナップサック問題の定式化を取り上げます。ナップサック問題とは、いくつかのアイテムとそれぞれのアイテムの重量と価値が与えられた場合に、重量の合計が \(W\) を超えないように、価値の合計が最大となるアイテムを取り出す問題です。

具体的な問題例として、4 アイテムの重量と価値と重量制限が次のように与えられるとします。

w = [55, 44, 26, 30]  # 重量
v = [50, 36, 24, 33]  # 価値
maxW = 100  # 最大の合計重量

アイテムを選ぶかどうかを表すバイナリ変数配列を生成し、価値の合計を最大化する目的関数を作成します。

from amplify import BinarySymbolGenerator

gen = BinarySymbolGenerator()  # BinaryPoly の変数ジェネレータを宣言
q = gen.array(4)  # 長さ 4 の変数配列を生成

# 目的関数 (符号を反転して最小化問題にする)
f = -(v * q).sum()

一方、重量の合計が maxW を超えないために制約条件を設定する必要があります。

from amplify.constraint import less_equal

# 重さ maxW 以下制約
g = less_equal((w * q).sum(), maxW)

from amplify import BinarySymbolGenerator, Solver
from amplify.constraint import less_equal
from amplify.client import FixstarsClient

# 価値を最大化
f = -(v * q).sum()

# 重さ100以下制約
maxW = 100
g = less_equal((w * q).sum(), maxW)

# モデル
m = f + g

client = FixstarsClient()
client.parameters.timeout = 1000

solver = Solver(client)
result = solver.solve(m)

次のようにして結果を確認できます。

>>> print(f"価値の合計: {-result[0].energy}")
>>> print(f"重量の合計: {(w * q.decode(result[0].values)).sum()}")
価値の合計: 93.0
重量の合計: 100.0

ペナルティ関数を用いた制約

構築したい制約条件がヘルパー関数の適用範囲外であったり、独自のペナルティ関数を設定したい場合には、penalty() 関数を用いることで制約条件オブジェクトを独自に作成できます。

独自の制約条件オブジェクトを作成する例として、次のような制約条件を考えます。ある2変数 \(q_0, q_1\) について、両方同時に \(1\) にならないような制約 (\(q_0 q_1 = 0\)) を与えたいとします。この制約条件は equal_to() でも実現できますが、あえてペナルティ関数を設定すると次のように与えられます。

\[g = q_0 q_1\]

ペナルティ関数 \(g\) は \(q_0 = q_1 = 1\) の時だけ \(g = 1\) となるペナルティが課せられます。その他の場合は \(g = 0\) です。

次のように、penalty() 関数を用いて制約条件オブジェクトを作成します。

from amplify import BinarySymbolGenerator
from amplify.constraint import penalty

gen = BinarySymbolGenerator()  # BinaryPoly の変数ジェネレータを宣言
q = gen.array(2)  # 長さ 2 の変数配列を生成
g = penalty(q[0] * q[1])  # ペナルティ関数 q[0] * q[1] を設定

>>> g
q_0 q_1 == 0

これまでと同様に、次のようにクライアントを設定し、作成した制約条件オブジェクトから論理模型を構築することで、実行結果を確認できます。

from amplify import BinaryQuadraticModel, Solver
from amplify.client import FixstarsClient

client = FixstarsClient()
client.token = "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx"
client.parameters.timeout = 1000
client.parameters.outputs.duplicate = True  # エネルギー値が同一の解を重複して出力する
client.parameters.outputs.num_outputs = 0

model = BinaryQuadraticModel(g)
solver = Solver(client)
result = solver.solve(model)

>>> for solution in result:
>>>     values = solution.values
>>>     energy = solution.energy
>>>     print(f"q = {q.decode(values)}, energy = {energy}")
q = [0. 0.], energy = 0.0
q = [0. 1.], energy = 0.0
q = [1. 0.], energy = 0.0

\(q_0 = q_1 = 1\) 以外の解が得られていることがわかります。この問題の場合は、ペナルティ関数のみから目的関数を構築したので、制約条件を満たす解のエネルギー値は \(0\) となります。

高度な制約条件の構築

penalty() 関数に与えるペナルティ関数 \(g\) は、全ての変数の組合せに対して \(g \ge 0\) であり \(g = 0\) の場合に制約が満たされているとみなされます。この充足判定の式を変更し、任意の値による等式や不等式に変更することが可能です。

等式の充足判定の指定では、充足条件を \(g = 0\) から任意の値 \(g = m\) に変更が可能です。例えばペナルティ関数 \(g\) の最小値 \(m\) が \(m \ne 0\) である場合、次のように等式を表すキーワード引数 eq に右辺を与えます。

from amplify import BinarySymbolGenerator
from amplify.constraint import penalty

gen = BinarySymbolGenerator()
q = gen.array(3) # 長さ 3 の変数配列を生成
g = q[0] * q[1] + q[1] * q[2] + q[2] * q[0] + 1 # ペナルティ関数 g (最小値 1)
p = penalty(g, eq = 1) # g = 1 で充足判定を行う

>>> p
q_0 q_1 + q_0 q_2 + q_1 q_2 + 1 == 1

一方で、不等式による充足判定の指定は、ペナルティ関数の持つ本来の制約条件より「弱い制約条件」を実現したり、「不等式制約」を緩和した問題を実現するために用いられます。例えば、なるべく \(g\) は小さい方が良いが、必ずしも最小値 \(g = m\) でなくて良いという状況を作ることが可能です。

充足条件を不等式に変更した制約条件オブジェクトは次のようにして構築できます。

from amplify import BinarySymbolGenerator
from amplify.constraint import penalty

gen = BinarySymbolGenerator()
q = gen.array(3) # 長さ 3 の変数配列を生成
f = -q[0] - 0.5 * q[1] - 0.5 * q[2]         # 目的関数 f
g = q[0] * q[1] + q[1] * q[2] + q[2] * q[0] # ペナルティ関数 g (最小値 0)
p = penalty(g, le = 1)    # g <= 1 で充足判定を行う

>>> p
q_0 q_1 + q_0 q_2 + q_1 q_2 <= 1

不等式を表すキーワード引数は le (less equal), ge (greater equal), lt (less than), gt (greater than) が指定できます。上記の例では、ペナルティ関数 \(g \ge 0\) に対して、(\(g = 0\) ではなく) \(g \le 1\) で充足判定を行うことを意味します。

>>> result = solver.solve(f + p)
>>> for sol in result:
>>>     print(f"q = {q.decode(sol.values)}, f = {sol.energy}, g = {g.replace_all(sol.values)}")
q = [1, 1, 0], f = -1.5, g = 1.0
q = [1, 0, 1], f = -1.5, g = 1.0
q = [1, 0, 0], f = -1.0, g = 0.0

制約の充足判定を \(g \le 1\) とした解を複数得ることが出来ました。

制約条件クラスの演算子とメソッド

制約条件へのラベル付け

制約条件オブジェクトにはラベルを付けることができます。これにより、ラベルを元にした制約条件のフィルタリングを行うことができます。 具体的なフィルタリングの使用例については 制約条件のフィルタリング を参照してください。

以下のように、制約条件オブジェクトの生成関数の最終引数または label キーワード引数にラベルを与えます。

from amplify import BinarySymbolGenerator
from amplify.constraint import penalty, equal_to

gen = BinarySymbolGenerator()
q = gen.array(2)  # 長さ 2 の変数配列を生成
g1 = penalty(q, label="penalty_label")
g2 = equal_to(q, 2, label="equal_to_label")

>>> g1
penalty_label: q_0 + q_1 == 0
>>> g2
equal_to_label: q_0 + q_1 == 2

プリント時にラベルが制約式とともに表示される他、ラベルのみを確認することもできます。

>>> g1.label
'penalty_label'

制約条件を満たすか調べる

制約条件オブジェクトに変数値を与えることで制約条件を満たしているかを返します。変数値は、辞書・リスト・関数オブジェクトとして与えることが出来ますが、制約条件に含まれる変数のインデックス全てに対して解決可能でなければなりません。

from amplify import BinarySymbolGenerator
from amplify.constraint import one_hot

gen = BinarySymbolGenerator()
q = gen.array(4)  # 長さ 4 の変数配列を生成
f = one_hot(q)

>>> f.is_satisfied({0: 1, 1: 0, 2: 0, 3: 0})
True
>>> f.is_satisfied([0, 1, 0, 0])
True
>>> f.is_satisfied(lambda i: 1 if i == 0 else 0)
True

>>> f.is_satisfied({0: 1, 1: 0, 2: 0, 3: 1})
False
>>> f.is_satisfied([1, 1, 1, 1])
False
>>> f.is_satisfied(lambda i: 0)
False

ペナルティ関数の取得

制約条件オブジェクトの内部表現であるペナルティ関数を取得します。

from amplify import BinarySymbolGenerator, BinaryQuadraticModel
from amplify.constraint import one_hot, less_equal

gen = BinarySymbolGenerator()
q = gen.array(4)  # 長さ 4 の変数配列を生成
f1 = one_hot(q)
f2 = less_equal(q, 2)

>>> f1
q_0 + q_1 + q_2 + q_3 == 1
>>> f1.penalty
2 q_0 q_1 + 2 q_0 q_2 + 2 q_0 q_3 + 2 q_1 q_2 + 2 q_1 q_3 + 2 q_2 q_3 - q_0 - q_1 - q_2 - q_3 + 1

高次多項式を与えた場合や不等式制約など、ペナルティ関数の構築に補助変数が必要な場合はそのままでは取得できません。これは補助変数のインデックスが未確定なためです。

>>> f2
q_0 + q_1 + q_2 + q_3 <= 2
>>> f2.penalty
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
RuntimeError: Function to publish ancillary variables is not given

この場合は BinaryQuadraticModel 等に与え論理模型に変換した後で、ペナルティ関数が確認出来ます。ただし論理模型に変換されているため、入力変数のインデックスと論理変数のインデックスの間に変数変換が行われている可能性があることに注意が必要です。

>>> m = BinaryQuadraticModel(f2)
>>> m.input_constraints
[(q_0 + q_1 + q_2 + q_3 <= 2, 1)]
>>> m.input_constraints[0][0].penalty
2 q_0 q_1 + 2 q_0 q_2 + 2 q_0 q_3 - 2 q_0 q_4 - 2 q_0 q_5 + 2 q_1 q_2 + 2 q_1 q_3 - 2 q_1 q_4 - 2 q_1 q_5 + 2 q_2 q_3 - 2 q_2 q_4 - 2 q_2 q_5 - 2 q_3 q_4 - 2 q_3 q_5 + 2 q_4 q_5 + q_0 + q_1 + q_2 + q_3 + q_4 + q_5

制約条件の重みの設定

制約条件オブジェクトは、BinaryQuadraticModel 等の論理模型に与えることを前提としていますが、この際に制約条件の重みを設定することが出来ます。重みは制約条件オブジェクト単体では実質意味を持ちませんが、論理模型に変換されるときに、「重み」 \(\times\) 「ペナルティ関数」が目的関数に足され、イジングマシンへの入力となります。つまり、目的関数である多項式オブジェクトや行列クラス、またその他の制約条件オブジェクトに対する相対的な重みとして扱われます。

from amplify import BinarySymbolGenerator
from amplify.constraint import one_hot

gen = BinarySymbolGenerator()
q = gen.array(4)  # 長さ 4 の変数配列を生成
g = 1 * one_hot(q)

>>> g
(q_0 + q_1 + q_2 + q_3 == 1, 1)

制約条件オブジェクトに設定すると、制約条件と重みのペアで表現されます(BinaryConstraintTerm)。重みの値は後から変更することも可能ですが、制約条件自体は変更できません。

>>> g[1] = 2
>>> g
(q_0 + q_1 + q_2 + q_3 == 1, 2)

制約条件のコンテナ化

BinaryQuadraticModel 等の論理模型に与える時に、複数の制約条件オブジェクトを BinaryConstraints クラスとしてコンテナ化することが可能です。

from amplify import BinarySymbolGenerator
from amplify.constraint import one_hot

gen = BinarySymbolGenerator()
q = gen.array(2, 2)  # 2x2 の変数配列を生成
g1 = one_hot(q[0])
g2 = 2 * one_hot(q[1])
g = g1 + g2

>>> g
[(q_0 + q_1 == 1, 1), (q_2 + q_3 == 1, 2)]

重みが設定されている場合もそうで無い場合も、区別無く加算してコンテナ化することが出来ます。この時、重みが設定されていない制約条件オブジェクトはデフォルト値として 1 が与えられます。

下記の様にして制約条件オブジェクトコンテナ全体に重みを乗算・除算することが出来ます。

>>> g *= 2
>>> g
[(q_0 + q_1 == 1, 2), (q_2 + q_3 == 1, 4)]

参考

具体的な制約条件のコンテナ化や重みの設定例については One-hot 制約の定式化 を参照してください。

乗算・除算演算子

制約条件オブジェクトに対する乗算と除算は制約の重みの値を変更します。

乗算・除算対象と返却値の関係は次の通りです。

数値型

乗算・除算対象クラス

返却クラス

int

BinaryConstraint BinaryConstraintTerm

BinaryConstraintTerm

BinaryIntConstraint BinaryIntConstraintTerm

BinaryIntConstraintTerm

BinaryConstraints

BinaryConstraints

BinaryIntConstraints

BinaryIntConstraints

float

BinaryConstraint BinaryConstraintTerm

BinaryConstraintTerm

BinaryIntConstraint BinaryIntConstraintTerm

BinaryConstraints

BinaryConstraints

BinaryIntConstraints

加算演算子

制約条件オブジェクトに対する加算は、加算対象によって振る舞いを変えます。対象が多項式オブジェクトや行列オブジェクトの場合には論理模型を構築します。一方、対象が制約条件オブジェクトの場合には制約条件コンテナを構築します。

加算対象と返却値の関係は次の通りです。

制約条件クラス

加算対象クラス

返却クラス

BinaryConstraint BinaryConstraintTerm

BinaryConstraint BinaryConstraintTerm

BinaryConstraints

BinaryIntConstraint BinaryIntConstraintTerm

BinaryPoly

BinaryQuadraticModel

BinaryMatrix

BinaryIntPoly

BinaryIntMatrix

BinaryIntConstraint BinaryIntConstraintTerm

BinaryConstraint BinaryConstraintTerm

BinaryConstraints

BinaryIntConstraint BinaryIntConstraintTerm

BinaryIntConstraints

BinaryPoly

BinaryQuadraticModel

BinaryMatrix

BinaryIntPoly

BinaryIntQuadraticModel

BinaryIntMatrix