Zenn
広島大学コンピューティングラボPublicationへの投稿
🕸️

QUBO式生成ツールの処理性能比較：QAPを用いたベンチマーク評価

Koji Nakano
に公開
C++
pyhon
数理最適化
量子アニーリング
qubo
tech

Quadratic Assignment Problem (QAP)を用いた性能比較

組合せ最適化問題を解く手法として，等価なQUBO（Quadratic Unconstrained Binary Optimization）式を生成し，それをQUBOソルバーで解くアプローチが有望視されており，これに対応するさまざまな定式化ツールが開発されています．

QUBOソルバーの性能比較は数多く行われていますが，ソルバーに入力するQUBO式そのものを生成する処理の性能比較は，これまでほとんど注目されてきませんでした．しかしながら，QUBO式の生成に過度な時間を要するようでは，ソルバーの性能を評価する以前に，実用上のボトルネックとなりかねません．

そこで本記事では，Quadratic Assignment Problem（QAP）を題材に，各ツールによるQUBO式の生成性能をベンチマークします．具体的には，以下の10種類のツールを用いてQUBO式の生成プログラムを実装し，それぞれの処理時間を比較します．なお，ソルバー機能を内蔵しているツールも含まれますが，本記事ではQUBO式の生成処理のみに注目し，ソルバー性能の比較は行いません．

ツール 開発・提供元 ベース言語 バックエンド 備考
QUBO++ 広島大学 C++ -
PyQBPP 広島大学 Python - QUBO++のPython版
Amplify Fixstars Python C++
JijModeling Jij Python Rust
PyQUBO リクルートコミュニケーションズ Python C++
dimod D-Wave Systems Python CPython
TYTAN Blueqat Python SymEngine
qubovert Joseph T. Iosue Python -
SymPy オープンソース開発者コミュニティ Python 汎用記号処理ツール
SymEngine オープンソース開発者コミュニティ C++ - SymPyのC++実装，Pythonラッパーあり

ここで言う「QUBO式の生成」とは，例えば以下のように式を展開・整理し，QUBO形式に適した形へ変換する記号処理のことを指します．

\begin{align*} (a+b-1)(a+b-1) &= 1 +a^2 +ba +ab +b^2 -a -b -a -b\\ &= 1 -a -b +2ab \end{align*}

この変形では，a^2 = a，b^2 = b といった二値変数の性質も用いており，QUBOの各項（定数項・1次項・2次項）の係数を抽出することが可能になります．
一方で，QUBO式を直接辞書形式で与えることも可能です．たとえば，次のように係数を直接指定する方法です：

{"1":1, "a":-1, "b":-1, "a*b":2}

しかし，このような形式は極めて単純な問題にしか適用できず，複雑な問題においては手作業で係数を導出することは現実的ではありません．開発効率や保守性を考慮すると，記号処理によるQUBO式の自動生成は実用上不可欠と言えます．

なお，各ツールについては，提供されている組み込み関数やプレースホルダ機能を活用し，可能な範囲で効率化を図りましたが，すべての機能を完全に把握しているわけではありません．また，マニュアルやドキュメントが不十分なケースも多く，各ツールにおいて最も効率的な記述方法になっていない可能性がある点については，あらかじめご承知おきください．

Quadratic Assignment Problem (QAP)とは

サイズsのQAPとは，s個の施設をs個の場所に1対1に対応付けて配置する問題です．施設間の物流量と場所間の距離がそれぞれ与えられており，配置後の物流量と距離の積の総和を最小化することを目的とする組合せ最適化問題です．

施設と場所を，それぞれ0からs-1までの整数で表します．施設iとjの間の物流量を\mathrm{flow}_{i,j}，場所kとlの間の距離を\mathrm{dist}_{k,l}と表します．1対1関数\sigma: \{0, 1, \dotsc, s-1\} \to \{0, 1, \dotsc, s-1\}を用いて，施設と場所の割当を表します．各施設iは場所\sigma(i)に割り当てられます．このときの目的関数は次の式で表されます．

\sum_{i=0}^{s-1}\sum_{j=0}^{s-1}\mathrm{flow}_{i,j}\mathrm{dist}_{\sigma(i),\sigma(j)}

この値を最小化するような1対1関数\sigmaを求める問題がQAPです．

下の図は，s=5の場合を示しています．

QAPをQUBO式として表す

s\times sの二値変数の２次元配列xを用いることで，1対1関数を表現することができます．各行と各列に1が一つだけ存在するという制約(one-hot制約と呼ばれる)を設定します．これは巡回セールスマン問題において都市の順列を表現する方法と同じです．ここで，上の図のように，x_{i,j}=1のとき，施設iが場所jに割り当てられている，すなわち\sigma(i)=jであるとみなします．この制約条件は，QUBO式として次のように表すことができます．

\begin{align*} \mathrm{constraint} &= \sum_{i=0}^{s-1}(1-\sum_{j=0}^{s-1} x_{i,j})^2+ \sum_{j=0}^{s-1}(1-\sum_{i=0}^{s-1} x_{i,j})^2 \end{align*}

この式は各行各列の合計が1であるとき，最小値0を取ります．ツールによっては，等号を使った制約条件の形式で次のように記述することもできます．

\begin{align*} \mathrm{constraint} &= \sum_{i=0}^{s-1}\left(\sum_{j=0}^{s-1} x_{i,j}=1\right)+ \sum_{j=0}^{s-1}\left(\sum_{i=0}^{s-1} x_{i,j}=1\right) \end{align*}

また，目的関数はx_{i,j}を用いて次のように表すことができます．

\begin{align*} \mathrm{objective} &=\sum_{i=0}^{s-1}\sum_{j=0}^{s-1}\sum_{k=0}^{s-1}\sum_{l=0}^{s-1} x_{i,k}x_{j,l}\mathrm{flow}_{i,j}\mathrm{dist}_{k,l} \end{align*}

ここで，x_{i,k}=x_{j,l}=1のとき，\sigma(i)=kかつ\sigma(j)=lとなるため，このQUBO式の値はQAPの目的関数の値と一致することがわかります．この制約条件と目的関数を，制約条件を優先するように重み付けして足し合わせれば，QAPに対応するQUBO式が得られます．制約条件を十分に優先するために，係数\mathrm{priority}は\mathrm{objective}の値を考慮して十分大きい値に設定します．

\begin{align*} \mathrm{QUBO} &= \mathrm{objective} + \mathrm{priority}\times \mathrm{constraint} \end{align*}

ベンチマークに用いるQAP

ランダムに生成した QAP（Quadratic Assignment Problem）を用いて，QUBO式の生成に要する時間（モデル構築時間）を評価します．距離行列 \mathrm{dist}_{i,j}\ (i \ne j) には，1から100までのランダムな整数を割り当てます．物流行列 \mathrm{flow}_{i,j} は，各施設 i について，異なる d 個の施設 j\ (j \ne i) に対してのみ，1から100までのランダムな整数を割り当て，それ以外は 0 とします．すなわち，物流量はスパースであり，各施設から出る物流はちょうど d 件となります．物流をグラフとみなすと，全ノードの出次数(degree)がちょうど d の有向グラフとなります．本評価では，d = 5 として処理時間計測を行いました．

QAPをQUBOに変換するプログラム

