PyQUBOの使い方

PyQUBOでは、直感的な記述で定式化をプログラムで構築することができます。

量子アニーリングやイジングマシンを利用する際、最大の壁となるのが「数式からQUBO（行列）への変換」です。この煩雑な作業を自動化し、Pythonのコードとして数式を記述できるツールがPyQUBOです。

本記事では、PyQUBOで何ができるのか、そして実際のコードを用いた実装の流れを解説します。

1. PyQUBOとは？

PyQUBOは、QUBO（Quadratic Unconstrained Binary Optimization）やイジングモデルを構築するためのPythonライブラリです。

通常、アニーリングマシンに問題を解かせるには、QUBOを定式化して、自身で数式を展開してあげる必要がありますが、PyQUBOを使うことで、以下のように数式に近い形で記述するだけで、自動的に行列を生成してくれます。

PyQUBOでできること（主な特徴）

1. 数式に近い記述: Pythonの演算子を使って、ハミルトニアン（目的関数＋制約条件）を自然に記述できます。
2. 制約条件の自動処理: 制約条件（例：AとBは同時に選べない）をペナルティ項として表現する際、Constraint関数を使うことで、コンパイル後（計算直前）にペナルティの強さを自由に調整できます。再コンパイルの手間がありません。
3. 変数の柔軟な定義: 0/1のBinary変数だけでなく、1/-1のSpin変数、さらには整数を表すための表現（One-hotエンコーディング等）もサポートしています。
4. ソルバーへの橋渡し: 作成したモデルは .to_qubo() や .to_ising() メソッドで、OpenJijやD-Wave Ocean SDK、D-Wave Leapなどが受け取れる形式へ簡単に変換できます。

---

2. 実践：PyQUBOによるモデリングの流れ

ここから、実際に作成したプログラムの解説に入ります。対象とする問題は巡回セールスマン問題(TSP)です。

$$\begin{align*} \text{min}\quad & \sum_{t=0}^{n-1}\sum_{i=0}^{n-1}\sum_{j=0}^{n-1} d_{i,j} x_{i,t} x_{j,t+1} \\ \text{s.t.}\quad & \sum_{i=0}^{n-1} x_{i,t} = 1 \quad (\forall t) \\ & \sum_{t=0}^{n-1} x_{i,t} = 1 \quad (\forall i) \end{align*}$$

• 目的関数: 総移動距離の最小化
• 制約: 各時刻に1拠点のみ訪れる & 各拠点は1回のみ訪れる

PyQUBOを使ったフローは主に以下の4ステップです。

1. 変数の定義: 数式に使う変数を用意する
2. ハミルトニアンの定義: 目的関数と制約条件を作る
3. コンパイル: モデルを計算可能な状態にする
4. QUBO変換: ソルバーに渡せる辞書型データにする

Step 1: ライブラリのインポートと変数の定義

まずは必要な変数を定義します。

# 必要なライブラリ:
# pip install pyqubo dwave-ocean-sdk numpy

import numpy as np
from pyqubo import Array, Constraint
from dwave.system.samplers import DWaveSampler
from dwave.system.composites import EmbeddingComposite

# 1. データの準備
distance_matrix = np.array([
    [0,  10, 20, 15],
    [10, 0,  15, 30],
    [20, 15, 0,  10],
    [15, 30, 10, 0]
])
num_cities = 4
penalty_strength = 50.0

# 2. PyQUBOによる数式モデルの構築
# 変数の作成: 4x4 のバイナリ変数配列 'x' を作成
# x[i, j] が「都市i が j番目」に対応します
x = Array.create('x', shape=(num_cities, num_cities), vartype='BINARY')

解説のポイント:

• Array.create('x', shape=N, vartype='BINARY') などを使って、配列変数を一括定義している箇所を説明します。
• Binary（0/1）なのか Spin（-1/+1）なのか、問題設定に合わせて選んでいることを伝えます。

Step 2: ハミルトニアン（数式）の記述

解きたい問題の目的関数と、守るべき制約条件を記述します。

# --- 目的関数の記述 ---
#  総距離の最小化
H_obj = 0.0
for t in range(num_cities):
    next_t = (t + 1) % num_cities
    for c1 in range(num_cities):
        for c2 in range(num_cities):
            if c1 != c2:
                dist = distance_matrix[c1][c2]
                # 都市c1(時点t) と 都市c2(時点next_t) が両方1なら距離を加算
                H_obj += dist * x[c1, t] * x[c2, next_t]

