OMMX AdapterでQUBOからサンプリングする#

ここでは巡回セールスマン問題を例として、問題をQUBOに変換しサンプリングを行う方法を説明します。

../_images/taraimawashi_businessman.png — たらい回しのイラスト（スーツ・男性）#

巡回セールスマン問題（TSP）は一人のセールスマンが複数の都市を順番に巡る方法を求める問題です。都市間の移動コストが与えられたときコストが最小になる経路を求めます。ここでは自己完結した例として、固定した乱数seedを使って10×10の領域に16都市を再現可能な形で生成します。

from random import Random

N = 16
rng = Random(42)
city_points = [
    (rng.uniform(0.0, 10.0), rng.uniform(0.0, 10.0))
    for _ in range(N)
]

都市の位置をプロットしてみましょう

%matplotlib inline
from matplotlib import pyplot as plt

x_coords, y_coords = zip(*city_points)
plt.scatter(x_coords, y_coords)
plt.xlabel('X Coordinate')
plt.ylabel('Y Coordinate')
plt.title('ランダム生成した都市配置')
plt.show()

コストとして単純に移動距離を考えましょう。\(i\)番目の都市と\(j\)番目の都市の距離 \(d(i, j)\)を計算しておきます。

def distance(x, y):
    return ((x[0] - y[0])**2 + (x[1] - y[1])**2)**0.5

# 各都市間の距離
d = [[distance(city_points[i], city_points[j]) for i in range(N)] for j in range(N)]

これを使って次のような最適化問題としてTSPを定式化します。まずある時刻 \(t\) に都市 \(i\) にいるかどうかをバイナリ変数 \(x_{t, i}\) で表します。このとき、以下の制約を満たすような \(x_{t, i}\) を求めます。するとセールスマンが移動する距離は次で与えられます：

\[ \sum_{t=0}^{N-1} \sum_{i, j = 0}^{N-1} d(i, j) x_{t, i} x_{(t+1 \% N), j} \]

ただし \(x_{t, i}\) は自由に取れるわけではなく、各時刻 \(t\) において一箇所の都市にしかいられないという制約と各都市について一度だけ訪れるという制約

\[ \sum_{i=0}^{N-1} x_{t, i} = 1, \quad \sum_{t=0}^{N-1} x_{t, i} = 1 \]

を満たす必要があります。これらを合わせてTSPは制約付き最適化問題として定式化できます

\[ \begin{aligned} \min \quad & \sum_{t=0}^{N-1} \sum_{i, j = 0}^{N-1} d(i, j) x_{t, i} x_{(t+1 \% N), j} \\ \text{s.t.} \quad & \sum_{i=0}^{N-1} x_{t, i} = 1 \quad (\forall t = 0, \ldots, N-1) \\ \quad & \sum_{t=0}^{N-1} x_{t, i} = 1 \quad (\forall i = 0, \ldots, N-1) \end{aligned} \]

これに対応する ommx.Instance は次のように作成できます

from ommx import DecisionVariable, Instance

x = [[
        DecisionVariable.binary(
            i + N * t,  # 決定変数のID
            name="x",           # 決定変数の名前、解を取り出すときに使う
            subscripts=[t, i])  # 決定変数の添字、解を取り出すときに使う
        for i in range(N)
    ]
    for t in range(N)
]

objective = sum(
    d[i][j] * x[t][i] * x[(t+1) % N][j]
    for i in range(N)
    for j in range(N)
    for t in range(N)
)
place_constraint = {
    t: (sum(x[t][i] for i in range(N)) == 1)
        .set_name("place")
        .add_subscripts([t])
    for t in range(N)
}
time_constraint = {
    i + N: (sum(x[t][i] for t in range(N)) == 1)
        .set_name("time")
        .add_subscripts([i])
    for i in range(N)
}

instance = Instance.from_components(
    decision_variables=[x[t][i] for i in range(N) for t in range(N)],
    objective=objective,
    constraints={**place_constraint, **time_constraint},
    sense=Instance.MINIMIZE
)

バイナリの決定変数の作成時 DecisionVariable.binary に追加した決定変数の名前と添字は後で得られたサンプルを解釈する際に使います。

OpenJijによるサンプリング#

ommx-openjij-adapter のinput classに属するのは、任意次数の多項式目的関数を持つバイナリ変数のみの制約なし最小化問題です。上で作成したTSPインスタンスには制約があるため、有限のペナルティ重みを指定して明示的に準備します。その後、prepared.input の Instance をAdapterへ渡し、得られたsampleを変換元モデルに対して明示的に評価します。

from ommx_openjij_adapter import (
    OMMXOpenJijSAAdapter,
    OpenJijPreparationConfig,
)

config = OpenJijPreparationConfig(
    uniform_penalty_weight=20.0,
)
prepared = OMMXOpenJijSAAdapter.prepare(instance, config=config)

prepared_samples = OMMXOpenJijSAAdapter.sample(
    prepared.input,
    num_reads=16,
)
sample_set = prepared.evaluate_source(prepared_samples)
sample_set.summary