QUBO式の生成に要する時間を評価するため，固定した疑似乱数を用いて，1以上100以下のランダムな値を \mathrm{dist}_{i,j} および \mathrm{flow}_{k,l} として2次元配列に保存します．このとき，\mathrm{flow}_{k,l} には，指定された d に基づいてスパースな物流量を設定します．また，C++とPythonで同一のQAPを生成できるように，組み込みの疑似乱数器は使用せず，両者で共通の生成法に基づくオリジナルの乱数生成器を用います．

QUBO式の生成は，以下の手順で行います．

  • 二値変数の2次元配列 x の生成
  • 制約条件に対応するQUBO式 constraint の生成
  • 目的関数に対応するQUBO式 objective の生成
  • 最終的な整理されたQUBO式 qubo の構築
    これらに要する処理時間を計測します．なお，言語仕様や都合や効率化のため，constraint や objective を明示的に作成せず，quboを直接構築する場合もあります．なお，乱数によって \mathrm{dist} や \mathrm{flow} を作成する処理時間は含めません．ここで「整理済み」とは，同一項は統合され，定数項・一次項・二次項の係数がなんらかの順序で順次取り出せる形式のデータ構造が得られている状態を指します．例えば，(a + b - 1)^2 を展開し整理した結果として，以下のように統合・整理された項が取得できる状況を想定しています．
  • 定数項：1
  • 一次項：-a，-b
  • 二次項：2ab
    ツールによって，これらの項は辞書やリストなどの形式で保持されます．いずれの場合も，これらの項を順次取り出せるデータ構造が得られた時点で，QUBO式の生成が完了したものと見なします．

各プログラムは，オプション-oが指定された場合には，生成されたQUBOの各項とその係数を辞書式順序で出力し，異なる実装間で同一の結果が得られることを確認できるようにします．

補足
以下のプログラムでは，辞書式順序として「変数のインデックス順」と「変数の文字列表現に基づく順」の両方が混在していますが，インデックスが1桁の整数である場合，両者は一致します．そのため，そのツールの都合のよい方法を採用しています．

SymPyによるプログラム

SymPyを用いて，QAPをQUBOに変換するプログラムです．

import time
import argparse
import sympy as sp

class sp_Binary(sp.Symbol):
    def _eval_power(self, _):
        return self

def fixed_random(n: int) -> int:
    if not hasattr(fixed_random, "seed"):
        fixed_random.seed = 123456789
    fixed_random.seed = (fixed_random.seed * 1664525 + 1013904223) & 0xFFFFFFFF
    return fixed_random.seed % n

def select(size: int, index: int, degree: int):
    work = [i for i in range(size) if i != index]
    result = []
    for _ in range(degree):
        k = fixed_random(len(work))
        result.append(work[k])
        work[k], work[-1] = work[-1], work[k]
        work.pop()
    return result

parser = argparse.ArgumentParser(description="Generate QUBO of a sample QAP.")
parser.add_argument("-s", "--size", type=int,
                    required=True, help="Size of the QAP.")
parser.add_argument("-d", "--degree", type=int,
                    required=True, help="Degree of the QAP.")
parser.add_argument("-o", "--output", action="store_true",
                    help="Output the obtained QUBO.")
args = parser.parse_args()
size = args.size
degree = args.degree
output = args.output

dist = [[fixed_random(100) + 1 if i != j else 0 for j in range(size)]
        for i in range(size)]
flow = [(i, j, fixed_random(100) + 1)
        for i in range(size) for j in select(size, i, degree)]

start_time = time.perf_counter()

x = [[sp_Binary(f"x[{i}][{j}]") for j in range(size)] for i in range(size)]
constraint = sum((sum(x[i]) - 1) ** 2 for i in range(size)) + \
    sum((sum(x[i][j] for i in range(size)) - 1) ** 2 for j in range(size))
objective = sum(
    f * x[i][k] * x[j][l] * dist[k][l]
    for i, j, f in flow
    for k in range(size)
    for l in range(size)
)
qubo = (objective + 10000 * constraint).expand().as_coefficients_dict()

elapsed_time = time.perf_counter() - start_time

constant = 0
linear = {}
quadratic = {}
for term, coeff in qubo.items():
    if term == 1:
        constant = coeff
    elif term.is_Symbol:
        linear[term] = coeff
    elif term.is_Mul and sum(arg.is_Symbol for arg in term.args) == 2:
        quadratic[term] = coeff
if output:
    print(constant)
    for key, value in sorted(linear.items(), key=str):
        print(f"{key} {value}")
    for key, value in sorted(quadratic.items(), key=str):
        print(f"{key} {value}")

print(f"{size},{degree},{len(linear)},{len(quadratic)},{elapsed_time}")

このSymPyプログラムの構成を簡単に説明します．

  • sp_Binary：二値変数クラスを定義しています．二値変数は何乗しても元の値と変わらないという性質（x^2 = x）を利用しています．
  • fixed_random(n: int)：固定された乱数列を生成します．
  • select(size: int, index: int, degree: int)：各 i(index) に対して，異なる j を d(degree) 個選びます．
  • dist : 二次元配列（リストのリスト）であり，ランダムな距離値を i \ne j の要素に割り当てています．
  • flow：スパースな物流量行列をタプル (i, j, f) (\mathrm{flow}_{i,j}=f)のリストとして格納しています．
  • sp_Binary のリストのリストを用いて，二値変数の2 次元配列 x を定義しています．
  • constraint：各行・各列の合計が 1 になる制約式です．
  • objective：QAP の目的関数に対応する式です．
  • qubo：制約と目的関数を合算した後，expand と as_coefficients_dict を用いて QUBO 式を展開し，係数付き項の辞書として取得しています．

このプログラムを qap_sympy.py として保存し，次のコマンドを実行します．

$ python3 qap_sympy.py -s 3 -d 1 -o

すると次の出力結果が得られます．

60000
x[0][0] -20000
x[0][1] -20000
x[0][2] -20000
x[1][0] -20000
x[1][1] -20000
x[1][2] -20000
x[2][0] -20000
x[2][1] -20000
x[2][2] -20000
x[0][0]*x[0][1] 20000
x[0][0]*x[0][2] 20000
x[0][0]*x[1][0] 20000
x[0][0]*x[1][1] 2338
x[0][0]*x[1][2] 6772
x[0][0]*x[2][0] 20000
x[0][0]*x[2][1] 2170
x[0][0]*x[2][2] 6510
x[0][1]*x[0][2] 20000
x[0][1]*x[1][0] 2398
x[0][1]*x[1][1] 20000
x[0][1]*x[1][2] 3052
x[0][1]*x[2][0] 2310
x[0][1]*x[2][1] 20000
x[0][1]*x[2][2] 1680
x[0][2]*x[1][0] 6622
x[0][2]*x[1][1] 4792
x[0][2]*x[1][2] 20000
x[0][2]*x[2][0] 6160
x[0][2]*x[2][1] 5740
x[0][2]*x[2][2] 20000
x[1][0]*x[1][1] 20000
x[1][0]*x[1][2] 20000
x[1][0]*x[2][0] 20000
x[1][1]*x[1][2] 20000
x[1][1]*x[2][1] 20000
x[1][2]*x[2][2] 20000
x[2][0]*x[2][1] 20000
x[2][0]*x[2][2] 20000
x[2][1]*x[2][2] 20000
3,1,9,9,30,0.00747319

定数項，1次項，2次項が係数とともに辞書式順序で出力されます．最後の行には，問題サイズや生成時間などの統計情報が記録されます．
以下のすべてのプログラムにおいて，同一の定数項，1次項，2次項が得られます．

SymEngineによるプログラム

