量子アルゴリズム：変分量子アルゴリズム

Note

今道隆司（2024年5月24日）

講演の原文の PDFをダウンロードする。これらは静的画像なので、いくつかのコード・スニペットは非推奨になるかもしれないことに注意してください。

この実験の実行にかかるQPU時間の目安は9分（Eagleプロセッサーでテスト）。

(このノートはオープン・プランの時間内に評価できないかもしれない）。量子コンピューティングのリソースを賢く使ってください)

1. はじめに

このチュートリアルでは、特に変分量子固有値ソルバー（VQE）と量子近似最適化アルゴリズム（QAOA）に焦点を当て、ハイブリッド量子古典アルゴリズムの概要を説明します。これらのアルゴリズムの主な目的は、パラメータ化された量子ゲートを持つ量子回路を用いて最適化問題に取り組むことである。

量子コンピューティングの進歩にもかかわらず、現在の量子デバイスにはノイズが存在するため、深層量子回路から意味のある結果を引き出すことは困難である。この課題を克服するため、VQEとQAOAは、量子計算で比較的短い量子回路を繰り返し実行し、古典計算で目標とするパラメトリック量子回路のパラメータを最適化するという、ハイブリッド量子古典アプローチを採用している。

QAOAは、様々なエラー緩和・抑制技術の適用により、ユーティリティ・スケールで目標問題に対する最適解を提供する可能性を持っている。 VQEには（量子化学のように）スケーラビリティの低いアプリケーションも多い。しかし、VQEを補完・補強するために、クリロフ部分空間対角化やサンプリングに基づく量子対角化（SQD）など、固有値に関連するアプローチが数多く登場している。 VQEを理解することは、登場した幅広い古典-量子ハイブリッドアルゴリズムを理解するための重要な第一歩である。

このモジュールでは、VQEとQAOAの基本的なコンセプトと実装について説明する。さらなるチュートリアルでは、これらのアルゴリズムをスケールアップするための高度なトピックとテクニックを探求する。

このノートブックを実行するには、以下のライブラリーが必要です。まだインストールしていない場合は、コメントを外して以下のセルを実行すればインストールできる。

# % pip install 'qiskit[visualization]' qiskit-ibm-runtime

2. 単純なハミルトニアンの最小固有値の計算

まずは、VQEがどのように機能するのかを確認するために、非常にシンプルなケースにVQEを適用してみる。 VQEを用いてパウリ $Z$ 行列の最小固有値を計算する。まずは一般的なパッケージをいくつかインポートする。

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()

Output:

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)

Output:

Eigenvalues: [-1.  1.]

変分量子アルゴリズムを使って固有値を求めるために、変分パラメータを取るゲートで回路を構成する：

# 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")

Output:

ある演算子（ $Z$ など）の期待値を推定したい場合は、Estimator を使用する。システムの状態を見たい場合は、サンプラーを使う。

sampler = StatevectorSampler()
estimator = StatevectorEstimator()

Samplerを使って、ビット列0と1のカウントをランダムなパラメータ値 [1, 2, 3] 。

# 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

Output:

{'0': 783, '1': 241}

Z の期待値は、確率 $\{0: p_0, 1: p_1\}$ を用いて $\langle Z \rangle = p_0 - p_1$ で計算できることがわかっている。

# compute the expectation value of Z based on the counts
(counts.get("0", 0) - counts.get("1", 0)) / sum(counts.values())

Output:

0.529296875

この回路は機能したが、選ばれたパラメータ値は、非常に低エネルギー（あるいは低固有値）状態には対応しなかった。得られた固有値は、最小値よりもかなり高い。 estimatorを使っても結果は同様である。

なお、エスティメーターは測定なしで量子回路を作る。

result = estimator.run([(qc_estimator, op, [1, 2, 3])]).result()
result[0].data.evs

Output:

array(0.54030231)

パラメータを検索し、最小固有値をもたらすものを見つける必要がある。変分形式のパラメータ値を受け取り、期待値 $\langle Z \rangle$ を返す関数を作る。

# 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's minimize 関数を適用して、Zの固有値の最小値を求めてみよう。

# minimize the cost function with scipy's minimize
min_result = minimize(cost, [0, 0, 0], method="COBYLA", tol=1e-8)
min_result

Output:

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()

Output:

{'0': 1, '1': 1023}

2.1 運動

$Z \otimes Z$ の最小固有値を VQE で計算する。

z2 = SparsePauliOp("ZZ")
print(z2)
print(z2.to_matrix())

