定式化ベンチマーク

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

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

-	Amplify	PyQUBO	dimod BQM (index)	dimod BQM (symbol)	dimod CQM*
記号演算	✅	✅	❌	✅	✅
目的関数	✅	✅	✅	✅	✅
制約条件	✅	✅**	❌***	❌***	✅*
高次多項式	✅	✅	❌	❌	❌
係数行列	✅	❌	❌	❌	❌
変数タイプ	B/S/I/R	B/S	B/S	B/S	B/S/I/R
対応マシン	さまざま	ユーザ次第	D-Wave のみ	D-Wave のみ	D-Wave のみ

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

*: QUBO 出力 (制約条件のペナルティ関数化) ができないためモデル作成のみ計測
**: ペナルティ関数を定義する必要あり
***: ペナルティ関数を目的関数に足すことで表現

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

ベンチマークコード

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,
    ]:
        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}")

ベンチマーク環境

CPU

Intel^R Core^TM i9-12900K
(E-Cores disabled)

OS

Ubuntu 22.04

Python 3.11

• amplify 1.0.0

• pyqubo 1.4.0

• dimod 0.12.14

ベンチマーク結果

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

定式化	モデル作成時間	QUBO 作成時間	合計時間
Amplify	1.309 ms	0.998 ms	2.307 ms 🏆
PyQUBO	266.8 ms	24.63 ms	291.5 ms (x126.3)
dimod BQM (index)	22.53 ms	41.58 ms	64.12 ms (x27.8)
dimod BQM (symbol)	1028 ms	58.48 ms	1087 ms (x471.2)
dimod CQM	895.9 ms	N/A ms	895.9 ms (x388.4)

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

定式化	モデル作成時間	QUBO 作成時間	合計時間
Amplify	0.063 s	0.077 s	0.140 s 🏆
PyQUBO	8.519 s	1.360 s	9.879 s (x70.7)
dimod BQM (index)	0.706 s	1.419 s	2.126 s (x15.2)
dimod BQM (symbol)	30.83 s	1.935 s	32.77 s (x234.6)
dimod CQM	26.93 s	N/A	26.93 s (x192.8)

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

定式化	モデル作成時間	QUBO 作成時間	合計時間
Amplify	2.992 s	2.946 s	5.938 s 🏆
PyQUBO	279.7 s	69.23 s	349.0 s (x58.8)
dimod BQM (index)	31.05 s	52.05 s	83.10 s (x14.0)
dimod BQM (symbol)	991.5 s	70.83 s	1062 s (x178.9)
dimod CQM	855.3 s	N/A	855.3 s (x144.0)