SymEngineを用いて，QAPをQUBOに変換するプログラムです．

import time
import argparse
from collections import defaultdict
import symengine as se

[省略]

start_time = time.perf_counter()

x = [[se.Symbol(f"x[{i}][{j}]") for j in range(size)] for i in range(size)]
constraint = sum((sum(x[i]) - 1) ** 2 for i in range(size)) + \
    sum((sum(x[i][j] for i in range(size)) - 1) ** 2 for j in range(size))
objective = sum(
    f * x[i][k] * x[j][l] * dist[k][l]
    for i, j, f in flow
    for k in range(size)
    for l in range(size)
)
qubo = se.expand(objective + 10000 * constraint).as_coefficients_dict()

constant = 0
linear = defaultdict(int)
quadratic = defaultdict(int)
for term, coeff in qubo.items():
    if term == 1:
        constant = coeff
    elif term.is_Mul and sum(arg.is_Symbol for arg in term.args) == 2:
        args = term.args
        if args[0].name < args[1].name:
            quadratic[(args[0], args[1])] += coeff
        else:
            quadratic[(args[1], args[0])] += coeff
    elif term.is_Symbol:
        linear[term] += coeff
    elif term.is_Pow and term.exp == 2 and term.base.is_Symbol:
        linear[term.base] += coeff

elapsed_time = time.perf_counter() - start_time

if output:
    print(constant)
    for key, value in sorted(linear.items(), key=str):
        print(f"{key} {value}")
    for key, value in sorted(quadratic.items(), key=str):
        print(f"{key[0]}*{key[1]} {value}")

print(f"{size},{degree},{len(linear)},{len(quadratic)},{elapsed_time}")

SymEngine は SymPy とある程度の互換性を持ち，かつ C++ で実装されているため，SymPy よりも高速に動作します．ただし，_eval_power などの一部のカスタマイズ機能は提供されていないため，式は単に記号として展開し，辞書形式でquboに保存しています．この qubo では，同一項の統合が行われていません．

そこで，後処理として，quboの辞書を読み出し，定数項を constant，一次項を linear，二次項を quadratic に分類・統合しています．そこで，公平な評価のため，この後処理も含めた全体の処理時間を測定対象としています．

qubovertによるプログラム

qubovertを用いて，QAPをQUBOに変換するプログラムです．

import time
import argparse
from collections import defaultdict
import qubovert as qv

[省略]

start_time = time.perf_counter()

x = [[qv.boolean_var(f"x[{i}][{j}]") for j in range(size)] for i in range(size)]

constraint = sum((sum(x[i]) - 1) ** 2 for i in range(size)) + \
    sum((sum(x[i][j] for i in range(size)) - 1) ** 2 for j in range(size))
objective = sum(
    f * x[i][k] * x[j][l] * dist[k][l]
    for i, j, f in flow
    for k in range(size)
    for l in range(size)
)

qubo = objective + 10000 * constraint

elapsed_time = time.perf_counter() - start_time

constant = 0
linear = {}
quadratic = {}
for term, coeff in qubo.items():
    if len(term) == 0:
        constant = coeff
    elif len(term) == 1:
        linear[term[0]] = coeff
    elif len(term) == 2:
        quadratic[term] = coeff

if output:
    print(constant)
    for key, value in sorted(linear.items(), key=str):
        print(f"{key} {value}")
    for key, value in sorted(quadratic.items(), key=str):
        print(f"{key[0]}*{key[1]} {value}")

print(f"{size},{degree},{len(linear)},{len(quadratic)},{elapsed_time}")

二値変数の定義にqv.boolean_varを用いている点以外は，SymPyのプログラムとほぼ同じです．

TYTANによるプログラム

TYTANを用いて，QAPをQUBOに変換するプログラムです．

import time
import argparse
import tytan

[省略]

start_time = time.perf_counter()

x = tytan.symbols_list([size, size], 'x{}_{}')
constraint = sum((sum(x[i]) - 1) ** 2 for i in range(size)) + sum(
    (sum(x[:, j]) - 1) ** 2 for j in range(size)
)
objective = sum(
    f * x[i, k] * x[j, l] * dist[k][l]
    for i, j, f in flow
    for k in range(size)
    for l in range(size)
)
qubo, offset = tytan.Compile(objective + 10000 * constraint).get_qubo()

elapsed_time = time.perf_counter() - start_time

array = qubo[0]
linear = {}
quadratic = {}
for i in range(size*size):
    if array[i, i] != 0:
        linear[i] = int(array[i, i])
    for j in range(i+1, size*size):
        if array[i, j] != 0:
            quadratic[(i, j)] = int(array[i, j])
if output:
    print(int(offset))
    for key, value in sorted(linear.items()):
        print(f"x[{key//size}][{key%size}] {value}")
    for key, value in sorted(quadratic.items()):
        print(
            f"x[{key[0]//size}][{key[0]%size}]*x[{key[1]//size}][{key[1]%size}] {value}")

print(f"{size},{degree},{len(linear)},{len(quadratic)},{elapsed_time}")

この Tytan プログラムは SymPy プログラムとほぼ同一であり，異なる部分のみを以下に説明します．

  • symbols_list を用いて，二値変数の 2 次元配列 x を定義しています．
  • スライス x[i, :] や x[:, j] を用いることで， constraint を簡潔に記述しています．
  • Compile と get_qubo 関数を用いて，QUBO 式 quboと定数項offset を求めています．

dimodによるプログラム

dimodを用いて，QAPをQUBOに変換するプログラムです．

import time
import argparse
import dimod

[省略]

start_time = time.perf_counter()

x = [[dimod.Binary(f"x[{i}][{j}]") for j in range(size)] for i in range(size)]
constraint = sum((sum(x[i]) - 1) ** 2 for i in range(size)) + \
    sum((sum(x[i][j] for i in range(size)) - 1) ** 2 for j in range(size))
objective = sum(
    f * x[i][k] * x[j][l] * dist[k][l]
    for i, j, f in flow
    for k in range(size)
    for l in range(size)
)
qubo = objective + 10000 * constraint

elapsed_time = time.perf_counter() - start_time

constant = int(qubo.offset)
linear = {}
quadratic = {}
for var, coeff in qubo.linear.items():
    linear[var] = int(coeff)
for vars, coeff in qubo.quadratic.items():
    if (vars[0] < vars[1]):
        quadratic[(vars[0], vars[1])] = int(coeff)
    else:
        quadratic[(vars[1], vars[0])] = int(coeff)

if output:
    print(constant)
    for key, value in sorted(linear.items(), key=str):
        print(f"{key} {value}")
    for key, value in sorted(quadratic.items(), key=str):
        print(f"{key[0]}*{key[1]} {value}")

print(f"{size},{degree},{len(linear)},{len(quadratic)},{elapsed_time}")

二値変数 x の定義に dimod.Binary を用いている点を除けば，この dimod プログラムは SymPy 版とほぼ同一です．

PyQUBOによるプログラム

PyQUBOを用いて，QAPをQUBOに変換するプログラムです．

import time
import argparse
import pyqubo as pq

[省略]

x = pq.Array.create("x", shape=(size, size), vartype="BINARY")
constraint = sum((sum(x[i, :]) - 1) ** 2 for i in range(size)) + sum(
    (sum(x[:, j]) - 1) ** 2 for j in range(size)
)
objective = sum(
    f * x[i, k] * x[j, l] * dist[k][l]
    for i, j, f in flow
    for k in range(size)
    for l in range(size)
)
model = objective + 10000 * constraint
qubo, offset = model.compile().to_qubo()

