定式化ベンチマーク

Python のライブラリとして提供される数理最適化モデルについて、Amplify SDK との比較として定式化のベンチマークを取得しています。ここでは QUBO ソルバーの実行を前提として、巡回セールスマン問題の定式化 を例にモデルの作成し QUBO として出力するまでの実行時間を計測しました。

ただし、Amplify SDK を含む各ライブラリはそれぞれのカバーする機能が定式化の方法によっても異なります。ここでは次のように定式化を行い、それぞれの機能と方針についてまとめました。

-	Amplify	PyQUBO	dimod BQM (index)	dimod BQM (symbol)	dimod CQM1	PyQBPP double
記号演算	✅	✅	❌	✅	✅	✅
配列演算・einsum	✅	❌5	❌	❌	❌	✅
目的関数	✅	✅	✅	✅	✅	✅
制約条件	✅	✅2	❌3	❌3	✅1	✅
ペナルティ関数の自動生成	✅	❌	✅6	✅6	❌1	✅
高次多項式	✅	✅	❌	❌	❌	✅
係数行列	✅	❌	❌	❌	❌	❌
変数タイプ	B/S/I/R	B/S	B/S	B/S	B/S/I/R	B4
変数の自動符号化	✅	❌7	❌	❌	❌1	❌8
対応マシン	さまざま	ユーザ次第	D-Wave のみ	D-Wave のみ	D-Wave のみ	さまざま
ファイル入出力	LP/QPLIB	❌	❌	❌	LP9	❌
型ヒント	✅	❌	✅	✅	✅	❌

B: Binary, S: Ising Spin, I: Integer, R: Real

dimod の列は dimod 単体の機能を示します。D-Wave Ocean SDK 全体では、別のパッケージが一部の機能を提供します。

1. QUBO 出力 (制約条件のペナルティ関数化) ができないためモデル作成のみ計測

2. ペナルティ関数を定義する必要あり

3. ペナルティ関数を目的関数に足すことで表現

4. 宣言できる変数はバイナリのみ。整数は onehot_to_int、スピンは binary_to_spin で表現

5. Array は変数のコンテナであり、配列同士の要素ごとの演算やブロードキャストはできない

6. 線形の等式・不等式制約に限り add_linear_equality_constraint などでペナルティ項を生成できる

7. 整数は符号化の方法をクラスで明示的に選ぶ必要がある。実数変数は扱えない

8. 整数変数のみ。binarize による明示的な変換が必要

9. LP ファイルのみ

Amplify

import amplify


def tsp_for_amplify(ncity: int, distances: np.ndarray, dmax: float):
    q = amplify.VariableGenerator().array("Binary", ncity + 1, ncity)
    q[-1, :] = q[0, :]

    # 目的関数
    objective: amplify.Poly = amplify.einsum(
        "ij,ki,kj->", distances, q[:-1], q[1:]
    )

    # 制約条件
    constraints: amplify.ConstraintList = amplify.one_hot(
        q[:-1], axis=1
    ) + amplify.one_hot(q[:-1], axis=0)

    return objective + dmax * constraints


class BenchTspAmplify:
    def create_model(self, ncity: int, distances: np.ndarray, dmax: float):
        self.model = tsp_for_amplify(ncity, distances, dmax)

    def to_qubo(self):
        self.model.to_unconstrained_poly()

PyQUBO

import pyqubo


def tsp_for_pyqubo(ncity: int, distances: np.ndarray, dmax: float):
    # from https://github.com/recruit-communications/pyqubo/blob/master/notebooks/TSP.ipynb
    # NOTE: https://github.com/recruit-communications/pyqubo/blob/master/benchmark/benchmark.py
    # is not valid for TSP

    x = pyqubo.Array.create("c", (ncity, ncity), "BINARY")

    # Constraint not to visit more than two cities at the same time.
    time_const = 0.0
    for i in range(ncity):
        # If you wrap the hamiltonian by Const(...), this part is recognized as constraint
        time_const += pyqubo.Constraint(
            (sum(x[i, j] for j in range(ncity)) - 1) ** 2, label=f"time{i}"
        )

    # Constraint not to visit the same city more than twice.
    city_const = 0.0
    for j in range(ncity):
        city_const += pyqubo.Constraint(
            (sum(x[i, j] for i in range(ncity)) - 1) ** 2, label=f"city{j}"
        )

    # distance of route
    feed_dict = {}

    distance = 0.0
    for i in range(ncity):
        for j in range(ncity):
            for k in range(ncity):
                # we set the constant distance
                distance += distances[i, j] * x[k, i] * x[(k + 1) % ncity, j]

    # Construct hamiltonian
    A = pyqubo.Placeholder("A")
    H = distance + A * (time_const + city_const)

    feed_dict["A"] = dmax

    # Compile model
    return H.compile(), feed_dict