Output:

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())

Output:

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]]

変分量子アルゴリズムを使って固有値を求めるために、変分パラメータを取るゲートで回路を構成する：

# 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")

Output:

ある演算子（ $Z \otimes Z$ など）の期待値を推定したい場合は、Estimator を使うことになる。システムの状態を見たい場合は、サンプラーを使う。

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

Output:

{'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())

Output:

-0.3828125

# 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

Output:

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

Output:

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()

Output:

{'01': 1024}

3. Qiskitパターンを用いた量子最適化

このハウツーでは、Qiskitパターンと量子近似最適化について学ぶ。 Qiskitパターンは、量子コンピューティングのワークフローを実装するための、直感的で反復可能なステップのセットです：

これらのパターンを組み合わせ最適化の文脈に適用し、ハイブリッド（量子・古典）反復法である量子近似最適化アルゴリズム（QAOA） を用いて最大切断問題を解く方法を示す。

なお、このQAOAパートは、量子近似最適化アルゴリズム・チュートリアルの「パート1：小規模QAOA」に基づいています。拡大する方法はチュートリアルをご覧ください。

3.1 (小規模) 最適化のためのQiskitパターン

このセクションでは、小規模な最大切断問題を用いて、量子コンピュータを用いて最適化問題を解くために必要な手順を説明します。

最大カット問題は、クラスタリング、ネットワーク科学、統計物理学など、さまざまな分野で応用されているが、解くのが困難な最適化問題（より具体的には、NP困難問題）である。このチュートリアルでは、辺で結ばれたノードからなるグラフを対象とし、辺を「切断」することでノードを2つの集合に分割し、切断される辺の数を最大化する方法を解説します。

この問題を量子アルゴリズムに適用する前に、背景を説明しておくと、まず関数 $f(x)$ の最小化について考えることで、最大切断問題がどのように古典的な組み合わせ最適化問題となるのかをよりよく理解できるでしょう

\min_{x\in \{0, 1\}^n}f(x),

ここで入力 $x$ は、グラフの各ノードに対応する成分を持つベクトルである。次に、これらの各成分を $0$ または $1$ （カットに含まれるか含まれないかを表す）のいずれかになるように制約する。この小規模な例では、 $n=5$ のノードを持つグラフを使用する。

ノードのペア $i,j$ 、対応するエッジ $(i,j)$ がカット内にあるかどうかを示す関数を書くことができる。例えば、関数 $x_i + x_j - 2 x_i x_j$ は、 $x_i$ と $x_j$ のどちらかが1の場合のみ1となり（これはエッジがカット内にあることを意味する）、それ以外は0となる。カットのエッジを最大化する問題は次のように定式化できる

\max_{x\in \{0, 1\}^n} \sum_{(i,j)} x_i + x_j - 2 x_i x_j,

の最小化として書き直すことができる

\min_{x\in \{0, 1\}^n} \sum_{(i,j)} 2 x_i x_j - x_i - x_j.

この場合の $f(x)$ の最小値は、カットが横切る辺の数が最大になるときである。ご覧の通り、量子コンピューティングに関連するものはまだ何もない。この問題を量子コンピューターが理解できるように定式化する必要がある。

$n=5$ のノードでグラフを作成し、問題を初期化する。

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)

Output:

3.2 ステップ1. 古典的な入力を量子問題にマッピングする

このパターンの最初のステップは、古典的な問題（グラフ）を量子回路と 演算子にマッピングすることである。そのためには、主に3つのステップを踏む必要がある：

一連の数学的再定式化を利用し、2次制約なし2値最適化（QUBO）問題表記を用いてこの問題を表現する。
最適化問題を、基底状態がコスト関数を最小化する解に対応するハミルトニアンとして書き直す。
量子アニーリングに似たプロセスで、このハミルトニアンの基底状態を準備する量子回路を作る。

注： QAOA手法では、最終的に、ハイブリッドアルゴリズムのコスト関数を表す演算子（ハミルトニアン ）と、問題の解の候補となる量子状態を表すパラメトライズ回路（アンサッツ ）を持ちたい。これらの候補状態からサンプリングし、コスト関数を用いて評価することができる。

グラフ → 最適化問題

マッピングの最初のステップは表記法の変更である：

\min_{x\in \{0, 1\}^n}x^T Q x,

ここで、 $Q$ は実数の $n\times n$ 行列、 $n$ はグラフのノード数、 $x$ は上で紹介したバイナリ変数のベクトル、 $x^T$ はベクトル $x$ の転置を示す。

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問題をハミルトニアン （ここでは系のエネルギーを表す行列）として再定式化することができる：