elapsed_time = time.perf_counter() - start_time

linear = {}
quadratic = {}
for key, value in qubo.items():
    if (key[0] == key[1]):
        linear[key] = int(value)
    else:
        if key[0] > key[1]:
            key = (key[1], key[0])
        quadratic[key] = int(value)

if output:
    print(int(offset))
    for key, value in sorted(linear.items()):
        print(key[0], value)
    for key, value in sorted(quadratic.items()):
        print(f"{key[0]}*{key[1]}", value)


print(f"{size},{degree},{len(linear)},{len(quadratic)},{elapsed_time}")

この PyQUBO プログラムも SymPy プログラムとほぼ同一であり，異なる部分のみを以下に説明します．

  • Array.create を用いて，二値変数の 2 次元配列 x を定義しています．
  • constraintはTytanと同一です．
  • compile と to_qubo 関数を用いて，QUBO 式 qubo と定数項 offset を求めています．

Amplifyによるプログラム

Amplifyを用いて，QAPをQUBOに変換するプログラムです．

import argparse
import time
import numpy as np
import amplify

[省略]

start_time = time.perf_counter()

gen = amplify.VariableGenerator()
x = gen.array("BINARY", (size, size), name="x")

constraint = amplify.one_hot(x, axis=0) + amplify.one_hot(x, axis=1)
objective = amplify.sum(f * amplify.einsum("k,l,kl->",
                        x[i], x[j], dist) for i, j, f in flow)

qubo = (objective + 10000 * constraint).to_unconstrained_poly().as_dict()

elapsed_time = time.perf_counter() - start_time

constant = 0
linear = {}
quadratic = {}
for key, value in sorted(qubo.items()):
    if len(key) == 0:
        constant = int(value)
    elif len(key) == 1:
        linear[key] = int(value)
    elif len(key) == 2:
        quadratic[key] = int(value)

if output:
    print(f"{constant}")
    for key, value in sorted(linear.items()):
        i, k = divmod(key[0], size)
        print(f"x[{i}][{k}] {value}")
    for key, value in sorted(quadratic.items()):
        i1, k1 = divmod(key[0], size)
        i2, k2 = divmod(key[1], size)
        print(f"x[{i1}][{k1}]*x[{i2}][{k2}] {value}")

print(f"{size},{degree},{len(linear)},{len(quadratic)},{elapsed_time}")

この Amplify プログラムも SymPy プログラムとほぼ同一であり，異なる部分のみを以下に説明します．

  • gen.array を用いて，二値変数の 2 次元配列 x を定義しています．
  • ampilfy.one_hot関数 を用いることで， constraint を簡潔に記述しています．
  • objective の計算には，Amplifyの einsum を用いて高速化しています．この einsum では，固定した i および j に対して，以下の式を計算します．C++で実装されたバックエンドが呼び出されるため，処理時間が大幅に短縮されます．
\sum_{k=0}^{s-1}\sum_{l=0}^{s-1}x_{i,k}x_{j,l}\mathrm{dist}_{k,l}
  • to_unconstrained_poly と as_dict 関数を用いて，QUBO 式 qubo を辞書形式で取得しています．

JijModelingによるプログラム

JijModelingを用いて，QAPをQUBOに変換するプログラムです．SymPyで用いたobjectiveを求めるsumを使った式をそのまま用いることができますが，2次項を式に一つづつ追加する処理が繰り返されるため，著しく処理速度が低下します．高速に処理するには，JijModelingが提供するjm.sum関数を用いて，効率的にモデルを記述する必要があります．

import time
import argparse
import jijmodeling as jm
import jijmodeling_transpiler as jmt

[省略]

start_time = time.perf_counter()

size_p = jm.Placeholder("size")
flow_size_p = jm.Placeholder("flow_size")
dist_p = jm.Placeholder("dist", shape=(size_p, size_p))
flow_p = jm.Placeholder("flow", shape=(flow_size_p, 3))
x = jm.BinaryVar("x", shape=(size_p, size_p))
i = jm.Element("i", belong_to=(0, size_p))
j = jm.Element("j", belong_to=(0, size_p))
k = jm.Element("k", belong_to=(0, size_p))
l = jm.Element("l", belong_to=(0, size_p))
f = jm.Element("f", belong_to=(0, flow_size_p))
problem = jm.Problem("QAP")
problem += jm.Constraint("row_sum",
                         jm.sum(j, x[i][j]) == 1, forall=i)
problem += jm.Constraint("column_sum",
                         jm.sum(i, x[i][j]) == 1, forall=j)
problem += jm.sum(f, jm.sum(k, jm.sum(l,
                  x[flow_p[f][0]][k] * x[flow_p[f][1]][l] * flow_p[f][2] * dist_p[k][l])))


compiled_model = jmt.core.compile_model(
    problem, {"size": size, "flow_size": len(flow), "dist": dist, "flow": flow})
qubo = jmt.core.pubo.transpile_to_pubo(compiled_model, normalize=False).get_qubo_dict(
    multipliers={'row_sum': 10000, 'column_sum': 10000})

elapsed_time = time.perf_counter() - start_time

constant = int(qubo[1])
linear = {}
quadratic = {}
x_map = {val: key for key, val in compiled_model.var_map.var_map['x'].items()}
for key, value in qubo[0].items():
    if (key[0] == key[1]):
        linear[x_map[key[0]]] = int(value)
    else:
        if (x_map[key[0]] < x_map[key[1]]):
            quadratic[(x_map[key[0]], x_map[key[1]])] = int(value)
        else:
            quadratic[(x_map[key[1]], x_map[key[0]])] = int(value)

if output:
    print(constant)
    for key, value in sorted(linear.items()):
        print(f"x[{key[0]}][{key[1]}] {value}")
    for key, value in sorted(quadratic.items()):
        print(
            f"x[{key[0][0]}][{key[0][1]}]*x[{key[1][0]}][{key[1][1]}] {value}")

print(f"{size},{degree},{len(linear)},{len(quadratic)},{elapsed_time}")

このJijModelingプログラムは，SymPyなどの一般的な数式処理ライブラリとは異なり，独特な記述スタイルを持っています．以下に主要な部分を説明します．

  • jm.Placeholderを用いて，以下のプレイスホルダを定義します．
    • size_p：sizeを与えます．
    • flow_size_p：flowの要素数を与えます．
    • dist_p：distを与えます．
    • flow_p：flowを与えます．
  • jm.BinaryVarを使って，二値変数の２次元配列xを定義します．
  • jm.Elementを使って，変数i,j,k,lの取る範囲を指定します．これらはxやdist_pのインデックスを指定するのに用います．
  • jm.Elementを使って，変数fの取る範囲を指定します．これはflow_pのインデックスを指定するのに用います．
  • jm.Problemを使って，QUBO問題全体を格納するモデルproblemを定義します．
  • jm.Constraintとjm.sumを用いて，各行・各列の合計が1になる制約をproblemに加えます．
  • jm.sumの3重ネストを用いて，2次項を加えています．
  • core.compile_modelを用いて，プレイスホルダに，際のsize，flowの要素数，dist，およびflowを割り当てた，compiled_modelを生成します．-
  • compiled_modelにtranspile_to_pubo，get_qubo_dict順に呼び出すことで，最終的なQUBO式を生成します．
  • compiled_model.var_map.var_mapは，変数xのインデックスをキーとし，その通し番号を値とする辞書です．求めたquboには，通し番号による一次項や二次項が格納されているので，x[0][0]などの表示を行うためには，この辞書を逆引きする必要があり，逆引き辞書をx_mapとしています．