sample() は SampleSet を返します。これは決定変数のサンプル値に加えて、評価した目的関数値と制約違反を保持します。SampleSet.summary はこの情報の要約を表示します。prepared.evaluate_source() が準備済み入力のsampleを変換元モデルに対して評価するため、その feasible 列は変換元の制約付き問題に対する実行可能性を示します。

config のペナルティ重みはOpenJij backend samplerのパラメータではなく、明示的な準備に対する指定です。有限ペナルティは実行可能なサンプルを得やすくしますが、すべてのサンプルが変換元の問題に対して実行可能になることを保証しません。

準備内容の確認#

check_preparation はインスタンスを変更せずに、変換元モデルと準備configを検査します。prepare は検査した変換を実行し、監査用レポートを prepared.report に保存します。

report = prepared.report
config_used = report.config
final = report.input_applicability
outcomes = {
    "source_membership": report.source_check.source_membership.is_member,
    "steps": [step.operation for step in report.steps],
    "preparation_failures": report.preparation_failures,
    "input_applicability": None if final is None else final.is_applicable,
}
config_used, outcomes

report.config は、正規化済みで実際に使われた不変のpreparation設定を記録します。これは設定の監査記録です。その他のfieldは、次のいずれか1つの終端状態を表します。

状態	`source_check`	`preparation_failures`	`input_applicability`
Source rejected	preparation source classの外	空	`None`
Phase rejected	accepted	ownerの `operation` を含む非空の値	`None`
Candidate rejected	accepted	空	non-applicable report
Success	accepted	空	applicable report

source_check は構造的なsource-class membershipです。operation availabilityと preparation policyは、そのoperationを実行するphaseが所有し、 preparation_failures に記録します。steps は終端状態までに完了したoperationの prefixです。独立したoutcomeや、合成された数学的guaranteeではありません。共通のpreparation policy、guarantee、自動選択は OMMX issue #1111 で扱います。このprototype が既定で使うのは、利用可能な厳密operationだけです。離散的なinteger slack近似には OpenJijPreparationConfig で allow_approximate_integer_slack=True を設定する必要があり、inequality_integer_slack_max_range の指定だけでは近似への同意になりません。同じConfigで uniform_penalty_weight または penalty_weights を明示的に設定すると finite-penalty preparationが選択されますが、制約付き入力をAdapterが直接または厳密にサポートするという意味ではありません。

変数boundから不等式が実行不可能だと証明できた場合、check_preparation と prepare はcore所有の InfeasibleDetected を送出します。従来の InfeasibleDetected importは同じexception objectへのalias として残ります。これはモデル自体の性質であり、Adapter applicabilityの失敗ではありません。

使用されるInteger変数ごとに最大53個の補助bitという条件を検査しますが、これは OMMXのInteger-to-Binary log encodingのavailability limitです。OpenJij Adapterのinput classの性質ではなく、serialized semanticsをforward compatibilityのためにreaderが安全に解釈できるかを管理する ommx.v2.Feature とも別物です。OpenJijが直接受け付ける Spin入力を含むSpin変数のサポートは OMMX issue #1082 で別途管理しています。

各制約条件毎のfeasibilityを見るには summary_with_constraints プロパティを使います。

sample_set.summary_with_constraints

より詳しい情報は SampleSet.decision_variables_df() 及び SampleSet.constraints_df() メソッドを使って取得できます。

sample_set.decision_variables_df().head(2)

sample_set.constraints_df().head(2)

得られたサンプルを取得するには SampleSet.extract_decision_variables メソッドを使います。これは ommx.DecisionVariables を作る時に登録した name と subscripts を使ってサンプルを解釈します。例えば sample_id=1 の x という名前の決定変数の値を取得するには次のようにすると dict[subscripts, value] の形で取得できます。

sample_id = 1
x = sample_set.extract_decision_variables("x", sample_id)
t = 2
i = 3
x[(t, i)]

\(x_{t, i}\)に対するサンプルが得れたのでこれをTSPのパスに変換します。これは今回の定式化自体に依存するので自分で処理を書く必要があります。

def sample_to_path(sample: dict[tuple[int, ...], float]) -> list[int]:
    path = []
    for t in range(N):
        for i in range(N):
            if sample[(t, i)] == 1:
                path.append(i)
    return path

これを表示してみましょう。まず元の問題に対してfeasibleであるサンプルのIDを取得します。

feasible_ids = sample_set.summary.query("feasible == True").index
feasible_ids

これらについて最適化された経路を表示しましょう

fig, axie = plt.subplots(3, 3, figsize=(12, 12))

for i, ax in enumerate(axie.flatten()):
    if i >= len(feasible_ids):
        break
    s = feasible_ids[i]
    x = sample_set.extract_decision_variables("x", s)
    path = sample_to_path(x)
    xs = [city_points[i][0] for i in path] + [city_points[path[0]][0]]
    ys = [city_points[i][1] for i in path] + [city_points[path[0]][1]]
    ax.plot(xs, ys, marker='o')
    ax.set_title(f"Sample {s}, objective={sample_set.objectives[s]:.2f}")

plt.tight_layout()
plt.show()

OMMX AdapterでQUBOからサンプリングする

目次

OMMX AdapterでQUBOからサンプリングする#

OpenJijによるサンプリング#

準備内容の確認#