H_C=\sum_{ij}Q_{ij}Z_iZ_j + \sum_i b_iZ_i.

QAOA問題からハミルトニアンへの再定式化ステップ

この方法でQAOA問題がどのように書き換えられるかを示すために、まず、バイナリ変数 $x_i$ を新しい変数セット $z_i\in\{-1, 1\}$ に置き換える

x_i = \frac{1-z_i}{2}.

ここで、 $x_i$ が $0$ であるならば、 $z_i$ は $1$ でなければならないことがわかる。 $x_i$ を最適化問題( $x^TQx$ )の $z_i$ に代入すると、等価な定式化が得られる。

x^TQx=\sum_{ij}Q_{ij}x_ix_j \\ =\frac{1}{4}\sum_{ij}Q_{ij}(1-z_i)(1-z_j) \\=\frac{1}{4}\sum_{ij}Q_{ij}z_iz_j-\frac{1}{4}\sum_{ij}(Q_{ij}+Q_{ji})z_i + \frac{n^2}{4}.

ここで、 $b_i=-\sum_{j}(Q_{ij}+Q_{ji})$ を定義し、プレファクターと定数（ $n^2$ ）の項を取り除くと、同じ最適化問題の2つの等価な定式化に到達する。

min_{x\in\{0,1\}^n} x^TQx\Longleftrightarrow \min_{z\in\{-1,1\}^n}z^TQz + b^Tz

ここで、 $b$ は $Q$ に依存する。 $z^TQz + b^Tz$ を得るために、1/4 の係数と $n^2$ の定数オフセットを取り除いた。

さて、問題の量子定式化を得るためには、 $z_i$ 変数を Pauli $Z$ 行列、例えば、 $2\times 2$ 形式の行列に昇格させます

Z_i = \begin{pmatrix}1 & 0 \\ 0 & -1\end{pmatrix}.

これらの行列を上記の最適化問題に代入すると、次のようなハミルトニアンが得られる

H_C=\sum_{ij}Q_{ij}Z_iZ_j + \sum_i b_iZ_i.

また、 $Z$ の行列は量子コンピュータの計算空間、つまりサイズ $2^n\times 2^n$ のヒルベルト空間に埋め込まれていることを思い出してください。したがって、 $Z_iZ_j$ のような用語は、 $2^n\times 2^n$ ヒルベルト空間に埋め込まれたテンソル積 $Z_i\otimes Z_j$ として理解する必要があります。例えば、5つの決定変数を持つ問題では、 $Z_1Z_3$ という用語は、 $I\otimes Z_3\otimes I\otimes Z_1\otimes I$ （ $I$ は $2\times 2$ の恒等行列）を意味すると理解される。

このハミルトニアンは コスト関数ハミルトニアン と呼ばれ、その基底状態が はコスト関数 $f(x)$ の解に対応するという性質を持っています。したがって、最適化問題を解くためには、 $H_C$ の基底状態（またはそれと重なりが大きい状態）を量子コンピュータに用意する必要があります。そして、この状態からサンプリングすると、高い確率で $\min~f(x)$ の解が得られる。

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)

Output:

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

ハミルトニアン → 量子回路

ハミルトニアン $H_C$ には、あなたの問題の量子論的定義が含まれている。これで、量子コンピューターから良い解をサンプリングするのに役立つ量子回路を作ることができる。 QAOAは量子アニーリングに着想を得ており、量子回路に演算子を交互に重ねて適用する。

一般的な考え方は、既知のシステムの基底状態（ $H^{\otimes n}|0\rangle$ ）から始めて、興味のあるコスト演算子の基底状態にシステムを誘導することである。これは、角度 $\gamma_1,...,\gamma_p$ と $\beta_1,...,\beta_p~$ を持つ演算子 $\exp\{-i\gamma_k H_C\}$ と $\exp\{-i\beta_k H_m\}$ を適用することによって行われる。

生成される量子回路は、 $\gamma_i$ と $\beta_i$ でパラメトリック化されているので、 $\gamma_i$ と $\beta_i$ の異なる値を試して、結果の状態からサンプリングすることができます。

今回は、 $\gamma_1$ と $\beta_1$ の2つのパラメーターを含む、1つのQAOAレイヤーの例を試してみる。

from qiskit.circuit.library import QAOAAnsatz