QUBO++によるプログラム

QUBO++を用いて，QAPをQUBOに変換するプログラムです．

#include <boost/program_options.hpp>
#include "qbpp.hpp"

namespace po = boost::program_options;

inline uint32_t fixed_random(uint32_t n) {
  static uint32_t seed = 123456789;
  seed = seed * 1664525 + 1013904223;
  return seed % n;
}

inline std::vector<uint32_t> select(uint32_t size, uint32_t index,
                                    uint32_t degree) {
  if (!(degree < size)) {
    throw std::runtime_error("degree must be less than size.");
  }

  std::vector<uint32_t> work;
  work.reserve(size - 1);
  for (uint32_t i = 0; i < size; i++) {
    if (i != index) work.push_back(i);
  }

  std::vector<uint32_t> result;
  result.reserve(degree);

  for (uint32_t j = 0; j < degree; j++) {
    uint32_t k = fixed_random(static_cast<uint32_t>(work.size()));
    result.push_back(work[k]);
    std::swap(work[k], work.back());
    work.pop_back();
  }

  return result;
}

int main(int argc, char **argv) {
  po::options_description desc("Program to generate a QUBO for a sample QAP.");

  desc.add_options()("help,h", "Display this help message.")(
      "size,s", po::value<uint32_t>(), "Size of QUBO.")(
      "degree,d", po::value<uint32_t>(), "Degree of flow graph.")(
      "output,o", "Output the resulting QUBO.");

  po::variables_map vm;
  try {
    po::store(po::parse_command_line(argc, argv, desc), vm);
  } catch (const std::exception &e) {
    std::cout << "Wrong arguments. Please use -h/--help option to see the "
                 "usage.\n";
  }
  po::notify(vm);

  if (vm.count("help") || vm["size"].empty() || vm["degree"].empty()) {
    std::cout << desc << std::endl;
  }

  const uint32_t size = vm["size"].as<uint32_t>();
  const uint32_t degree = vm["degree"].as<uint32_t>();

  std::vector<std::vector<uint32_t>> dist(size, std::vector<uint32_t>(size, 0));
  for (uint32_t i = 0; i < size; i++) {
    for (uint32_t j = 0; j < size; j++) {
      if (i != j) dist[i][j] = fixed_random(100) + 1;
    }
  }
  std::vector<std::tuple<uint32_t, uint32_t, uint32_t>> flow;
  for (uint32_t i = 0; i < size; i++) {
    std::vector<uint32_t> selected = select(size, i, degree);
    for (const auto &j : selected) {
      flow.emplace_back(i, j, fixed_random(100) + 1);
    }
  }

  auto start = qbpp::get_time();

  auto x = qbpp::var("x", size, size);
  auto qubo = (qbpp::sum(qbpp::vector_sum(x) == 1) +
               qbpp::sum(qbpp::vector_sum(qbpp::transpose(x)) == 1)) *
              10000;
  for (const auto &[i, j, f] : flow) {
    for (uint32_t k = 0; k < size; k++) {
      for (uint32_t l = 0; l < size; l++) {
        qubo += x[i][k] * x[j][l] * f * dist[k][l];
      }
    }
  }

  qubo.simplify_as_binary();

  auto elapsed_time = qbpp::get_time() - start;

  size_t linear = 0;
  size_t quadratic = 0;
  if (vm.count("output")) std::cout << qubo.constant() << std::endl;
  for (const auto &term : qubo.terms()) {
    if (term.var_count() == 1) {
      ++linear;
      if (vm.count("output"))
        std::cout << term.vars()[0] << " " << term.coeff() << std::endl;
    } else if (term.var_count() == 2) {
      ++quadratic;
      if (vm.count("output"))
        std::cout << term.vars()[0] << "*" << term.vars()[1] << " "
                  << term.coeff() << std::endl;
    }
  }
  std::cout << size << "," << degree << "," << linear << "," << quadratic << ","
            << elapsed_time << std::endl;
}

このQUBO++プログラムの構成を簡単に説明します．

  • fixed_random(uint32_t n)：オリジナルの疑似乱数生成器であり，固定されたシードから一様なランダム整数（0以上n未満）を生成します．
  • select(uint32_t size, uint32_t index, uint32_t degree)：各 i に対して，\mathrm{flow}_{i,j} に値を割り当てる d 個の j をランダムに選ぶ関数です．
  • 二重の for ループで，\mathrm{dist}_{i,j} および \mathrm{flow}_{i,j} にランダムな値を設定しています．
  • qbpp::var により，2次元の二値変数 x を定義しています．
  • qbpp::vector_sumで2次元の二値変数xの各行の和を求めたベクトルを返しています．そしてそれらが1と一致するときに最小値0となる式を返し，qbpp::sumでその和を求めています．
  • qbpp::transposeで2次元の二値変数xを転置しているので，2次元の二値変数xの各列の和が1と一致するときに最小値が0となるとなる式を返します．
  • 四重の for ループを用いて目的関数を計算しquboに代入しています．このとき，\mathrm{flow}_{i,j} = 0 の場合は計算をスキップするようにしています．
  • 各行各列の合計が1となる制約条件をquboに加えています．
  • simplify_as_binary() 関数を用いて，重複する項の統合およびソートを行い，式を整理します．

PyQBPPによるプログラム

PyQBPPは，本評価を行うにあたり，QUBO++と互換性をもたせることを目的として急遽開発した，Pythonベースの軽量なツールです．QUBO++と同じアルゴリズムをPythonで実装した場合に，どの程度の性能が得られるかに関心があり，試作的に開発しました．そのため，現時点ではQUBO++の機能の一部のみをサポートしていますが，構文や動作の互換性を保っています．

以下に，PyQBPPを用いてQAPをQUBOに変換するプログラムを示します．

import time
import argparse
import pyqbpp as qbpp

[省略]

start_time = time.perf_counter()

x = qbpp.var("x", size, size)

constraint = sum(qbpp.vector_sum(x) == 1) + \
    sum(qbpp.vector_sum(qbpp.transpose(x)) == 1)

objective = sum(
    x[i][k] * x[j][l] * f * dist[k][l]
    for i, j, f in flow
    for k in range(size)
    for l in range(size)
)

qubo = objective + 10000 * constraint

qubo.simplify_as_binary()

elapsed_time = time.perf_counter() - start_time

constant = qubo.const
linear = {}
quadratic = {}
for term in qubo.terms:
    if len(term.vars) == 1:
        linear[f"{term.vars[0]}"] = term.coeff
    elif len(term.vars) == 2:
        quadratic[f"{term.vars[0]}*{term.vars[1]}"] = term.coeff
if output:
    print(constant)
    for key, value in sorted(linear.items(), key=str):
        print(f"{key} {value}")
    for key, value in sorted(quadratic.items(), key=str):
        print(f"{key} {value}")

print(f"{size},{degree},{len(linear)},{len(quadratic)},{elapsed_time}")

QUBO++の書式に従い，qbpp.vector_sumとqbpp.transposeを用いて制約constraintを定義しています．また，simplify_as_binary関数を呼び出すことで，式中の項の整理を行っています．

性能評価

