量子アルゴリズム：変分量子アルゴリズム
Takashi Imamichi（2024年5月24日）
元の講義資料はPDFでダウンロードできます。なお、コードスニペットは静的な画像のため、一部が非推奨となっている場合があります。
この実験を実行するおおよそのQPU時間は9分です（Eagleプロセッサーでテスト済み）。
（このノートブックは、Open Planの制限時間内に評価が完了しない場合があります。量子コンピューティングのリソースを大切に使ってください。）
1. はじめに
量子コンピューティングの進歩にもかかわらず、現在の量子デバイスにはノイズが存在するため、深い量子Circuitから有意義な結果を得ることは困難です。この課題を克服するために、VQEとQAOAはハイブリッド量子古典アプローチを採用しています。このアプローチでは、比較的浅い量子Circuitを量子計算で繰り返し実行し、目標とするパラメーター化された量子Circuitのパラメーターを古典計算で最適化します。
QAOAは、さまざまなエラー緩和・抑制技術の適用により、ユーティリティ規模で対象問題の最適解を提供できる可能性があります。VQEは量子化学など多くの応用がありますが、スケーラビリティには課題があります。しかし、Krylov部分空間対角化やサンプリングベースの量子対角化（SQD）など、VQEを補完・拡張する固有値関連のアプローチが数多く登場しています。VQEを理解することは、これまでに登場した幅広い古典量子ハイブリッドアルゴリズムを理解するための重要な第一歩です。
様々な誤り緩和・抑制技術の適用により、QAOAはユーティリティスケールで対象問題に対する最適解を提供できる可能性を秘めています。VQEには（量子化学など）多くの応用がありますが、それらの拡張性は劣ります。しかし、VQEを補完・拡張するために、クリロフ部分空間対角化やサンプリングベースの量子対角化（SQD）など、固有値に関連する様々なアプローチが登場しています。VQEを理解することは、これまで登場してきた幅広い古典・量子ハイブリッド・アルゴリズムを理解するための重要な第一歩です。
本モジュールでは、VQEとQAOAの基本概念と実装について説明します。今後のチュートリアルでは、これらのアルゴリズムをスケールアップするための高度なトピックと技術を解説します。 このノートブックを実行するには、以下のライブラリが環境にインストールされている必要があります。 まだインストールしていない場合は、以下のセルのコメントを外して実行することでインストールできます。
# Added by doQumentation — installs the packages this notebook needs if they are missing
import importlib.util
_needed = {"matplotlib": "matplotlib", "numpy": "numpy", "qiskit": "qiskit", "qiskit_ibm_runtime": "qiskit-ibm-runtime", "rustworkx": "rustworkx", "scipy": "scipy"}
_missing = [pip for module, pip in _needed.items()
if importlib.util.find_spec(module) is None]
# One at a time, so a package that fails to install does not block the others
for _pip in _missing:
%pip install -q {_pip}
if not _missing:
print("\u2713 All required packages are installed")
# % pip install 'qiskit[visualization]' qiskit-ibm-runtime
2. 単純なハミルトニアンの最小固有値の計算
まず、VQEがどのように機能するかを理解するために、非常に単純なケースに適用してみます。VQEを使ってパウリ 行列の最小固有値を計算します。最初に、いくつかの一般的なパッケージをインポートします。
import numpy as np
from qiskit.circuit import ParameterVector, QuantumCircuit
from qiskit.primitives import StatevectorEstimator, StatevectorSampler
from qiskit.quantum_info import SparsePauliOp
from scipy.optimize import minimize
次に、対象の演算子を定義し、行列形式で確認します。
op = SparsePauliOp("Z")
op.to_matrix()
array([[ 1.+0.j, 0.+0.j],
[ 0.+0.j, -1.+0.j]])
固有値は古典的な方法で容易に求めることができるので、結果を検証できます。ユーティリティ規模にスケールアップするにつれて、これは困難になる可能性があります。ここではnumpyを使用します。
# compute eigenvalues with numpy
result = np.linalg.eigh(op.to_matrix())
print("Eigenvalues:", result.eigenvalues)
Eigenvalues: [-1. 1.]
変分量子アルゴリズムを使って固有値を求めるために、変分パラメーターを持つGateで構成されたCircuitを構築します。
# define a variational form
param = ParameterVector("a", 3)
qc = QuantumCircuit(1, 1)
qc.u(param[0], param[1], param[2], 0)
qc_estimator = qc.copy()
qc.measure(0, 0)
qc.draw("mpl")
演算子（ など）の期待値を推定したい場合はEstimatorを使用します。システムの状態を確認したい場合はSamplerを使用します。
sampler = StatevectorSampler()
estimator = StatevectorEstimator()
Samplerを使って、ランダムなパラメーター値 [1, 2, 3] でビット列0と1のカウントを計算できます。
# compute counts of bitstrings with random parameter values by Sampler
result = sampler.run([(qc, [1, 2, 3])]).result()
counts = result[0].data.c.get_counts()
counts
{'0': 783, '1': 241}
確率 を用いて によりZの期待値を計算できます。
# compute the expectation value of Z based on the counts
(counts.get("0", 0) - counts.get("1", 0)) / sum(counts.values())
0.529296875
このCircuitは機能しましたが、選択したパラメーター値は非常に低エネルギー（または低固有値）の状態に対応していませんでした。得られた固有値は最小値よりもかなり高い値です。Estimatorを使用した場合も同様の結果が得られます。
EstimatorはBitの測定を含まない量子Circuitを受け取ることに注意してください。
result = estimator.run([(qc_estimator, op, [1, 2, 3])]).result()
result[0].data.evs
array(0.54030231)
パラメーターを探索して、最低の固有値をもたらすものを見つける必要があります。 変分形式のパラメーター値を受け取り、期待値 を返す関数を作成します。
# define a cost function to look for the minimum eigenvalue of Z
def cost(x):
result = sampler.run([(qc, x)]).result()
counts = result[0].data.c.get_counts()
expval = (counts.get("0", 0) - counts.get("1", 0)) / sum(counts.values())
# the following line shows the trajectory of the optimization
print(expval, counts)
return expval
SciPyの minimize 関数を適用して、Zの最小固有値を求めましょう。
# minimize the cost function with scipy's minimize
min_result = minimize(cost, [0, 0, 0], method="COBYLA", tol=1e-8)
min_result
1.0 {'0': 1024}
0.494140625 {'0': 765, '1': 259}
0.466796875 {'0': 751, '1': 273}
0.564453125 {'0': 801, '1': 223}
-0.4296875 {'1': 732, '0': 292}
-0.984375 {'1': 1016, '0': 8}
-0.8984375 {'1': 972, '0': 52}
-0.990234375 {'1': 1019, '0': 5}
-0.892578125 {'1': 969, '0': 55}
-0.986328125 {'1': 1017, '0': 7}
-0.861328125 {'1': 953, '0': 71}
-1.0 {'1': 1024}
-0.982421875 {'1': 1015, '0': 9}
-0.99609375 {'1': 1022, '0': 2}
-0.986328125 {'1': 1017, '0': 7}
-1.0 {'1': 1024}
-0.990234375 {'1': 1019, '0': 5}
-0.998046875 {'1': 1023, '0': 1}
-0.99609375 {'1': 1022, '0': 2}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-0.998046875 {'1': 1023, '0': 1}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-0.998046875 {'1': 1023, '0': 1}
-0.998046875 {'1': 1023, '0': 1}
-0.998046875 {'1': 1023, '0': 1}
-1.0 {'1': 1024}
-0.99609375 {'1': 1022, '0': 2}
-1.0 {'1': 1024}
-0.99609375 {'1': 1022, '0': 2}
-0.998046875 {'1': 1023, '0': 1}
-0.998046875 {'1': 1023, '0': 1}
-0.99609375 {'1': 1022, '0': 2}
-0.998046875 {'1': 1023, '0': 1}
-1.0 {'1': 1024}
-0.998046875 {'1': 1023, '0': 1}
-0.998046875 {'1': 1023, '0': 1}
-0.99609375 {'1': 1022, '0': 2}
-1.0 {'1': 1024}
-0.998046875 {'1': 1023, '0': 1}
-1.0 {'1': 1024}
-0.998046875 {'1': 1023, '0': 1}
-0.998046875 {'1': 1023, '0': 1}
-1.0 {'1': 1024}
-0.998046875 {'1': 1023, '0': 1}
-0.998046875 {'1': 1023, '0': 1}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-0.998046875 {'1': 1023, '0': 1}
-0.994140625 {'1': 1021, '0': 3}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
message: Optimization terminated successfully.
success: True
status: 1
fun: -1.0
x: [ 3.182e+00 1.338e+00 1.664e-01]
nfev: 63
maxcv: 0.0
# check counts of bitstrings with the optimal parameters
result = sampler.run([(qc, min_result.x)]).result()
result[0].data.c.get_counts()
{'0': 1, '1': 1023}
2.1 演習
VQEを使って の最小固有値を計算してください。
z2 = SparsePauliOp("ZZ")
print(z2)
print(z2.to_matrix())
SparsePauliOp(['ZZ'],
coeffs=[1.+0.j])
[[ 1.+0.j 0.+0.j 0.+0.j 0.+0.j]
[ 0.+0.j -1.+0.j 0.+0.j 0.+0.j]
[ 0.+0.j 0.+0.j -1.+0.j 0.+0.j]
[ 0.+0.j 0.+0.j 0.+0.j 1.+0.j]]
# compute eigenvalues with numpy
# define a variational form
# qc = ...
# compute counts of bitstrings with a random parameter values by Sampler
# result = sampler.run(...)
# result
# compute the expectation value of ZZ based on the counts
# verify the expectation value of ZZ with Estimator
# define a cost function to look for the minimum eigenvalue of ZZ
# def cost(x):
# expval = ...
# return expval
# minimize the cost function with scipy's minimize
# min_result = minimize(cost, [...], method="COBYLA", tol=1e-8)
# min_result
# check counts of bitstrings with the optimal parameter values
# result = sampler.run(qc, min_result.x).result()
# result
演習の解答
対象の演算子を定義し、行列形式で確認します。
z2 = SparsePauliOp("ZZ")
print(z2)
print(z2.to_matrix())
SparsePauliOp(['ZZ'],
coeffs=[1.+0.j])
[[ 1.+0.j 0.+0.j 0.+0.j 0.+0.j]
[ 0.+0.j -1.+0.j 0.+0.j 0.+0.j]
[ 0.+0.j 0.+0.j -1.+0.j 0.+0.j]
[ 0.+0.j 0.+0.j 0.+0.j 1.+0.j]]
変分量子アルゴリズムを使って固有値を求めるために、変分パラメーターを持つGateで構成されたCircuitを構築します。
# define a variational form
param = ParameterVector("a", 6)
qc = QuantumCircuit(2, 2)
qc.u(param[0], param[1], param[2], 0)
qc.u(param[3], param[4], param[5], 1)
qc_estimator = qc.copy()
qc.measure([0, 1], [0, 1])
qc.draw("mpl")
演算子（ など）の期待値を推定したい場合はEstimatorを使用します。システムの状態を確認したい場合はSamplerを使用します。
sampler = StatevectorSampler()
estimator = StatevectorEstimator()
# compute counts of bitstrings with random parameter values by Sampler
result = sampler.run([(qc, [1, 2, 3, 4, 5, 6])]).result()
counts = result[0].data.c.get_counts()
counts
{'10': 661, '11': 203, '01': 47, '00': 113}
# compute the expectation value of ZZ based on the counts
(
counts.get("00", 0)
- counts.get("01", 0)
- counts.get("10", 0)
+ counts.get("11", 0)
) / sum(counts.values())
-0.3828125
このCircuitは機能しましたが、選択したパラメーター値は非常に低エネルギー（または低固有値）の状態に対応していませんでした。得られた固有値は最小値よりもかなり高い値です。Estimatorを使用した場合も同様の結果が得られます。
# verify the expectation value of ZZ with Estimator
result = estimator.run([(qc_estimator, z2, [1, 2, 3, 4, 5, 6])]).result()
result[0].data.evs
array(-0.35316516)
パラメーターを探索して、最低の固有値をもたらすものを見つける必要があります。
# define a cost function to look for the minimum eigenvalue of ZZ
def cost(x):
result = sampler.run([(qc, x)]).result()
counts = result[0].data.c.get_counts()
expval = (
counts.get("00", 0)
- counts.get("01", 0)
- counts.get("10", 0)
+ counts.get("11", 0)
) / sum(counts.values())
print(expval, counts)
return expval
# minimize the cost function with scipy's minimize
min_result = minimize(cost, [0, 0, 0, 0, 0, 0], method="COBYLA", tol=1e-8)
min_result
1.0 {'00': 1024}
0.578125 {'00': 808, '01': 216}
0.5234375 {'00': 780, '01': 244}
0.548828125 {'00': 793, '01': 231}
0.3515625 {'00': 637, '10': 164, '11': 55, '01': 168}
0.3359375 {'00': 638, '11': 46, '10': 174, '01': 166}
0.283203125 {'00': 602, '10': 181, '01': 186, '11': 55}
-0.087890625 {'01': 414, '00': 184, '10': 143, '11': 283}
0.236328125 {'10': 27, '11': 623, '01': 364, '00': 10}
-0.0625 {'11': 261, '01': 403, '00': 219, '10': 141}
0.248046875 {'01': 366, '11': 628, '00': 11, '10': 19}
-0.0625 {'10': 145, '11': 254, '01': 399, '00': 226}
0.228515625 {'01': 373, '11': 609, '00': 20, '10': 22}
0.0546875 {'11': 376, '10': 273, '01': 211, '00': 164}
-0.447265625 {'01': 731, '10': 10, '11': 267, '00': 16}
-0.71484375 {'01': 871, '11': 99, '00': 47, '10': 7}
-0.46484375 {'01': 741, '00': 253, '10': 9, '11': 21}
-0.87890625 {'01': 962, '00': 39, '11': 23}
-0.640625 {'00': 176, '01': 837, '11': 8, '10': 3}
-0.88671875 {'01': 966, '00': 41, '11': 17}
-0.994140625 {'01': 1021, '11': 3}
-0.91796875 {'01': 982, '11': 35, '00': 7}
-0.994140625 {'01': 1021, '11': 2, '00': 1}
-0.939453125 {'01': 993, '00': 31}
-0.990234375 {'01': 1019, '11': 5}
-0.90234375 {'01': 974, '00': 21, '11': 29}
-0.98046875 {'01': 1014, '11': 10}
-0.994140625 {'01': 1021, '00': 3}
-0.990234375 {'01': 1019, '11': 4, '00': 1}
-0.98828125 {'01': 1018, '11': 6}
-0.990234375 {'01': 1019, '11': 4, '00': 1}
-0.994140625 {'01': 1021, '11': 2, '00': 1}
-0.99609375 {'01': 1022, '11': 2}
-0.998046875 {'01': 1023, '00': 1}
-0.99609375 {'01': 1022, '00': 2}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '00': 1}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '00': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '00': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '00': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '00': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-0.99609375 {'01': 1022, '00': 1, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '00': 1}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-0.99609375 {'01': 1022, '11': 1, '00': 1}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '00': 1}
-0.994140625 {'01': 1021, '00': 3}
-0.998046875 {'01': 1023, '00': 1}
-0.99609375 {'01': 1022, '11': 2}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '00': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-0.99609375 {'01': 1022, '11': 2}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
message: Optimization terminated successfully.
success: True
status: 1
fun: -0.998046875
x: [ 3.167e+00 6.940e-01 1.033e+00 -2.894e-02 8.933e-01
1.885e+00]
nfev: 128
maxcv: 0.0
message: Optimization terminated successfully.
success: True
status: 1
fun: -0.99609375
x: [ 3.098e+00 -5.402e-01 1.091e+00 -1.004e-02 3.615e-01
6.913e-01]
nfev: 115
maxcv: 0.0
numpyから得られた最小値に非常に近い固有値が得られました。
# check counts of bitstrings with the optimal parameters
result = sampler.run([(qc, min_result.x)]).result()
result[0].data.c.get_counts()
{'01': 1024}
3. Qiskitパターンによる量子最適化
このハウツーでは、Qiskitパターンと量子近似最適化について学びます。Qiskitパターンとは、量子コンピューティングのワークフローを実装するための直感的で再現可能な一連の手順です。
このパターンを組合せ最適化のコンテキストに適用し、ハイブリッド（量子古典）反復手法である**量子近似最適化アルゴリズム（QAOA）を使って最大カット（Max-Cut）**問題を解く方法を示します。
このQAOA部分は、量子近似最適化アルゴリズムチュートリアルの「パート1：小規模QAOA」に基づいています。スケールアップの方法については、そのチュートリアルを参照してください。
3.1 （小規模）最適化問題のQiskitパターン
このパートでは、小規模なマックスカット問題を用いて、量子コンピューターを使って最適化問題を解くために必要な手順を説明します。
マックスカット問題は解くのが難しい（より具体的にはNP困難問題である）最適化問題であり、クラスタリング、ネットワーク科学、統計物理学など様々な応用があります。このチュートリアルでは、エッジで接続されたノードから成るグラフを考え、カットされるエッジの数が最大になるように、ノードを2つの集合に「カット」することで分割することを目指します。
この問題を量子アルゴリズムにマッピングする前に少し背景を説明すると、まず関数の最小化を考えることで、マックスカット問題がどのように古典的な組合せ最適化問題になるかをよりよく理解できます。
ここで、入力はグラフの各ノードに対応する成分を持つベクトルです。そして、これらの各成分をまたは（それぞれカットに含まれるか含まれないかを表す）に制約します。この小規模な例では、個のノードを持つグラフを使用します。
ノードのペアの関数を書いて、対応するエッジがカットに含まれるかどうかを示すことができます。例えば、関数は、またはのどちらか一方だけが1の場合（つまりエッジがカットに含まれる場合）にのみ1となり、それ以外の場合は0になります。カットに含まれるエッジを最大化する問題は、次のように定式化できます。
これは次のような最小化の形に書き換えることができます。
この場合、の最小値は、カットによって横断されるエッジの数が最大となるときです。ご覧のとおり、ここまでは量子コンピューティングに関連するものは何もありません。この問題を、量子コンピューターが理解できる形に再定式化する必要があります。 個のノードを持つグラフを作成して、問題を初期化しましょう。
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import rustworkx as rx
from rustworkx.visualization import mpl_draw
n = 5
graph = rx.PyGraph()
graph.add_nodes_from(range(1, n + 1))
edge_list = [
(0, 1, 1.0),
(0, 2, 1.0),
(1, 2, 1.0),
(1, 3, 1.0),
(2, 4, 1.0),
(3, 4, 1.0),
]
graph.add_edges_from(edge_list)
pos = rx.spring_layout(graph, seed=2)
mpl_draw(graph, node_size=600, pos=pos, with_labels=True, labels=str)
3.2 ステップ1：古典的入力を量子問題にマッピングする
パターンの最初のステップは、古典的な問題（グラフ）を量子 Circuit および 演算子 にマッピングすることです。これには主に3つの手順があります。
-
一連の数学的再定式化を利用して、この問題を二次制約なし二値最適化（QUBO）問題の表記法で表現します。
-
最適化問題をハミルトニアンとして書き直します。このハミルトニアンの基底状態がコスト関数を最小化する解に対応します。
-
量子アニーリングに類似した処理によって、このハミルトニアンの基底状態を準備する量子 Circuit を作成します。
注意： QAOAの手法では、最終的にハイブリッドアルゴリズムのコスト関数を表す演算子（ハミルトニアン）と、問題の候補解を表す量子状態を表すパラメータ化された Circuit（Ansatz）が必要です。これらの候補状態からサンプリングし、コスト関数を使って評価することができます。
グラフ → 最適化問題
マッピングの最初のステップは表記法の変換です。以下はQUBO表記法で問題を表現したものです。
ここで は の実数行列、 はグラフのノード数、 は上で導入した二値変数のベクトル、 はベクトル の転置を表します。
Problem name: maxcut
Minimize
2*x_1*x_2 + 2*x_1*x_3 + 2*x_2*x_3 + 2*x_2*x_4 + 2*x_3*x_5 + 2*x_4*x_5 - 2*x_1
- 3*x_2 - 3*x_3 - 2*x_4 - 2*x_5
Subject to
No constraints
Binary variables (5)
x_1 x_2 x_3 x_4 x_5
最適化問題 → ハミルトニアン
続いてQUBO問題をハミルトニアン（ここではシステムのエネルギーを表す行列）として再定式化できます。
QAOA問題からハミルトニアンへの再定式化手順
QAOA問題をこのように書き直せることを示すために、まず二値変数 を新しい変数 に次の変換で置き換えます。
ここで、 が なら は でなければならないことがわかります。最適化問題（）内の を で置き換えると、同等の定式化が得られます。
ここで と定義し、係数と定数項 を除去すると、同一の最適化問題の2つの同値な定式化に到達します。
ここで は に依存します。 を得るために、最適化において役割を果たさない の係数と定数オフセット を除去したことに注意してください。
次に、問題の量子定式化を得るために、 変数をパウリ 行列（次の形の 行列）に昇格させます。
これらの行列を上の最適化問題に代入すると、以下のハミルトニアンが得られます。
また、 行列は量子コンピュータの計算空間（つまり のサイズのヒルベルト空間）に埋め込まれることを思い出してください。したがって、 などの項は のヒルベルト空間に埋め込まれたテンソル積 として理解すべきです。たとえば、5つの決定変数を持つ問題では、 という項は （ は 単位行列）を意味します。
このハミルトニアンはコスト関数ハミルトニアンと呼ばれます。その基底状態がコスト関数 を最小化する解に対応するという性質を持ちます。 したがって、最適化問題を解くには、量子コンピュータ上で の基底状態（あるいは高い重複度を持つ状態）を準備する必要があります。この状態からサンプリングすることで、高い確率で の解が得られます。
def build_max_cut_operator(graph: rx.PyGraph) -> tuple[SparsePauliOp, float]:
sp_list = []
constant = 0
for s, t in graph.edge_list():
w = graph.get_edge_data(s, t)
sp_list.append(("ZZ", [s, t], w / 2))
constant -= 1 / 2
return SparsePauliOp.from_sparse_list(
sp_list, num_qubits=graph.num_nodes()
), constant
cost_hamiltonian, constant = build_max_cut_operator(graph)
print("Cost Function Hamiltonian:", cost_hamiltonian)
print("Constant:", constant)
Cost Function Hamiltonian: SparsePauliOp(['IIIZZ', 'IIZIZ', 'IIZZI', 'IZIZI', 'ZIZII', 'ZZIII'],
coeffs=[0.5+0.j, 0.5+0.j, 0.5+0.j, 0.5+0.j, 0.5+0.j, 0.5+0.j])
Constant: -3.0
ハミルトニアン → 量子Circuit
ハミルトニアン には問題の量子的定義が含まれています。次に、量子コンピュータから良い解をサンプリングするための量子 Circuit を作成できます。QAOAは量子アニーリングからインスピレーションを受けており、量子 Circuit 内で交互に演算子の層を適用します。
一般的なアイデアは、既知のシステムの基底状態 から始め、目的のコスト演算子の基底状態へシステムを導くことです。これは、角度 および を持つ演算子 と を交互に適用することで実現されます。
生成する量子回路はとによってパラメータ化されているため、との異なる値を試し、得られた状態からサンプリングすることができます。
この場合、2つのパラメータとを含む1つのQAOA層を使用する例を試してみます。
from qiskit.circuit.library import QAOAAnsatz
circuit = QAOAAnsatz(cost_operator=cost_hamiltonian, reps=1)
circuit.measure_all()
circuit.draw("mpl")
circuit.decompose(reps=3).draw("mpl", fold=-1)
circuit.parameters
ParameterView([ParameterVectorElement(β[0]), ParameterVectorElement(γ[0])])
3.3 ステップ2：量子ハードウェア実行に向けてCircuitを最適化する
上記の Circuit には量子アルゴリズムを考える上で有用な抽象化が含まれていますが、ハードウェア上ではそのままでは実行できません。QPU上で実行するためには、パターンのトランスパイル（Circuit最適化）ステップを構成する一連の操作を Circuit に施す必要があります。
Qiskit ライブラリには、幅広い Circuit 変換に対応する一連のトランスパイルパスが用意されています。目的に応じて Circuit が最適化されていることを確認する必要があります。
トランスパイルには以下のようないくつかのステップが含まれます。
-
Circuit 内の qubit（決定変数など）をデバイス上の物理 qubit に初期マッピングする。
-
量子 Circuit 内の命令を、Backend が理解できるハードウェアネイティブな命令に展開（アンロール）する。
-
Circuit 内で相互作用する qubit を、互いに隣接した物理 qubit にルーティングする。
-
ダイナミカルデカップリングによってノイズを抑制する単一 qubit Gate を追加し、エラー抑制を行う。
トランスパイルの詳細についてはドキュメントを参照してください。
以下のコードは、抽象的な Circuit を変換・最適化して、Qiskit IBM® Runtime サービスを通じてクラウド経由でアクセス可能なデバイスで実行できる形式に変換します。
なお、実際の量子コンピュータに送信する前に「ローカルテストモード」を使ってプログラムをローカルでテストすることができます。 ローカルテストモードの詳細についてはドキュメントを参照してください。
from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
# Use a quantum device
service = QiskitRuntimeService()
backend = service.least_busy(min_num_qubits=127)
# backend = service.backend("ibm_kingston")
# You can test your programs locally with a fake backend (local testing mode)
# backend = FakeBrisbane()
print(backend)
# Create pass manager for transpilation
pm = generate_preset_pass_manager(optimization_level=3, backend=backend)
candidate_circuit = pm.run(circuit)
candidate_circuit.draw("mpl", fold=False, idle_wires=False)
service = QiskitRuntimeService(channel="ibm_quantum_platform")
<IBMBackend('ibm_strasbourg')>
3.4 ステップ3：IBM Quantumプリミティブを使って実行する
QAOAのワークフローでは、最適なQAOAパラメータを反復的な最適化ループで探索します。このループでは一連の Circuit 評価を実行し、古典的な最適化器を使って最適な および パラメータを見つけます。この実行ループは以下の手順で実行されます。
-
初期パラメータを定義する
-
最適化ループと Circuit のサンプリングに使用するプリミティブを含む新しい
Sessionをインスタンス化する -
最適なパラメータセットが見つかったら、最終的な分布を得るために Circuit を最後に1回実行し、後処理ステップで使用する
初期パラメータを持つCircuitの定義
任意に選んだパラメータから開始します。
initial_gamma = np.pi
initial_beta = np.pi / 2
init_params = [initial_gamma, initial_beta]
Backend と実行プリミティブの定義
IBM Quantum プリミティブを使って IBM® Backend と通信します。2つのプリミティブは Sampler と Estimator であり、量子コンピューター上で実行したい測定の種類によって選択します。 の最小化には Estimator を使用します。コスト関数の測定は単に の期待値を取ることだからです。
実行
プリミティブは、量子デバイス上でワークロードをスケジュールするための様々な実行モードを提供しており、QAOAのワークフローはセッション内で反復的に実行されます。
Samplerベースのコスト関数をSciPyの最小化ルーチンに組み込んで、最適なパラメータを見つけることができます。
def cost_func_estimator(params, ansatz, hamiltonian, estimator):
# transform the observable defined on virtual qubits to
# an observable defined on all physical qubits
isa_hamiltonian = hamiltonian.apply_layout(ansatz.layout)
pub = (ansatz, isa_hamiltonian, params)
job = estimator.run([pub])
results = job.result()[0]
cost = results.data.evs
objective_func_vals.append(cost)
return cost
from qiskit_ibm_runtime import Session, EstimatorV2
from scipy.optimize import minimize
objective_func_vals = [] # Global variable
with Session(backend=backend) as session:
# If using qiskit-ibm-runtime<0.24.0, change `mode=` to `session=`
estimator = EstimatorV2(mode=session)
estimator.options.default_shots = 1000
# Set simple error suppression/mitigation options
estimator.options.dynamical_decoupling.enable = True
estimator.options.dynamical_decoupling.sequence_type = "XY4"
estimator.options.twirling.enable_gates = True
estimator.options.twirling.num_randomizations = "auto"
result = minimize(
cost_func_estimator,
init_params,
args=(candidate_circuit, cost_hamiltonian, estimator),
method="COBYLA",
tol=1e-2,
)
print(result)
message: Optimization terminated successfully.
success: True
status: 1
fun: -0.6557925874481715
x: [ 2.873e+00 9.414e-01]
nfev: 21
maxcv: 0.0
最適化器は、コストを削減し、回路のより良いパラメータを見つけることができました。
plt.figure(figsize=(12, 6))
plt.plot(objective_func_vals)
plt.xlabel("Iteration")
plt.ylabel("Cost")
plt.show()
回路の最適なパラメータが見つかったら、これらのパラメータを割り当てて、最適化されたパラメータで得られた最終的な分布をサンプリングできます。ここでSamplerプリミティブを使用する必要があります。なぜなら、それはグラフの最適なカットに対応するビット列測定の確率分布だからです。
注： これは、コンピューター内で量子状態を準備し、それを測定することを意味します。測定によって状態は単一の計算基底状態に収縮します。例えば、010101110000...は、初期の最適化問題（タスクに応じてまたは）に対する候補解に対応します。
optimized_circuit = candidate_circuit.assign_parameters(result.x)
optimized_circuit.draw("mpl", fold=False, idle_wires=False)
from qiskit_ibm_runtime import SamplerV2
# If using qiskit-ibm-runtime<0.24.0, change `mode=` to `backend=`
sampler = SamplerV2(mode=backend)
# Set simple error suppression/mitigation options
sampler.options.dynamical_decoupling.enable = True
sampler.options.dynamical_decoupling.sequence_type = "XY4"
sampler.options.twirling.enable_gates = True
sampler.options.twirling.num_randomizations = "auto"
pub = (optimized_circuit,)
job = sampler.run([pub], shots=int(1e4))
counts_int = job.result()[0].data.meas.get_int_counts()
counts_bin = job.result()[0].data.meas.get_counts()
shots = sum(counts_int.values())
final_distribution_int = {key: val / shots for key, val in counts_int.items()}
final_distribution_bin = {key: val / shots for key, val in counts_bin.items()}
print(final_distribution_int)
{12: 0.0652, 31: 0.0089, 4: 0.0085, 13: 0.0731, 26: 0.0256, 28: 0.0246, 17: 0.0405, 25: 0.0591, 20: 0.031, 15: 0.0221, 8: 0.017, 21: 0.0371, 14: 0.0461, 16: 0.0229, 19: 0.0723, 23: 0.0199, 22: 0.0478, 18: 0.0708, 24: 0.0165, 6: 0.0525, 7: 0.0155, 5: 0.0245, 3: 0.0231, 29: 0.0121, 30: 0.0062, 10: 0.0363, 1: 0.0097, 9: 0.042, 27: 0.0094, 11: 0.0349, 0: 0.0129, 2: 0.0119}
3.5 ステップ4：後処理と古典形式での結果の返却
後処理のステップでは、サンプリング出力を解釈して元の問題の解を返します。今回は、最適なカットを決定するために最も確率の高いビット列に注目します。問題の対称性により4つの解が存在しえますが、サンプリングによってそのうちの1つがわずかに高い確率で返されます。ただし、以下に示す分布のプロットを見ると、4つのビット列が他と比べて明らかに高い確率を持っていることが分かります。
# auxiliary functions to sample most likely bitstring
def to_bitstring(integer, num_bits):
result = np.binary_repr(integer, width=num_bits)
return [int(digit) for digit in result]
keys = list(final_distribution_int.keys())
values = list(final_distribution_int.values())
most_likely = keys[np.argmax(np.abs(values))]
most_likely_bitstring = to_bitstring(most_likely, len(graph))
most_likely_bitstring.reverse()
print("Result bitstring:", most_likely_bitstring)
Result bitstring: [1, 0, 1, 1, 0]
import matplotlib.pyplot as plt
matplotlib.rcParams.update({"font.size": 10})
final_bits = final_distribution_bin
values = np.abs(list(final_bits.values()))
top_4_values = sorted(values, reverse=True)[:4]
positions = []
for value in top_4_values:
positions.append(np.where(values == value)[0])
fig = plt.figure(figsize=(11, 6))
ax = fig.add_subplot(1, 1, 1)
plt.xticks(rotation=45)
plt.title("Result Distribution")
plt.xlabel("Bitstrings (reversed)")
plt.ylabel("Probability")
ax.bar(list(final_bits.keys()), list(final_bits.values()), color="tab:grey")
for p in positions:
ax.get_children()[p[0].item()].set_color("tab:purple")
plt.show()
最適カットの可視化
最適なビット列から、元のグラフ上でそのカットを可視化できます。
colors = ["tab:grey" if i == 0 else "tab:purple" for i in most_likely_bitstring]
mpl_draw(graph, node_size=600, pos=pos, with_labels=True, labels=str, node_color=colors)
そして、カットの値を計算します。ノイズの影響により解は最適ではありませんが（最適解のカット値は5です）、結果を確認してみましょう。
from typing import Sequence
def evaluate_sample(x: Sequence[int], graph: rx.PyGraph) -> float:
assert len(x) == len(
list(graph.nodes())
), "The length of x must coincide with the number of nodes in the graph."
return sum(
x[u] * (1 - x[v]) + x[v] * (1 - x[u]) for u, v in list(graph.edge_list())
)
cut_value = evaluate_sample(most_likely_bitstring, graph)
print("The value of the cut is:", cut_value)
The value of the cut is: 5
以上で、小規模なQAOAチュートリアルは完了です。 ユーティリティスケールでのQAOAの適用方法については、量子近似最適化アルゴリズムチュートリアルの「Part 2: スケールアップ！」で学ぶことができます。
# Check Qiskit version
import qiskit
qiskit.__version__
'2.0.2'