circuit = QAOAAnsatz(cost_operator=cost_hamiltonian, reps=1)
circuit.measure_all()
circuit.draw("mpl")

Output:

circuit.decompose(reps=3).draw("mpl", fold=-1)

Output:

circuit.parameters

Output:

ParameterView([ParameterVectorElement(β[0]), ParameterVectorElement(γ[0])])

3.3 ステップ2. 量子ハードウェア実行のための回路最適化

上記の回路は、量子アルゴリズムを考えるのに便利な抽象化された一連の機能を含んでいるが、ハードウェア上で実行することは不可能である。 QPU上で動作させるためには、回路はパターンのトランスパイルステップまたは回路最適化ステップを構成する一連の操作を受ける必要がある。

Qiskitライブラリは、幅広い回路変換に対応する一連のトランスピレーション・パスを提供します。回路が目的に応じて最適化されていることを確認する必要がある。

移籍には、次のようないくつかのステップがある：

回路内の量子ビット（決定変数など）をデバイス上の物理量子ビットに初期マッピングする。
量子回路の命令を、バックエンドが理解できるハードウェアネイティブな命令にアンロールする。
相互作用する回路内のあらゆる量子ビットを、互いに隣接する物理量子ビットにルーティングする。
動的デカップリングによるノイズ抑制のための単一量子ビットゲートの追加によるエラー抑制。

トランスピレーションの詳細については、ドキュメントをご覧ください。

以下のコードは、 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)

Output:

  service = QiskitRuntimeService(channel="ibm_quantum_platform")

<IBMBackend('ibm_strasbourg')>

3.4 ステップ3。 IBM Quantum プリミティブを使用して実行する

QAOAワークフローでは、最適なQAOAパラメータは反復最適化ループで求められる。このループでは、一連の回路評価が実行され、古典的なオプティマイザを使用して最適な $\beta_k$ および $\gamma_k$ パラメータが求められる。この実行ループは以下のステップを経て実行される：

初期パラメータの定義
最適化ループと回路のサンプリングに使用したプリミティブを含む新しい Session をインスタンス化する
最適なパラメータのセットが見つかったら、回路を最後に1回実行し、後処理ステップで使用する最終分布を得る。

初期パラメータで回路を定義する

任意に選んだパラメータでスタートする。

initial_gamma = np.pi
initial_beta = np.pi / 2
init_params = [initial_gamma, initial_beta]

バックエンドと実行プリミティブを定義する

IBM Quantum プリミティブを使用して、 IBM® バックエンドと連携します。 2つのプリミティブは「サンプラー」と「エスティメーター」であり、どのプリミティブを選択するかは、量子コンピュータ上でどのような測定を実行したいかによって決まります。 $H_C$ を最小化するには、コスト関数の値が単に $\langle H_C \rangle$ の期待値であるため、Estimator を使用します。

実行

プリミティブは、量子デバイス上でワークロードをスケジューリングするための様々な実行モードを提供し、QAOAワークフローはセッション内で繰り返し実行される。

サンプラーベースのコスト関数を 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)

Output:

 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()

Output:

回路の最適なパラメータを見つけたら、これらのパラメータを割り当て、最適化されたパラメータで得られた最終分布をサンプリングすることができます。グラフの最適カットに対応するビット列測定の確率分布であるため、ここで Sampler プリミティブを使用する。

注：これは、コンピューター内に量子状態（ $\psi$ ）を用意し、それを測定することを意味する。測定は、状態を単一の計算基礎状態（例えば、 010101110000... ）に折り畳む。これは、最初の最適化問題に対する解の候補 $x$ （タスクによっては $\max f(x)$ または $\min f(x)$ ）に対応する。

optimized_circuit = candidate_circuit.assign_parameters(result.x)
optimized_circuit.draw("mpl", fold=False, idle_wires=False)

Output:

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)

Output:

{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)

Output:

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()

Output:

最適なカットを可視化する

最適なビット列から、このカットを元のグラフ上に可視化することができる。

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)

Output:

そしてカットの価値を計算する。ノイズのため最適解ではない（最適解のカット値は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)

Output:

The value of the cut is: 5

これで小規模のQAOAチュートリアルは終了。 QAOAをユーティリティ・スケールで適応させる方法は、"パート2：スケールアップ！"で学ぶことができる量子近似最適化アルゴリズムチュートリアルの

# Check Qiskit version
import qiskit

qiskit.__version__

Output:

'2.0.2'

このページは役に立ちましたか？

バグや誤字の報告、またはコンテンツの要求はGitHubで行ってください。