class BenchTspPyQubo:
    def create_model(self, ncity: int, distances: np.ndarray, dmax: float):
        self.model, self._feed_dict = tsp_for_pyqubo(ncity, distances, dmax)

    def to_qubo(self):
        self.model.to_qubo(index_label=False, feed_dict=self._feed_dict)

dimod BQM (index)

import dimod


def tsp_for_dimod_bqm(ncity: int, distances: np.ndarray, dmax: float):
    bqm = dimod.BinaryQuadraticModel(ncity * ncity, dimod.BINARY)

    # 目的関数
    for n in range(ncity):
        for i in range(ncity):
            for j in range(ncity):
                bqm.add_quadratic(
                    n * ncity + i,
                    ((n + 1) % ncity) * ncity + j,
                    distances[i, j],
                )

    # 行に対する制約
    for n in range(ncity):
        left = [(n * ncity + i, 1) for i in range(ncity)]
        bqm.add_linear_equality_constraint(left, dmax, -1)

    # 列に対する制約
    for i in range(ncity):
        left = [(n * ncity + i, 1) for n in range(ncity)]
        bqm.add_linear_equality_constraint(left, dmax, -1)

    return bqm


class BenchTspDimodBQM:
    def create_model(self, ncity: int, distances: np.ndarray, dmax: float):
        self.model = tsp_for_dimod_bqm(ncity, distances, dmax)

    def to_qubo(self):
        self.model.to_qubo()

dimod BQM (symbol math)

import dimod


def tsp_for_dimod_bqm_sym(
    ncity: int, distances: np.ndarray, dmax: float
) -> dimod.BinaryQuadraticModel:
    bqm = dimod.BinaryQuadraticModel(ncity * ncity, dimod.BINARY)
    vars = [
        [dimod.Binary(f"{n},{i}") for i in range(ncity)] for n in range(ncity)
    ]

    # 目的関数
    for n in range(ncity):
        for i in range(ncity):
            for j in range(ncity):
                bqm += distances[i, j] * vars[n][i] * vars[(n + 1) % ncity][j]

    # 行に対する制約
    for n in range(ncity):
        bqm += dmax * (sum(vars[n][i] for i in range(ncity)) - 1) ** 2

    # 列に対する制約
    for i in range(ncity):
        bqm += dmax * (sum(vars[n][i] for n in range(ncity)) - 1) ** 2

    return bqm  # type: ignore


class BenchTspDimodBQMSym:
    def create_model(self, ncity: int, distances: np.ndarray, dmax: float):
        self.model = tsp_for_dimod_bqm_sym(ncity, distances, dmax)

    def to_qubo(self):
        self.model.to_qubo()

dimod CQM

import dimod


def tsp_for_dimod_cqm(ncity: int, distances: np.ndarray, dmax: float):
    cqm = dimod.ConstrainedQuadraticModel()
    vars = [
        [dimod.Binary(f"{n},{i}") for i in range(ncity)] for n in range(ncity)
    ]

    # 目的関数
    obj = 0.0
    for n in range(ncity):
        for i in range(ncity):
            for j in range(ncity):
                obj += distances[i, j] * vars[n][i] * vars[(n + 1) % ncity][j]
    cqm.set_objective(obj)

    # 行に対する制約
    for n in range(ncity):
        cqm.add_constraint(sum(vars[n]) == 1)

    # 列に対する制約
    for i in range(ncity):
        cqm.add_constraint(sum(vars[n][i] for n in range(ncity)) == 1)

    return cqm


class BenchTspDimodCQM:
    def create_model(self, ncity: int, distances: np.ndarray, dmax: float):
        self.model = tsp_for_dimod_cqm(ncity, distances, dmax)

    def to_qubo(self):
        pass

PyQBPP

import pyqbpp.double as qbpp