# --- 制約条件の記述 ---
# PyQUBOでは (合計 - 1)**2 のように直感的に書けます
# Constraint() で囲むことで、後で制約が守られているかチェックしやすくなります

H_const = 0.0

#  縦の制約 (Order制約)
# 「各順番 j において、選ばれる都市 i の合計は1つ」
for j in range(num_cities):
    # sum(x[i, j] for i in range(...)) は列の合計
    # ※PyQUBOのSumではなく、Python標準のsum()を使います
    col_sum = sum(x[i, j] for i in range(num_cities))
    H_const += Constraint((col_sum - 1)**2, label=f"Order_{j}")

#  横の制約 (City制約)
# 「各都市 i において、選ばれる順番 j の合計は1回」
for i in range(num_cities):
    # sum(x[i, j] for j in range(...)) は行の合計
    row_sum = sum(x[i, j] for j in range(num_cities))
    H_const += Constraint((row_sum - 1)**2, label=f"City_{i}")

# 全ハミルトニアン（総エネルギー式）の作成
# 制約項にはペナルティ係数を掛けます
H = H_obj + penalty_strength * H_const

解説のポイント:

• 目的関数: 「何を最小化（最大化）したいか」をコードと対比させます。
• 制約条件: Constraint(式, label="...") を使っている場合、ここがPyQUBOのハイライトです。「制約を守らせるための項」であることを説明します。
• Placeholder を使っている場合は、「後から値を代入できる変数」として紹介すると親切です。

Step 3: コンパイルとQUBO/Isingへの変換

記述した数式をモデルオブジェクトにコンパイルし、ソルバーが理解できる形式（QUBOまたはIsing）に変換します。

# 3. QUBOへの変換 & 実行
print("数式をコンパイル中...")
model = H.compile()
# D-Waveが理解できる辞書形式のQUBOデータに変換
Q_dict, offset = model.to_qubo()

解説のポイント:

• .compile() することで、数式が行列計算の準備状態になることを説明します。
• feed_dict で Placeholder に具体的な数値（ペナルティ係数など）を与えている場合、ここで「制約の強さを調整している」と解説します。
• 返り値の qubo （辞書型）と offset （定数項）が、ソルバーへの入力になります。

Step 4: ソルバーで実行してみる

PyQUBOで書いたものをそのままD-Waveマシンで回してみます。

# --- D-Wave設定 ---
# 注意: 'xxxx' の部分はご自身のAPI Token に書き換えてください
#token = 'xxxx'
sampler = EmbeddingComposite(DWaveSampler(token=token))

print("量子プロセッサでアニーリングを実行中...")
sampleset = sampler.sample_qubo(Q_dict, num_reads=1000)

# 最良の解を取得
best_sample = sampleset.first.sample

# ==========================================
# 4. 結果表示
# ==========================================
print("\n--- 結果解析 ---")

# PyQUBOのモデルを使って、生の解(0/1)を扱いやすい形にデコードします
decoded_sample = model.decode_sample(best_sample, vartype='BINARY')

# 制約が破れていないかチェック
broken_constraints = decoded_sample.constraints(only_broken=True)

if not broken_constraints:
    print("制約OK! 最適ルートが見つかりました。")

    # 辞書から直接ルートを復元（.array()メソッドはエラー回避のため使いません）
    route = [-1] * num_cities

    # decoded_sample.sample は {'x[0][0]': 1, 'x[0][1]': 0, ...} という形式の辞書です
    for i in range(num_cities): # 都市
        for j in range(num_cities): # 順番
            # 変数名 'x[i][j]' の値が 1 かどうかを確認
            key_name = f"x[{i}][{j}]"
            if decoded_sample.sample[key_name] == 1:
                route[j] = i

    print("ルート:", route)

    # 総距離の計算と表示
    total_dist = 0
    print("\n--- 移動詳細 ---")
    for i in range(num_cities):
        u, v = route[i], route[(i+1)%num_cities]
        dist = distance_matrix[u][v]
        total_dist += dist
        print(f"{u} -> {v}: {dist}km")
    print("-" * 20)
    print(f"総距離: {total_dist} km")

else:
    print("制約違反があります。")
    print("破れた制約:", broken_constraints.keys())

参考文献

[1] : PyQUBO, https://pyqubo.readthedocs.io/en/latest/, 閲覧日 2025/12/18