SymPy，SymEngine, qubovert, TYTAN, PyQBPP, dimod, PyQUBO，JijModeling, Amplify，およびQUBO++について，d=5に固定し，施設数および場所数 s を変化させながら，QAPのQUBO式を生成するのに要する処理時間 t を計測しました．各条件で10回実行し，その平均処理時間を求めています．tが1000秒を超える程度の n まで測定を行いました．実行環境および使用条件は以下の通りです．

  • 4 x AMD EPYC 7702 64-Core Processor @ 1.50GHz
  • 2.0 Tbyte memory
  • Ubuntu 20.04.6
  • g++ 13.1.0
  • Python 3.12.10
  • sympy 1.14.0
  • symengine 0.14.1
  • tytan 0.1.1
  • qubovert 1.2.5
  • dimod 0.12.20
  • pyqubo 1.5.0
  • jijmodeling 1.12.4
  • jijmodeling-transpiler 0.7.3
  • amplify 1.3.1
  • qubo++ 2025.5.12
  • pyqbpp 2025.5.11

QAPにおける二値変数の個数は，s \times sの行列を用いるため n = s^2 となります．また，d = 5 とした場合，項の総数 T はおおよそ 5s^3 を超えたものになります．ここでは，1秒あたりに処理できる項の個数である T/t（Terms Per Second: TPS）を性能指標として用います．T/t が s の増加に伴って減少する場合，QUBO式の生成に要する時間が項数に対して線形以上の割合で増加していることを示しており，望ましくない挙動と評価されます．

補足
処理時間 t は有効数字3桁に丸めて表示していますが，TPS（T/t）は元の処理時間データに基づいて計算しているため，丸めた値を用いて再計算すると一致しない場合があります．

SymPyの性能

SymPyについては，s = 6から15までの範囲で計測を行いました．SymPyは汎用的な記号処理ツールであるため，1つの項を追加する際にも，既存の式の長さに依存した時間がかかると考えられます．そのため，項の数が増えるにつれてTPS（Terms Per Second）は急激に低下します．

例えば，s = 15のとき，項の総数は約1万5千ですが，TPSはわずか10.6であり，1秒間にわずか約10個の項しか生成できていないことになります．この結果から，SymPyではごく小規模なQUBO問題しか扱えないことが明らかです．

もっとも，SymPyは豊富な機能を備えた汎用記号処理ツールであり，ごく小さなQUBO問題を試験的に扱う場合に限っては，選択肢の一つとみなすこともできます．

s(facilities) n(variables) T(terms) t(sec) T/t(TPS)
6 36 666 3.11 214
7 49 1,225 7.4 165
8 64 1,912 15.7 121
9 81 2,889 34.1 84.8
10 100 4,150 70.7 58.7
11 121 5,951 164 36.3
12 144 7,668 342 22.4
13 169 10,309 508 20.3
14 196 12,936 932 13.9
15 225 15,765 1480 10.6

SymEngineの性能

SymEngineは，SymPyほどではありませんが，nが増加するにつれてTPSが急激に悪化します．特に，s = 24（変数数576，項数約7万）を超えると，1秒あたりの項数生成速度が10^3を大幅に下回り，処理時間が急激に増加する傾向が見られます．

これは，SymEngineがSymPyより高速なC++ベースの記号処理エンジンであるにもかかわらず，記号展開と項整理の過程において，式の項数の増加に比例して処理オーバーヘッドが増大されるためと考えられます．

このことから，SymEngineはSymPyに比べれば数段高速ですが，QUBO式の生成において大規模問題に十分耐えうる性能があるとは言いがたく，実用的な項数上限は数万程度までと考えられます．

s(facilities) n(variables) T(terms) t(sec) T/t(TPS)
6 36 666 0.056 11.9 \times 10^3
8 64 1,912 0.31 6.17 \times 10^3
10 100 4,150 1.25 3.33 \times 10^3
12 144 7,668 2.87 2.67 \times 10^3
14 196 12,936 6.44 2.01 \times 10^3
16 256 20,416 13.9 1.47 \times 10^3
18 324 28,782 27.2 1.06 \times 10^3
20 400 40,680 56.4 0.721 \times 10^3
22 484 56,848 126 0.45 \times 10^3
24 576 73,440 305 0.24 \times 10^3
26 676 94,276 676 0.139 \times 10^3
28 784 119,476 1380 0.0867 \times 10^3

qubovertの性能

qubovertは，SymPyよりおおよそ2倍程度高速ですが，SymPyと同様に，項数の増加に伴ってTPS（Terms Per Second）は急激に低下します．たとえば，s=17のときにはTPSがわずか22.1まで落ち込み，実用的な規模のQUBO式の生成には適さないことが分かります．

s(facilities) n(variables) T(terms) t(sec) T/t(TPS)
6 36 666 2.36 283
7 49 1,225 4.97 246
8 64 1,912 9.99 191
9 81 2,889 21.4 135
10 100 4,150 39.1 106
11 121 5,951 75.4 78.9
12 144 7,668 128 59.8
13 169 10,309 209 49.3
14 196 12,936 324 39.9
15 225 15,765 491 32.1
16 256 20,416 760 26.9
17 289 24,769 1120 22.1

TYTANの性能

TYTANは記号処理エンジンとしてSymEngineを利用しているため，処理時間の傾向はSymEngineとほぼ同様になります．ただし，SymEngineを直接呼び出す場合に比べて，TYTANでは一貫してわずかに処理時間が長くなっています．

これは，TYTANがSymEngineをラップする形で機能を提供しているため，内部でのオーバーヘッド（ラッパー層の処理や，Pythonインタフェースの抽象化）が影響していると考えられます．実際，TPS（1秒あたりに生成できる項数）を比較すると，全体的にSymEngineよりも10〜30％程度低下しています．

s(facilities) n(variables) T(terms) t(sec) T/t(TPS)
6 36 666 0.111 5.99 \times 10^3
8 64 1,912 0.427 4.48 \times 10^3
10 100 4,150 1.32 3.14 \times 10^3
12 144 7,668 3 2.56 \times 10^3
14 196 12,936 6.73 1.92 \times 10^3
16 256 20,416 15.1 1.35 \times 10^3
18 324 28,782 29.9 0.962 \times 10^3
20 400 40,680 62.7 0.649 \times 10^3
22 484 56,848 140 0.407 \times 10^3
24 576 73,440 327 0.225 \times 10^3
26 676 94,276 709 0.133 \times 10^3
28 784 119,476 1440 0.0832 \times 10^3

dimodの性能

dimodは，D-Wave Systems社が開発した量子アニーラー向けのソフトウェアツールであり，QUBOやIsingモデルの定式化・変換・操作をPythonで行うためのライブラリです．もともと量子アニーラー上で実行するためのモデル構築を主眼としているため，扱う二値変数の規模も数百〜1000個程度を想定しており，大規模問題への対応や記号展開の最適化にはあまり注力されていないと考えられます．

実際に，s=30を超えるあたりからTPS（1秒あたりの生成項数）は急激に低下し，s = 60ではTPSが500を下回る結果となっています．これは，内部処理が主にPythonレベルで完結しており，データ構造や操作の最適化が限定的であるためと推測されます．

それでも，SymEngineなどの記号処理系よりは高いパフォーマンスを示しており，数十万項までの中規模QUBO定式化においては比較的実用的な速度を保っています．一方で，100万項を超える規模のモデルを構築する用途には不向きであることが，この結果からも明らかです．

