はじめに
量子アニーリング(D-Wave等)やイジングマシンで最適化問題を解いている方は、多くの場合すでに解きたい問題をQUBO行列の形に定式化済みだと思います。この資産をゲート方式の量子コンピュータ、特にQAOA(Quantum Approximate Optimization Algorithm)に持ち込んで解いてみたい、というのは自然な発想のはずです。
そこで今回は、「手元にあるQUBO行列を、そのままQiskitのQAOAに移植して解くにはどうすればいいか」を実際に手を動かして確認した記録を書きます。QAOAそのものの収束性能を研究するような記事ではなく、あくまで移植する側の実務的な視点で、
- QUBO行列からQiskitのコードへの変換は何行で済むのか
- アニーリング由来のQUBOをそのまま持ち込むときに踏みがちな罠は何か
- 移植した結果をそのまま信用していいのか、移植後に何を確認すべきか
を確認します。
TL;DR
- QUBO行列からQiskitへの移植は本質的に数行で済みます。行列の対角成分を線形項、非対角成分の和($Q_{ij}+Q_{ji}$)を二次項として
QuadraticProgramに登録するだけです - アニーリング由来のQUBOは制約なし・バイナリ変数のみなので、それ自体が既にQUBO形式です。
QuadraticProgramToQuboのような追加のQUBO変換処理は不要で、むしろ無駄な二重変換になります - 移植した直後の結果をそのまま信用せず、
reps(層数)や古典最適化器を振って安定して厳密解に近づくかを確認したほうがよさそうです。今回の例ではrepsを1→6に増やすと確率は36%→57%→55%→87%→84%→97%と改善しましたが単調ではなく、reps=4〜6はいずれもmaxiter=200の打ち切り時点でもまだ収束途中でした
環境
Windows + Python 3.10.7 + venvで検証しました。バージョンは以下のとおりです。
qiskit==1.2.4
qiskit-aer==0.15.1
qiskit-algorithms==0.3.1
qiskit-optimization==0.6.1
matplotlib==3.10.9
numpy==2.2.6
最初は最新のqiskit==2.5.1を入れましたが、この環境(Python 3.10.7)ではimport qiskitの時点で以下のエラーが出てクラッシュしました。
File ".../qiskit/passmanager/multistage_passmanager.py", line 216, in MultiStagePassManager
callback: Callback[Any] | None = None,
File "...\_collections_abc.py", line 483, in __getitem__
raise TypeError(f'Too {"many" if item_len > param_len else "few"}'
TypeError: Too few arguments for collections.abc.Callable[...]; actual 1, expected 2
qiskit.passmanager内のクラス定義がロード時にtyping.Callableのジェネリック引数を評価する際に起きている型ヒント関連の不具合で、この手のPython 3.10系との相性問題は珍しくありません。今回は原因を深追いせず、qiskit==1.2.4 + qiskit-algorithms==0.3.1 + qiskit-optimization==0.6.1という少し前の組み合わせに固定することで解決しました。移植先の環境を新規に作る場合は、同じ問題に当たる可能性があることを念頭に置いておくとよさそうです。
なおqiskit==1.2.4環境でも、後述のSampler()を使うたびに以下のDeprecationWarningが出続けます。
DeprecationWarning: The class ``qiskit.primitives.sampler.Sampler`` is deprecated as of qiskit 1.2.
It will be removed no earlier than 3 months after the release date. All implementations of the
`BaseSamplerV1` interface have been deprecated in favor of their V2 counterparts. The V2 alternative
for the `Sampler` class is `StatevectorSampler`.
以前書いたQiskit 2.0 チートシートでも触れたとおり、QiskitのprimitivesはV1(Sampler/Estimator)からV2(StatevectorSamplerなど)への移行が進んでいます。ただ今回使うqiskit-algorithms==0.3.1のQAOA(実体はSamplingVQE)はV1のBaseSamplerインタフェースを前提に実装されているため、あえて非推奨のSamplerをそのまま使っています。動作に支障はありませんが、qiskit-algorithmsが新しいV2対応版に上がったタイミングで書き換えが必要になりそうです。
QUBO行列をそのままQiskitに持ち込む
アニーリング側で扱うQUBOは、$n \times n$ 行列 $Q$ と
$$
\text{minimize} \ \ x^{\mathsf T} Q x, \ \ x_i \in {0, 1}
$$
という形で定式化されているはずです。今回の検証では例として以下の3変数QUBOを使います。
$$
Q = \begin{pmatrix} 1 & -2 & 0 \\ -2 & 1 & -1 \\ 0 & -1 & 1 \end{pmatrix}
$$
これをQiskit側のQuadraticProgramに変換するコードは以下です。行列を受け取ってバイナリ変数を登録し、対角成分を線形項、非対角成分の和を二次項として渡しているだけで、それ以上のことは何もしていません。
import numpy as np
from qiskit_optimization import QuadraticProgram
def build_qubo_from_matrix(Q: np.ndarray) -> QuadraticProgram:
n = Q.shape[0]
qp = QuadraticProgram(name="qubo")
for i in range(n):
qp.binary_var(name=f"x{i}")
linear = {f"x{i}": float(Q[i, i]) for i in range(n)}
quadratic = {}
for i in range(n):
for j in range(i + 1, n):
coeff = float(Q[i, j] + Q[j, i])
if coeff != 0.0:
quadratic[(f"x{i}", f"x{j}")] = coeff
qp.minimize(linear=linear, quadratic=quadratic)
return qp
これをQAOAで解く部分も含めた全体は以下のとおりです。
from qiskit_algorithms import QAOA
from qiskit_algorithms.optimizers import COBYLA
from qiskit_algorithms.utils import algorithm_globals
from qiskit.primitives import Sampler
from qiskit_optimization.algorithms import MinimumEigenOptimizer
def solve_with_qaoa(qp: QuadraticProgram, reps: int = 2, seed: int = 42):
algorithm_globals.random_seed = seed
sampler = Sampler()
optimizer = COBYLA(maxiter=200)
qaoa = QAOA(sampler=sampler, optimizer=optimizer, reps=reps)
meo = MinimumEigenOptimizer(qaoa)
return meo.solve(qp)
移植の要点は「QuadraticProgramに手元のQUBO行列をそのまま線形項・二次項として登録し、MinimumEigenOptimizer.solve(qp)に渡すだけ」という点に尽きます。制約の追加もペナルティ項の調整も要りません。手元のQUBO行列さえあれば、ここまでは数行で終わります。
実際に組まれる量子回路を見る
MinimumEigenOptimizer.solve(qp)はQUBOをそのまま量子回路に渡しているわけではなく、内部で$x_i = (1-Z_i)/2$という置き換えによってイジングハミルトニアンに変換してから回路を組んでいます。この変換の中身はQuadraticProgram.to_ising()で直接確認できます。今回のQUBOに対して実行すると以下が返ってきました。
SparsePauliOp(['IIZ', 'IZI', 'IZZ', 'ZZI'],
coeffs=[ 0.5+0.j, 1. +0.j, -1. +0.j, -0.5+0.j])
offset: 0.0
このコスト演算子をqiskit.circuit.library.QAOAAnsatz(cost_operator=cost_op, reps=1)に渡して展開すると、qiskit_algorithms.QAOAが内部で実際に組み立てているのと同じ構造の回路が得られます。reps=1の場合の回路は以下のとおりです。
回路は大きく3段構成になっています。
- 全qubitに
U2(0, π)(アダマールゲート相当)をかけて一様重ね合わせ状態を作る初期化 - コスト演算子由来の
RZ(1qubit項)とZZ結合(2qubit項)。パラメータγ(ガンマ)が入り、係数はQUBOの線形項・二次項をそのまま反映しています - 全qubitに
RX(2β)をかけるミキサー層(標準のXミキサー)
reps=2にすると、2.と3.のブロックがもう一段複製され、新たにγ[1], β[1]というパラメータが追加されます。
後述する「repsを1〜6まで振る」という検証は、まさにこの2.+3.のブロックを何段重ねるかという話だったことが、図から見て取れます。層を増やすほど回路は深くなり、古典最適化器が調整すべきパラメータ数($2 \times \text{reps}$個)も比例して増えます。
移植でやりがちな罠: QuadraticProgramToQuboは不要
実装の途中、MinimumEigenOptimizer.solve()に渡す前にQuadraticProgramToQubo().convert(qp)でQUBO形式に変換するコードを一度書いてしまいました。アニーリング側のツール(制約付きモデルをBQM/QUBOに変換する処理を挟むことが多い)に慣れていると、「Qiskit側でも念のためQUBOへの変換を挟んでおくべきでは」という発想が自然に出てきますが、これは今回のケースでは無駄な処理でした。
QuadraticProgramToQuboは本来、制約付きの問題(等式・不等式制約を持つQuadraticProgram)をペナルティ法で制約なしのQUBOへ変換するためのコンバータです。一方、アニーリングで既にQUBO行列に落とし込んである問題は、その時点で制約なし・バイナリ変数のみという条件を満たしています。つまりbuild_qubo_from_matrix()が作るQuadraticProgramは、変換するまでもなく最初から既にQUBOそのものになっています。ここにQuadraticProgramToQuboを挟むと「QUBO→QuadraticProgram(実質QUBO)→QuadraticProgramToQubo→QUBO」と、変換して同じ形式に戻すだけの二度手間になってしまいます。手元のQUBO行列を渡す場合、MinimumEigenOptimizer.solve()にはQuadraticProgramをそのまま渡せば十分です。 制約付きの問題(不等式制約などを含む、より一般的なQuadraticProgram)を扱うときだけ、QuadraticProgramToQuboが意味を持ちます。
移植した結果は、そのまま信用していいか
ここからが移植する側にとって一番重要な部分です。アニーリングでは「問題(QUBO)を投げれば、そこそこの低エネルギー状態のサンプルが返ってくる」という使い方ができますが、QAOAはそうではありません。QAOA回路には$\beta, \gamma$という変分パラメータがあり、これを古典最適化ループで訓練しながらコスト(ハミルトニアンの期待値)を下げていきます。つまり移植した瞬間からアニーリングと同じ感覚で結果を信用してよいわけではなく、「ちゃんと収束したか」「デフォルトのハイパーパラメータで十分か」を移植後に確認する必要があります。これを、全探索による厳密解と突き合わせながら確認しました。
まず厳密解を確保する
3変数程度の小さいQUBOなら、全探索($2^3=8$通り)で厳密解が一瞬で求まります。実運用の大きな問題ではこの方法は使えませんが、移植したパイプライン自体が正しく動いているかを検証する目的では、小さいテストケースで厳密解と突き合わせるのは有効な手段です。
import itertools
def exact_solve(qp: QuadraticProgram):
n = qp.get_num_binary_vars()
best_x, best_fval = None, np.inf
for bits in itertools.product([0, 1], repeat=n):
x = np.array(bits, dtype=float)
fval = qp.objective.evaluate(x)
if fval < best_fval:
best_fval = fval
best_x = x
return best_x, best_fval
実行結果は以下のとおりです。
x = [1. 1. 1.], fval = -3.0
MinimumEigenOptimizer.solve()が返すresult.xもこの厳密解と一致しました(x=[1,1,1], fval=-3.0)。ただしこれは要注意で、result.xは「サンプルした中で一番良かったビット列」を採用しているだけなので、たとえ低い確率でしか厳密解が出ていなくても、たまたま一度サンプルに含まれていればresult.xは厳密解を返してしまいます。移植後にresult.xだけを見て「解けた」と判断するのは危険で、実際にどれくらいの確率で厳密解が出ているのかを見る必要があります。
学習曲線を見る: ちゃんと収束しているか
qiskit_algorithms.QAOA(実体はSamplingVQE)にはcallback引数があり、callback(eval_count, params, mean, metadata)という形で各反復の期待値meanを受け取れます。これで最適化ループが実際に収束しているかを可視化できます。
cost_history = []
def callback(eval_count, params, mean, metadata):
cost_history.append(float(np.real(mean)))
qaoa = QAOA(sampler=sampler, optimizer=optimizer, reps=reps, callback=callback)
reps=2, COBYLA(maxiter=200)で実行した学習曲線が以下です。
72反復で収束し、最終的な期待値は<H>=-2.3482でした。厳密解の-3.0(赤破線)には届いていません。このときサンプル分布を見ると、厳密解x=[1,1,1]を測定する確率は**57.26%**でした。result.xはたまたまこの厳密解を返しましたが、それは「57%の確率で出る」という比較的手堅い当たりであって、100%ではありません。期待値$\langle H \rangle$が厳密値まで下がりきっていない場合、result.xが正解でも、それは何回か繰り返せば外れることもある「当たりくじ」に過ぎない、というのが移植者にとっての教訓になります。
reps(層数)を振ってみる
repsはQAOA回路の層数(コスト演算子とミキサー演算子を繰り返す回数)で、QAOA側で新たに触ることになるハイパーパラメータです。ただし「アニーリングにはハイパーパラメータがない」わけではない点は注意しておきたいところです。OpenJijのSASampler/SQASamplerにもnum_sweeps・num_reads・beta_min/beta_max(SAの逆温度スケジュール)、SQAならさらにgamma(横磁場強度)・trotter(トロッター数)といったパラメータがあり、D-Wave Ocean SDKにもnum_reads・annealing_time・chain_strength・anneal_scheduleがあります。JijZept SolverもOpenJijをアダプター経由で利用しつつtime_limit_secのような独自パラメータを持ち、Fujitsu Digital Annealerにも実行回数や温度スケジュール関連の設定項目があります。つまり移植したからといってハイパーパラメータ調整から解放されるわけではなく、チューニングすべきパラメータの種類が変わる、と捉えたほうが実態に近そうです。ここでは移植後にまず触ることになるrepsについて、COBYLA・maxiter=200固定でreps=1〜6を比較しました。
| reps | 収束までの反復数 | 最終<H>
|
厳密解を引く確率 |
|---|---|---|---|
| 1 | 39 | -1.3598 | 35.76% |
| 2 | 72 | -2.3482 | 57.26% |
| 3 | 200(打ち切り) | -2.0973 | 54.75% |
| 4 | 200(打ち切り) | -2.8465 | 87.46% |
| 5 | 200(打ち切り) | -2.4955 | 83.52% |
| 6 | 200(打ち切り) | -2.9536 | 97.06% |
reps=1が最も浅く、期待値も確率も最も悪い結果でした。reps=2は72反復で明確に収束し(グラフでも早期にプラトーに達しています)、reps=1より改善しています。reps=3はmaxiter=200まで反復してもグラフ上ではかなり早い段階(50反復あたり)で下降が緩やかになり、reps=2より良い期待値には届きませんでした。パラメータ空間が広がった分、COBYLA(局所探索の古典最適化器)がより浅い局所解に留まった可能性はありますが、これは今回1シード・1回の実行だけでの観測であり、シードや初期値を変えても同じ傾向になるかは未検証です(あくまで「そう見える」という記録です)。reps=4,5,6はいずれもグラフを見るとmaxiter=200で打ち切られた時点でもまだ右肩下がりが続いており、収束しきる前に反復上限に達しています。それでもreps=6は最終的な期待値<H>=-2.9536(厳密値-3.0にかなり近い)、確率97.06%まで到達しており、6つの中で最も良い結果でした。reps=5がreps=4よりわずかに悪い(83.52% < 87.46%)のは、層が増えて探索空間が広がった分だけmaxiter=200という同じ予算では相対的に収束が追いつきにくくなっている、という見方と整合します。
厳密解を引き当てる確率をまとめて棒グラフにすると以下のようになります。
reps=1→2で明確に改善、reps=3でやや後退、reps=4で大きく改善、reps=5でわずかに後退、reps=6でさらに改善、という非単調な挙動が見て取れます。「repsを増やせば増やすほど良くなる」と単純には言えず、repsとmaxiter(古典最適化にかけられる予算)はセットで振る必要がある、というのが移植時に気をつけるべきポイントです。移植したばかりのデフォルト設定(reps=2, maxiter=200程度)を鵜呑みにせず、自分の問題でrepsを振って確率が安定するかを確認したほうがよさそうです。
Ansatz(mixer/initial_state)もデフォルトのままとは限らない
ここまではQAOA(sampler=sampler, optimizer=optimizer, reps=reps)という、mixerやinitial_stateを指定しないデフォルト構成で通しました。この場合qiskit_algorithmsは標準的なQAOAアンザッツ、つまりQUBOから作ったコスト演算子と、全qubitに一様に効くXミキサー(各qubitに$X$回転を掛けるだけの標準ミキサー)を組み合わせた回路を自動的に組み立てます。
ただしインストール済みのqiskit_algorithms.QAOAのソース(qaoa.py)を見ると、コンストラクタはmixer: QuantumCircuit | BaseOperatorとinitial_state: QuantumCircuit | Noneという引数を受け付けており、どちらも明示的に渡せばデフォルトを上書きできることが確認できます。つまりAnsatzもデフォルトのXミキサー一択ではなく、選択肢があります。例えば実際のアニーリング案件では「ちょうど1つだけ1にする」といった制約をQUBO側でペナルティ項として表現することがありますが、そのぶんペナルティの重みという新たなハイパーパラメータが増えてしまいます。QAOA側では、そうした制約を尊重する向きしか遷移させないmixer(いわゆるXYミキサー等)を自作して渡すことで、ペナルティに頼らず制約を回路レベルで表現するというアプローチも取れます。今回の検証では標準のXミキサーのままrepsと古典最適化器だけを振っており、mixer/initial_stateのカスタマイズによる効果までは検証していませんが、移植後にデフォルトの収束が思わしくない場合に触れる余地がある、ということは覚えておいて損はなさそうです。
古典最適化器を振ってみる(COBYLA vs SPSA)
もう一つ、移植時に選択を迫られるのが古典最適化器です。reps=2固定で、COBYLAとSPSA(どちらもmaxiter=200)を比較しました。
| optimizer | 収束(反復)までの評価回数 | 最終<H>
|
実行時間 |
|---|---|---|---|
| COBYLA | 72 | -2.3482 | 0.54秒 |
| SPSA | 451 | -1.7103 | 3.31秒 |
COBYLAは72回の評価で滑らかに収束したのに対し、SPSA(勾配を確率的に近似する手法)はmaxiter=200の指定でも実際には451回コスト関数を評価しており、しかもグラフを見ると全体的にノイズが大きく、200反復消化してもCOBYLAの水準まで下がりきっていません。SPSAは本来ノイズの多い実機環境(コスト関数の評価自体に測定誤差が乗る状況)で威力を発揮する設計の最適化器なので、今回のような理想的なシミュレータ(Sampler()によるノイズなしサンプリング)でCOBYLAに劣るのはある意味自然な結果と言えます。シミュレータでまず動作確認したいだけなら、まずCOBYLAで試すのが無難そうです。
まとめ
- QUBO行列からの移植は
QuadraticProgramへの登録だけで完結する数行の作業です。既に制約なし・バイナリ変数のみのQUBOにQuadraticProgramToQuboのような追加変換を挟むのは無駄な二重変換になります - QAOAは古典最適化ループで変分パラメータを訓練する仕組みなので、移植直後の
result.xをそのまま信用するのは危険です。デフォルト設定(reps=2, COBYLA,maxiter=200)では厳密解を引き当てる確率は57%止まりでしたが、repsを6まで振ると97%まで改善しました。ただし単調ではなく、maxiterとセットで確認する必要があります - 古典最適化器(ノイズなしシミュレータではCOBYLAがSPSAより明確に有利)やAnsatz(
mixer/initial_stateでデフォルトのXミキサー以外も選べる)も、移植後に確認・調整すべきポイントです
注意事項
QPU実機で実行する場合は、ローカルPCとQPU実機の間の通信遅延が反復ごとに生じるため、このスクリプトでは現実的な時間で実行できません。対策として、Qiskit Runtimeを使用する必要があります。
Qiitaでの解説記事もあります。
参考文献
- Qiskit 2.0 チートシート — Qiskit v1.x→v2.0の変更点やprimitives(Sampler/Estimator)のV1→V2移行について
- Qiskit Optimization: Minimum Eigen Optimizer チュートリアル(公式) — QUBOをイジングハミルトニアンに変換し、QAOAやNumPyMinimumEigensolverなどで解く流れの公式チュートリアル
-
QuadraticProgram APIリファレンス(公式) —
binary_varやminimize(linear=, quadratic=)など、今回使ったメソッドの仕様 -
QAOA APIリファレンス(公式) —
reps,optimizer,sampler,callbackなど各パラメータの仕様 -
OpenJij SASampler ソースコード(公式GitHub) —
num_sweeps,num_reads,beta_min/beta_maxなど、SAサンプラーが持つパラメータ -
OpenJij SQASampler ソースコード(公式GitHub) —
beta,gamma,trotterなど、SQAサンプラーが持つパラメータ -
D-Wave QPU Solver Parameters(公式) —
num_reads,annealing_time,chain_strength,anneal_scheduleなど、D-Wave Ocean SDKのQPUパラメータ