def tsp_for_pyqbpp(ncity: int, distances: np.ndarray, dmax: float):
    variables = qbpp.var("x", shape=(ncity, ncity))

    # Constraints
    model = (
        qbpp.sum(qbpp.constrain(qbpp.vector_sum(variables, axis=1), equal=1))
        + qbpp.sum(qbpp.constrain(qbpp.vector_sum(variables, axis=0), equal=1))
    ) * dmax

    # Objective function
    distance = qbpp.array(distances.ravel().tolist(), shape=(ncity, ncity))
    successors = qbpp.concat([variables[1:], variables[:1]], axis=0)
    model += qbpp.einsum("jk,ij,ik->", distance, variables, successors)

    return model


class BenchTspPyQbpp:
    def create_model(self, ncity: int, distances: np.ndarray, dmax: float):
        self.model = tsp_for_pyqbpp(ncity, distances, dmax)

    def to_qubo(self):
        self.model.simplify_as_binary()

ベンチマークコード

import time


def make_distance(ncity: int) -> tuple[np.ndarray, float]:
    rng = np.random.default_rng(12345)
    x = rng.random(ncity)
    y = rng.random(ncity)
    distances = (
        (x[:, np.newaxis] - x[np.newaxis, :]) ** 2
        + (y[:, np.newaxis] - y[np.newaxis, :]) ** 2
    ) ** 0.5
    dmax: float = np.max(distances)  # type: ignore
    return distances, dmax


for ncity in [32, 100, 317]:
    distances, dmax = make_distance(ncity)

    for bench_class in [
        BenchTspAmplify,
        BenchTspPyQubo,
        BenchTspDimodBQM,
        BenchTspDimodBQMSym,
        BenchTspDimodCQM,
        BenchTspPyQbpp,
    ]:
        bench = bench_class()
        start = time.time()
        bench.create_model(ncity, distances, dmax)
        end = time.time()
        t1 = end - start
        start = time.time()
        bench.to_qubo()
        end = time.time()
        t2 = end - start
        print(f"{t1} {t2}")

ベンチマーク結果

PyQBPP (double) は 10,000 バイナリ変数を超える範囲を計測していません。

ベンチマーク環境

CPU

12th Gen Intel(R) Core(TM) i9-12900K
(E-Cores disabled)

OS

Linux-6.8.0-137-generic-x86_64-with-glibc2.43

Python 3.12.12

• amplify 1.7.0

• amplify 1.0.5

• pyqubo 1.5.0

• dimod 0.12.22

• pyqbpp 2026.8.16

サイズごとの定式化時間

各点は定式化の合計時間です。下にある点ほど高速です。

32 都市 (1,024 バイナリ変数)

定式化	モデル作成時間	QUBO 作成時間	合計時間
Amplify v1.7	0.24 ms	0.31 ms	0.55 ms 🏆
PyQUBO	142.60 ms	23.28 ms	165.88 ms (x301.3)
dimod BQM (index)	12.80 ms	29.18 ms	41.97 ms (x76.2)
dimod BQM (symbol)	611.76 ms	41.36 ms	653.12 ms (x1186.3)
dimod CQM	531.95 ms	N/A	531.95 ms (x966.2)
PyQBPP (double)	3.34 ms	4.94 ms	8.28 ms (x15.0)

100 都市 (10,000 バイナリ変数)

定式化	モデル作成時間	QUBO 作成時間	合計時間
Amplify v1.7	3.27 ms	5.35 ms	8.61 ms 🏆
PyQUBO	4.820 s	1.430 s	6.250 s (x725.6)
dimod BQM (index)	431.44 ms	1.038 s	1.469 s (x170.6)
dimod BQM (symbol)	18.425 s	1.447 s	19.872 s (x2307.0)
dimod CQM	15.433 s	N/A	15.433 s (x1791.7)
PyQBPP (double)	57.64 ms	61.36 ms	119.00 ms (x13.8)

317 都市 (100,489 バイナリ変数)

定式化	モデル作成時間	QUBO 作成時間	合計時間
Amplify v1.7	104.85 ms	180.43 ms	285.28 ms 🏆
PyQUBO	166.365 s	66.841 s	233.206 s (x817.5)
dimod BQM (index)	23.740 s	39.400 s	63.140 s (x221.3)
dimod BQM (symbol)	590.709 s	52.804 s	643.513 s (x2255.8)
dimod CQM	532.464 s	N/A	532.464 s (x1866.5)
PyQBPP (double)	N/A	N/A	N/A