s(facilities) n(variables) T(terms) t(sec) T/t(TPS)
10 100 4,150 0.661 6.28 \times 10^3
15 225 15,765 2.26 6.97 \times 10^3
20 400 40,680 5.39 7.55 \times 10^3
25 625 82,225 13.7 6.01 \times 10^3
30 900 147,930 31.9 4.63 \times 10^3
35 1,225 236,845 77.3 3.07 \times 10^3
40 1,600 357,280 201 1.78 \times 10^3
45 2,025 500,985 429 1.17 \times 10^3
50 2,500 708,100 867 0.816 \times 10^3
60 3,600 1,224,900 2630 0.466 \times 10^3

PyQBPPの性能

PyQBPPは急遽開発されたにもかかわらず，SymEngineをはるかに上回る処理性能を示しており，dimodと同等のパフォーマンスを実現しています．特に小規模〜中規模の問題においては，1秒あたり数千項の生成が可能であり，記号処理ベースのツールとしてはある程度高速です．

ただし，nが増加するにつれてTPS（Terms Per Second）は徐々に低下していきます．これは，式の統合・整理処理においてPython標準のsort関数を使用しており，その計算コストがデータ量の増加とともに非線形に増大するためと考えられます．

PyQBPPの最大の利点は，QUBO++とほぼ同一の構文でQUBO式をPython上で記述できる点にあります．これにより，まずはPython上で小規模なデータを用いて開発・デバッグを行い，そのまま大規模な実行環境ではQUBO++に置き換えることで，記述を変更することなく高速処理が可能となります．

このように，PyQBPPはQUBO++の前段階として，記述の試作・検証やモデル設計のテストを効率的に行うための軽量なプロトタイピングツールとして非常に有用です．

s(facilities) n(variables) T(terms) t(sec) T/t(TPS)
10 100 4,150 0.123 33.6 \times 10^3
15 225 15,765 1 15.8 \times 10^3
20 400 40,680 4.41 9.22 \times 10^3
25 625 82,225 13.9 5.93 \times 10^3
30 900 147,930 36.2 4.09 \times 10^3
35 1,225 236,845 106 2.23 \times 10^3
40 1,600 357,280 295 1.21 \times 10^3
45 2,025 500,985 699 0.716 \times 10^3
50 2,500 708,100 1450 0.489 \times 10^3

JijModelingの性能

JijModelingでは，TPS（Terms Per Second）が比較的高い水準で安定しており，中規模程度までのQUBO式を効率よく構築できるツールと言えます．特に問題サイズが大きくなるにつれてもTPSが大きく低下せず，処理時間がほぼ線形に伸びる傾向にあることから，一定のスケーラビリティを持っていると評価できます．

この性能は，jm.Constraintやjm.sumなどの専用の組み込み関数によって問題構造を抽象的に記述し，その後にjmt.core.compile_modelを通じてRustベースのバックエンドで高速に展開・変換される設計に支えられていると考えられます．

なお，jm.Constraintやjm.sumの処理時間は問題サイズに依存せず常に約0.05秒程度であり，なんらかのモデルの構築が行われているわけではなく，compile_modelの実行時にQUBO式をRustで生成するための操作手順を記録している構文テンプレートを構成していると推測されます．つまり，JijModelingは，Python上で操作手順の構築を行い，後段のRust実装によって実際にQUBO式を構築する遅延評価型の設計を採用していると考えられます．

このような構成により，Pythonの制約を回避しながら，中規模までのQUBO式を比較的効率的に生成できる点がJijModelingの特徴です．

s(facilities) n(variables) T(terms) t(sec) T/t(TPS)
10 100 4,150 0.314 13.2 \times 10^3
20 400 40,680 1.66 24.5 \times 10^3
30 900 147,930 4.12 35.9 \times 10^3
40 1,600 357,280 8.82 40.5 \times 10^3
50 2,500 708,100 16.6 42.7 \times 10^3
60 3,600 1,224,900 27.9 43.9 \times 10^3
70 4,900 1,961,050 44.2 44.4 \times 10^3
80 6,400 2,976,800 66.2 45 \times 10^3
90 8,100 4,253,400 94.2 45.2 \times 10^3
100 10,000 5,841,100 128 45.7 \times 10^3
120 14,400 10,096,080 222 45.5 \times 10^3
140 19,600 16,171,400 356 45.5 \times 10^3
160 25,600 24,142,720 536 45.1 \times 10^3
180 32,400 34,507,800 759 45.5 \times 10^3
200 40,000 47,242,800 1050 44.9 \times 10^3

PyQUBOの性能

PyQUBOは，JijModelingと比較して一貫して約2倍のTPS（Terms Per Second）を達成しており，中規模程度までのQUBO式の生成において高い性能を示しています．特筆すべきは，この性能が特殊な組み込み関数を直接使うことなく，通常のPython構文による素直な記述だけで得られている点です．

具体的には，Python上で演算子のオーバーロードを活用してQUBO式を構築し，内部ではC++で実装されたバックエンドが効率的に式展開・項の整理を行っています．そのため，ユーザーは複雑なDSLを学ばずに直感的なコードでQUBOを記述でき，なおかつ比較的高速に変換処理が実行されます．

また，問題サイズが増加してもTPSが安定しており，s = 240においても1秒あたり7万項以上を維持しています．このことから，中規模程度までのQUBO定式化においては，実用的でバランスの取れたツールであると評価できます．

s(facilities) n(variables) T(terms) t(sec) T/t(TPS)
10 100 4,150 0.118 35.3 \times 10^3
20 400 40,680 0.824 49.4 \times 10^3
30 900 147,930 2.34 63.2 \times 10^3
40 1,600 357,280 4.71 75.8 \times 10^3
50 2,500 708,100 8.39 84.4 \times 10^3
60 3,600 1,224,900 14.4 85.1 \times 10^3
70 4,900 1,961,050 23.3 84.2 \times 10^3
80 6,400 2,976,800 35.5 83.9 \times 10^3
90 8,100 4,253,400 49.8 85.4 \times 10^3
100 10,000 5,841,100 68.4 85.4 \times 10^3
120 14,400 10,096,080 119 84.7 \times 10^3
140 19,600 16,171,400 194 83.2 \times 10^3
160 25,600 24,142,720 292 82.7 \times 10^3
180 32,400 34,507,800 426 80.9 \times 10^3
200 40,000 47,242,800 603 78.3 \times 10^3
240 57,600 81,967,680 1140 71.9 \times 10^3

Amplifyの性能

Amplifyは，PyQUBOと比較して約10倍のTPS（Terms Per Second）を達成しており，中規模のQUBO式を高速に構築できるツールです．1秒あたり数十万項を安定して生成でき，記号処理ベースのツールとしては非常に高いスループットを誇ります．

この高性能の要因は，Amplifyが独自に提供するeinsum関数にあると推測されます．これはNumPyのeinsumに似た構文を持ちながら，二値変数を含むテンソル演算を記号的に処理し，内部的にQUBO項へ効率的に変換する仕組みを備えているためと思われます．ただし，n\geq 320でTPSが明らかに下っていく傾向があり，さらなる大規模化において実用上の障害になる可能性があります．

Pythonベースのツールとしては最高の性能を持ち，Pythonの利用が前提となる環境においては，最適な選択肢の一つと言えます．

s(facilities) n(variables) T(terms) t(sec) T/t(TPS)
10 100 4,150 0.0112 369 \times 10^3
20 400 39,594 0.0549 721 \times 10^3
30 900 144,550 0.176 823 \times 10^3
40 1,600 354,744 0.435 815 \times 10^3
50 2,500 703,145 0.913 770 \times 10^3
60 3,600 1,217,893 1.48 825 \times 10^3
70 4,900 1,953,177 2.47 792 \times 10^3
80 6,400 2,971,706 3.53 842 \times 10^3
90 8,100 4,240,557 4.76 891 \times 10^3
100 10,000 5,830,924 6.91 844 \times 10^3
120 14,400 10,084,329 11.3 893 \times 10^3
140 19,600 16,159,975 19.1 846 \times 10^3
160 25,600 24,135,352 29.2 826 \times 10^3
180 32,400 34,486,171 40.8 846 \times 10^3
200 40,000 47,182,047 59.5 793 \times 10^3
240 57,600 81,908,535 102 803 \times 10^3
280 78,400 130,756,655 165 793 \times 10^3
320 102,400 194,035,468 274 707 \times 10^3
360 129,600 277,989,966 401 693 \times 10^3
400 160,000 380,629,077 615 618 \times 10^3
450 202,500 541,891,917 937 578 \times 10^3
500 250,000 745,000,093 1390 537 \times 10^3

QUBO++の性能

QUBO++は，すべての問題サイズにおいてAmplifyよりも高いTPS（Terms Per Second）を達成しており，問題規模が大きくなるにつれてその差はさらに顕著になります．とくにn \geq 180の領域では10倍以上の差が生じ，n \geq 450では20倍以上に達します．

小規模な問題では，TPSは相対的に低くなる傾向がありますが，これはデータ量が少なすぎて，QUBO++が採用しているマルチスレッド処理の並列化効果が十分に発揮されないためと考えられます．一方で，問題サイズの増加に伴い，並列処理の効果が顕著になり，TPSは上昇し続ける傾向を示しています．処理時間を測定した最大規模（n=1400）においても，TPSの減少は見られません．

QUBO++の大きな利点は，式の展開，項の整理，および統合処理のすべてがC++で実装されており，Pythonインタプリタによるオーバーヘッドが一切存在しない点にあります．さらに，記号処理は二値変数に特化して最適化されており，oneTBBを用いたマルチスレッド処理によって高い並列性能も実現しています．このような設計により，QUBO++は極めて高いスループットと優れたスケーラビリティを備えており，本評価におけるすべてのツールの中で最も高い性能を発揮しています．

QUBO式の構築性能が，システム全体の処理効率に直接影響を及ぼすような用途においては，QUBO++が最も有力な選択肢となるでしょう．

s(facilities) n(variables) T(terms) t(sec) T/t(TPS)
10 100 4,150 0.0078 0.532 \times 10^6
20 400 40,680 0.034 1.19 \times 10^6
30 900 147,930 0.0572 2.59 \times 10^6
40 1,600 357,280 0.0878 4.07 \times 10^6
50 2,500 708,100 0.131 5.41 \times 10^6
60 3,600 1,224,900 0.22 5.57 \times 10^6
70 4,900 1,961,050 0.314 6.25 \times 10^6
80 6,400 2,976,800 0.446 6.68 \times 10^6
90 8,100 4,253,400 0.612 6.95 \times 10^6
100 10,000 5,841,100 0.779 7.5 \times 10^6
120 14,400 10,096,080 1.29 7.83 \times 10^6
140 19,600 16,171,400 1.85 8.76 \times 10^6
160 25,600 24,142,720 2.95 8.19 \times 10^6
180 32,400 34,507,800 4.05 8.51 \times 10^6
200 40,000 47,242,800 4.97 9.5 \times 10^6
240 57,600 81,967,680 7.35 11.2 \times 10^6
280 78,400 130,851,280 10.4 12.6 \times 10^6
320 102,400 194,156,480 15.7 12.4 \times 10^6
360 129,600 278,124,840 24.5 11.3 \times 10^6
400 160,000 380,806,000 32.4 11.8 \times 10^6
450 202,500 542,100,600 41.3 13.1 \times 10^6
500 250,000 745,257,000 58.8 12.7 \times 10^6
550 302,500 993,114,100 71.3 13.9 \times 10^6
600 360,000 1,289,887,200 95.4 13.5 \times 10^6
700 490,000 2,050,657,000 145 14.1 \times 10^6
800 640,000 3,061,129,600 220 13.9 \times 10^6
900 810,000 4,361,049,900 330 13.2 \times 10^6
1000 1,000,000 5,985,010,000 439 13.6 \times 10^6
1200 1,440,000 10,347,850,800 777 13.3 \times 10^6
1400 1,960,000 16,416,986,600 1160 14.2 \times 10^6

QUBO生成ツール別TPS（Terms Per Second）の比較

以下のグラフは，QAPのサイズ（施設配置数 s）に対するTPS（Terms Per Second，すなわち1秒間に生成される項数）を，各ツールごとに比較したものです．縦軸・横軸ともに対数スケールで表示しており，横軸の1目盛はサイズの2倍，縦軸の1目盛はTPSが10倍であることを意味します．

また，グラフには処理時間が10秒，100秒，1000秒となる等時間線を破線で表示しています．各ツールの曲線とこれらの破線が交差する点は，「その処理時間で対応できる最大の問題サイズ」と「そのときのTPSの目安」を示します．例えば，1000秒の等時間線とQUBP++は，s=1280ぐらいで交差しており，このことから，QUBO++はs=1280のQAPを1000秒程度で生成できると推測されます．

各ツールの性能傾向

  • QUBO++ は，sの増加とともにTPSが緩やかに上昇し，s \geq 240 では10^7TPSを超えるなど，他のツールを大きく引き離す性能を示しています．大規模な問題にも余裕をもって対応可能です．
  • Amplify は，すべてのサイズにおいておおむね 3 \times 10^5 ～ 9 \times 10^5 TPS を維持しており，Pythonベースのツールとしては最高の性能を発揮しています．
  • PyQUBO と JijModeling はともに 10^4 ～ 10^5 TPS の範囲で安定しており，PyQUBOの方が約2倍高い性能を示しています．中規模までの問題においては，実用的なツールです．
  • dimod と PyQBPP は類似した性能を示し，sの増加とともにTPSが急激に悪化します．この傾向から，実用的な上限は s = 40 ～ 50 程度と考えられます．
  • TYTAN は SymEngine をバックエンドに使用しているため，性能はほぼ同じです．SymEngineを直接用いた方がオーバーヘッドが少ない分，やや高速です．いずれも dimod や PyQUBO に比べて10倍以上遅く，TPSの低下も顕著なため，s = 20 ～ 25 が実用限界と推測されます．
  • SymPy および qubovert は性能傾向が近く，TPSはサイズの増加とともに急激に低下します．TYTANやSymEngineと比較してもさらに10倍以上低速であり，実用的に扱えるのは s = 15 以下のごく小規模な問題に限られます．

まとめ

QUBO式の生成時間を可能な限り短縮したい場合，QUBO++ が最も有力な選択肢となります．一方で，Pythonベースのツールを使用する必要がある場合には，Amplify が最適といえるでしょう．ただし，einsum のような高速化のための組み込み関数が適用しにくい問題や，それらの学習コストを避けたい場合には，PyQUBO を用いるのが現実的です．

それ以外のツールを選択する理由としては，そのツールが提供するソルバーとの連携が主目的である場合に限られると考えられます．

QUBO++情報

広島大学コンピューティングラボ
広島大学コンピューティングラボPublication

広島大学のコンピューティング関連分野で研究を行う研究者が情報を発信します． 主に，量子化学計算や数理最適化に関する技術や研究成果についての記事を掲載しています．

Discussion