コスト関数

このレッスンでは、 コスト関数を評価する方法を学ぶ：

まず、 IBM Quantum のプリミティブについて学びます
コスト関数 $C(\vec\theta)$ を定義する。これは問題固有の関数で、オプティマイザが最小化（または最大化）する問題のゴールを定義する
IBM Quantum のプリミティブを用いて測定戦略を策定し、速度と精度のバランスを最適化する

プリミティブ

古典的であれ量子的であれ、すべての物理系は異なる状態で存在することができる。例えば、道路を走る車は、その状態を特徴づける一定の質量、位置、速度、加速度を持つことができる。同様に、量子系も異なる構成や状態を持つことができるが、測定や状態変化の扱い方が古典系とは異なる。これにより、 重ね合わせや エンタングルメントといった量子力学特有の性質が生まれる。車の状態を速度や加速度などの物理的特性で表現できるように、量子系の状態も数学的対象である観測量を用いて表現することができる。

量子力学では、状態は正規化された複素列ベクトル、すなわちケット ( $|\psi\rangle$ )で表され、観測量はケットに作用するエルミート線形作用素( $\hat{H}=\hat{H}^{\dagger}$ )である。観測値の固有ベクトル( $|\lambda\rangle$ )は固有状態と呼ばれる。ある観測量をその固有状態( $|\lambda\rangle$ )の一つについて測定すると、対応する固有値( $\lambda$ )が読み出される。

量子システムの測定方法と測定可能なものに疑問をお持ちの場合、Qiskitは2つの基本的かつ基本的な操作またはデータ型。 QiskitにはSamplerとEstimatorプリミティブがあり、複雑なアルゴリズムのワークロードを簡単に構築するためのビルディングブロックとして機能します。：

Sampler:量子状態 $|\psi\rangle$ が与えられたとき、このプリミティブは計算可能な基底状態の確率を求める。
Estimator:量子観測量 $\hat{H}$ と状態 $|\psi\rangle$ が与えられたとき、このプリミティブは $\hat{H}$ の期待値を計算する。

サンプラープリミティブ

Sampler プリミティブは、状態 $|\psi\rangle$ を準備する量子回路が与えられた場合、計算基礎から各可能な状態 $|k\rangle$ を得る確率を計算します。次のように計算します

p_k = |\langle k | \psi \rangle|^2 \quad \forall k \in \mathbb{Z}_2^n \equiv \{0,1,\cdots,2^n-1\},

ここで、 $n$ は量子ビット数であり、 $k$ は、あらゆる可能な出力バイナリ文字列の整数表現 $\{0,1\}^n$ （つまり、整数の基数 $2$ ）である。

IBM Quantum は、量子デバイス上でこの回路を複数回実行 Sampler し、各実行ごとに測定を行い、復元されたビット列から確率分布を再構築します。実行するラン（またはショット ）の回数が増えれば増えるほど、結果はより正確になりますが、それにはより多くの時間と量子リソースが必要となります。

しかし、可能な出力の数は、量子ビットの数 $n$ （つまり、 $2^n$ ）とともに指数関数的に増加するため、密な確率分布を捉えるためには、ショットの数も同様に指数関数的に増加する必要がある。したがって、 Sampler は、 スパースな確率分布に対してのみ効率的である。ここで、ターゲット状態 $|\psi\rangle$ は、計算基底状態の線形結合として表現可能でなければならず、項の数は量子ビット数に応じて最大でも多項式に増加する：

|\psi\rangle = \sum^{\text{Poly}(n)}_k w_k |k\rangle.

また、 Sampler 、可能な全状態のサブセットを表す回路のサブセクションから確率を検索するように構成することもできる。

推定器プリミティブ

Estimator プリミティブは、量子状態 $|\psi\rangle$ に対する観測値 $\hat{H}$ の期待値を計算する。ここで、観測値の確率は $p_\lambda = |\langle\lambda|\psi\rangle|^2$ のように表すことができ、 $|\lambda\rangle$ は観測値 $\hat{H}$ の固有状態である。期待値は、状態 $|\psi\rangle$ の測定のすべての可能な結果 $\lambda$ (つまり、観測値の固有値) の平均として定義され、対応する確率で重み付けされます：

\langle\hat{H}\rangle_\psi := \sum_\lambda p_\lambda \lambda = \langle \psi | \hat{H} | \psi \rangle

しかし、観測量の本質基底がわからない場合が多いため、その期待値を計算できるとは限らない。 IBM Quantum は、観測量を、その固有基底が既知である他の観測量の組み合わせに分解することで、実在する量子デバイスにおける期待値を推定するために、複雑な代数的手法 Estimator を用いる。

もっと簡単に言えば、 Estimator 、測定方法を知らない観測可能なものは、量子コンピュータで量子状態を表現・操作するために一般的に使用される行列の集合で、恒等行列と3つのパウリ行列（X、Y、Z）で構成される。と呼ばれる、より単純で測定可能な観測可能なものに分解される。

どの演算子も、 $4^n$ パウリ演算子の組み合わせとして表現できる。

\hat{P}_k := \sigma_{k_{n-1}}\otimes \cdots \otimes \sigma_{k_0} \quad \forall k \in \mathbb{Z}_4^n \equiv \{0,1,\cdots,4^n-1\}, \\

条件

\hat{H} = \sum^{4^n-1}_{k=0} w_k \hat{P}_k

ここで、 $n$ は量子ビット数、 $k \equiv k_{n-1} \cdots k_0$ は $k_l \in \mathbb{Z}_4 \equiv \{0, 1, 2, 3\}$ （つまり、整数の底 $4$ ）、 $(\sigma_0, \sigma_1, \sigma_2, \sigma_3) := (I, X, Y, Z)$ 。

この分解を行った後、 Estimator 、（元の回路から）各観測値 $\hat{P}_k$ に対して新しい回路 $V_k|\psi\rangle$ を導出し、計算基底におけるパウリ観測値を効果的に対角化し、それを測定する。パウリ観測量を簡単に測定できるのは、 $V_k$ 。

$\hat{P}_{k}$ それぞれについて、 Estimator 量子デバイス上で対応する回路を複数回実行し、計算基礎で出力状態を測定し、各可能な出力 $j$ を得る確率 $p_{kj}$ を計算する。次に、各出力に対応する $P_k$ の固有値 $\lambda_{kj}$ $j$ を探し、 $w_k$ を乗算し、すべての結果を合計して、与えられた状態 $|\psi\rangle$ に対する観測値 $\hat{H}$ の期待値を求める。

\langle\hat{H}\rangle_\psi = \sum_{k=0}^{4^n-1} w_k \sum_{j=0}^{2^n-1}p_{kj} \lambda_{kj},

$4^n$ のパウリの期待値を計算するのは現実的ではない（つまり、計算量が指数関数的に増加する）ため、 $w_k$ の大部分がゼロである場合（つまり、密なパウリ分解ではなく疎なパウリ分解の場合）にのみ、効率的に処理 Estimator できる。形式的には、この計算を効率的に解くためには、ゼロでない項の数が量子ビット数 $n$ に対して多項式級以下に増加しなければならない、と述べられる： $\hat{H} = \sum^{\text{Poly}(n)}_k w_k \hat{P}_k$ 。

読者は、 Sampler で説明したように、確率 1つまたは複数の物事について、いくつかの測定を行うプロセス。も効率的である必要がある、という暗黙の前提に気づくかもしれない

\langle\hat{H}\rangle_\psi = \sum_{k}^{\text{Poly}(n)} w_k \sum_{j}^{\text{Poly}(n)}p_{kj} \lambda_{kj}.

期待値を計算するためのガイド付き例

1量子ビットの状態 $|+\rangle := H|0\rangle = \frac{1}{\sqrt{2}}(|0\rangle + |1\rangle)$ 、観測可能であるとする

\begin{aligned} \hat{H} & = \begin{pmatrix} -1 & 2 \\ 2 & 1 \\ \end{pmatrix}\\[1mm] & = 2X - Z \end{aligned}

その理論的な期待値は、 $\langle\hat{H}\rangle_+ = \langle+|\hat{H}|+\rangle = 2$ となる。

この観測値を測定する方法を知らないので、期待値を直接計算することはできず、 $\langle\hat{H}\rangle_+ = 2\langle X \rangle_+ - \langle Z \rangle_+$ として表現し直す必要がある。 $\langle+|X|+\rangle = 1$ , と $\langle+|Z|+\rangle = 0$ .

$\langle X \rangle_+$ と $\langle Z \rangle_+$ を直接計算する方法を見てみよう。 $X$ と $Z$ は交わらない（つまり、同じ固有ベーシスを共有しない）ので、同時に測定することはできず、したがって補助回路が必要となる：

from qiskit import QuantumCircuit
from qiskit.quantum_info import SparsePauliOp

# The following code will work for any other initial single-qubit state and observable
original_circuit = QuantumCircuit(1)
original_circuit.h(0)

H = SparsePauliOp(["X", "Z"], [2, -1])

aux_circuits = []
for pauli in H.paulis:
    aux_circ = original_circuit.copy()
    aux_circ.barrier()
    if str(pauli) == "X":
        aux_circ.h(0)
    elif str(pauli) == "Y":
        aux_circ.sdg(0)
        aux_circ.h(0)
    else:
        aux_circ.id(0)
    aux_circ.measure_all()
    aux_circuits.append(aux_circ)

original_circuit.draw("mpl")

Output:

# Auxiliary circuit for X
aux_circuits[0].draw("mpl")

Output:

# Auxiliary circuit for Z
aux_circuits[1].draw("mpl")

Output:

Sampler を使って手動で計算を行い、その結果を Estimator で確認することができる：

from qiskit.primitives import StatevectorSampler, StatevectorEstimator
from qiskit.result import QuasiDistribution
import numpy as np


## SAMPLER
shots = 10000
sampler = StatevectorSampler()
job = sampler.run(aux_circuits, shots=shots)

# Run the sampler job and step through results
expvals = []
for index, pauli in enumerate(H.paulis):
    data_pub = job.result()[index].data
    bitstrings = data_pub.meas.get_bitstrings()
    counts = data_pub.meas.get_counts()
    quasi_dist = QuasiDistribution(
        {outcome: freq / shots for outcome, freq in counts.items()}
    )

    # Use the probabilities and known eigenvalues of Pauli operators to estimate
    # the expectation value.
    val = 0

    if str(pauli) == "X":
        val += -1 * quasi_dist.get(1, 0)
        val += 1 * quasi_dist.get(0, 0)

    if str(pauli) == "Y":
        val += -1 * quasi_dist.get(1, 0)
        val += 1 * quasi_dist.get(0, 0)

    if str(pauli) == "Z":
        val += 1 * quasi_dist.get(0, 0)
        val += -1 * quasi_dist.get(1, 0)

    expvals.append(val)

# Print expectation values

print("Sampler results:")
for pauli, expval in zip(H.paulis, expvals):
    print(f"  >> Expected value of {str(pauli)}: {expval:.5f}")

total_expval = np.sum(H.coeffs * expvals).real
print(f"  >> Total expected value: {total_expval:.5f}")

# Use estimator for comparison
observables = [
    *H.paulis,
    H,
]  # Note: run for individual Paulis as well as full observable H

estimator = StatevectorEstimator()
job = estimator.run([(original_circuit, observables)])
estimator_expvals = job.result()[0].data.evs

# Print results
print("Estimator results:")
for obs, expval in zip(observables, estimator_expvals):
    if obs is not H:
        print(f"  >> Expected value of {str(obs)}: {expval:.5f}")
    else:
        print(f"  >> Total expected value: {expval:.5f}")

Output:

Sampler results:
  >> Expected value of X: 1.00000
  >> Expected value of Z: 0.00420
  >> Total expected value: 1.99580
Estimator results:
  >> Expected value of X: 1.00000
  >> Expected value of Z: 0.00000
  >> Total expected value: 2.00000

数学的厳密性（任意）

$\hat{H}$ , $|\psi\rangle = \sum_\lambda a_\lambda |\lambda\rangle$ の固有状態の基底に関して $|\psi\rangle$ を表現すると、以下のようになる：

\begin{aligned} \langle \psi | \hat{H} | \psi \rangle & = \bigg(\sum_{\lambda'}a^*_{\lambda'} \langle \lambda'|\bigg) \hat{H} \bigg(\sum_{\lambda} a_\lambda | \lambda\rangle\bigg)\\[1mm] & = \sum_{\lambda}\sum_{\lambda'} a^*_{\lambda'}a_{\lambda} \langle \lambda'|\hat{H}| \lambda\rangle\\[1mm] & = \sum_{\lambda}\sum_{\lambda'} a^*_{\lambda'}a_{\lambda} \lambda \langle \lambda'| \lambda\rangle\\[1mm] & = \sum_{\lambda}\sum_{\lambda'} a^*_{\lambda'}a_{\lambda} \lambda \cdot \delta_{\lambda, \lambda'}\\[1mm] & = \sum_\lambda |a_\lambda|^2 \lambda\\[1mm] & = \sum_\lambda p_\lambda \lambda\\[1mm] \end{aligned}

対象となる観測量 $\hat{H}$ の固有値や固有状態がわからないため、まずその対角化について検討する必要がある。 $\hat{H}$ がエルミートとは、それ自身の共役転置と等しい正方行列、または自己共役である線形作用素のことである。であることから、 $\hat{H}=V^\dagger \Lambda V$ を満たすユニタリ変換 $V$ が存在する。ここで、 $\Lambda$ は対角固有値行列であるため、 $j\neq k$ ならば $\langle j | \Lambda | k \rangle = 0$ となり、 $\langle j | \Lambda | j \rangle = \lambda_j$ となる。

これは期待値が次のように書き換えられることを意味する：

\begin{aligned} \langle\psi|\hat{H}|\psi\rangle & = \langle\psi|V^\dagger \Lambda V|\psi\rangle\\[1mm] & = \langle\psi|V^\dagger \bigg(\sum_{j=0}^{2^n-1} |j\rangle \langle j|\bigg) \Lambda \bigg(\sum_{k=0}^{2^n-1} |k\rangle \langle k|\bigg) V|\psi\rangle\\[1mm] & = \sum_{j=0}^{2^n-1} \sum_{k=0}^{2^n-1}\langle\psi|V^\dagger |j\rangle \langle j| \Lambda |k\rangle \langle k| V|\psi\rangle\\[1mm] & = \sum_{j=0}^{2^n-1}\langle\psi|V^\dagger |j\rangle \langle j| \Lambda |j\rangle \langle j| V|\psi\rangle\\[1mm] & = \sum_{j=0}^{2^n-1}|\langle j| V|\psi\rangle|^2 \lambda_j\\[1mm] \end{aligned}

\langle\psi|\hat{H}|\psi\rangle = \sum_{j=0}^{2^n-1} p_j \lambda_j.

確率は、 $|\psi\rangle$ ではなく、 $V |\psi\rangle$ 状態から取られることに注意することが非常に重要である。これが、行列 $V$ が絶対に必要な理由である。

行列 $V$ と固有値 $\Lambda$ をどうやって求めるのか、不思議に思うかもしれない。もしすでに固有値がわかっていれば、変分アルゴリズムの目的はこれらの固有値を求めることなので、量子コンピュータを使う必要はないだろう $\hat{H}$ 。

幸いなことに、これを回避する方法があります。任意の $2^n \times 2^n$ 行列は、 $n$ のパウリ行列と単位行列の $4^n$ のテンソル積の線形結合として表すことができ、これらはすべてエルミートかつユニタリであり、 $V$ および $\Lambda$ が既知です。これが、 IBM Quantum が内部で実行している処理であり Estimator 、任意のオブジェクト Operator を.に分解しています SparsePauliOp。

以下は使用可能な演算子である：

\begin{array}{c|c|c|c} \text{Operator} & \sigma & V & \Lambda \\[1mm] \hline I & \sigma_0 = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} & V_0 = I & \Lambda_0 = I = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} \\[4mm] X & \sigma_1 = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} & V_1 = H =\frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix} & \Lambda_1 = \sigma_3 = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \\[4mm] Y & \sigma_2 = \begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix} & V_2 = HS^\dagger =\frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix}\cdot \begin{pmatrix} 1 & 0 \\ 0 & -i \end{pmatrix} = \frac{1}{\sqrt{2}} \begin{pmatrix} 1 & -i \\ 1 & i \end{pmatrix}\quad & \Lambda_2 = \sigma_3 = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \\[4mm] Z & \sigma_3 = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} & V_3 = I & \Lambda_3 = \sigma_3 = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \end{array}

そこで、 $\hat{H}$ 、ポーリスとアイデンティティに関して書き直してみよう：

\hat{H} = \sum_{k_{n-1}=0}^3... \sum_{k_0=0}^3 w_{k_{n-1}...k_0} \sigma_{k_{n-1}}\otimes ... \otimes \sigma_{k_0} = \sum_{k=0}^{4^n-1} w_k \hat{P}_k,

ここで $k = \sum_{l=0}^{n-1} 4^l k_l \equiv k_{n-1}...k_0$ は $k_{n-1},...,k_0\in \{0,1,2,3\}$ （つまりベース $4$ ）、 $\hat{P}_{k} := \sigma_{k_{n-1}}\otimes ... \otimes \sigma_{k_0}$ ：

\begin{aligned} \langle\psi|\hat{H}|\psi\rangle & = \sum_{k=0}^{4^n-1} w_k \sum_{j=0}^{2^n-1}|\langle j| V_k|\psi\rangle|^2 \langle j| \Lambda_k |j\rangle \\[1mm] & = \sum_{k=0}^{4^n-1} w_k \sum_{j=0}^{2^n-1}p_{kj} \lambda_{kj}, \\[1mm] \end{aligned}

ここで、 $V_k := V_{k_{n-1}}\otimes ... \otimes V_{k_0}$ および $\Lambda_k := \Lambda_{k_{n-1}}\otimes ... \otimes \Lambda_{k_0}$ であり、次の関係が成り立つ： $\hat{P_k}=V_k^\dagger \Lambda_k V_k$ 。

コスト関数

一般に、コスト関数は問題のゴールと、そのゴールに対して試行状態がどの程度うまくいっているかを記述するために使われる。この定義は、化学、機械学習、金融、最適化など、さまざまな例に適用できる。

系の基底状態を見つける簡単な例を考えてみよう。我々の目的は、エネルギーを表す観測値（ハミルトニアン $\hat{\mathcal{H}}$ ）の期待値を最小化することである：

\min_{\vec\theta} \langle\psi(\vec\theta)|\hat{\mathcal{H}}|\psi(\vec\theta)\rangle

Estimator を使って期待値を評価し、この値をオプティマイザに渡して最小化することができる。最適化が成功すれば、最適なパラメータ値のセット $\vec\theta^*$ が返され、そこから提案された解の状態 $|\psi(\vec\theta^*)\rangle$ を構築し、観測された期待値を $C(\vec\theta^*)$ として計算することができる。

私たちが考えている限られた状態のセットに対してのみ、コスト関数を最小化することができることに注目してほしい。このことから、2つの可能性が考えられる：

私たちのアナザッツは、探索空間全体の解の状態を定義していません ：もしそうであれば、オプティマイザーは解を見つけることができないので、探索空間をより正確に表現できるかもしれない他のアサッツを試してみる必要があります。
我々のオプティマイザーは、この有効解を見つけることができない ：最適化はグローバルに定義することもローカルに定義することもできる。これが何を意味するかは、後のセクションで探ることにしよう。

結局のところ、古典的な最適化ループを実行することになるが、コスト関数の評価は量子コンピューターに頼ることになる。この観点から、最適化を純粋に古典的な取り組みと考えることもできる。そこでは、オプティマイザがコスト関数を評価する必要があるたびに、特定の機能を果たすが、その内部構造が不明な仮想的な装置またはソフトウェア・コンポーネント。ユーザーはブラックボックスの入力と出力しか知らず、ブラックボックスがどのように入力を処理して出力を生成するかは知らない。。

def cost_func_vqe(params, circuit, hamiltonian, estimator):
    """Return estimate of energy from estimator

    Parameters:
        params (ndarray): Array of ansatz parameters
        ansatz (QuantumCircuit): Parameterized ansatz circuit
        hamiltonian (SparsePauliOp): Operator representation of Hamiltonian
        estimator (Estimator): Estimator primitive instance

    Returns:
        float: Energy estimate
    """
    pub = (circuit, hamiltonian, params)
    cost = estimator.run([pub]).result()[0].data.evs
    return cost

from qiskit.circuit.library import TwoLocal

observable = SparsePauliOp.from_list([("XX", 1), ("YY", -3)])

reference_circuit = QuantumCircuit(2)
reference_circuit.x(0)

variational_form = TwoLocal(
    2,
    rotation_blocks=["rz", "ry"],
    entanglement_blocks="cx",
    entanglement="linear",
    reps=1,
)
ansatz = reference_circuit.compose(variational_form)

theta_list = (2 * np.pi * np.random.rand(1, 8)).tolist()
ansatz.decompose().draw("mpl")

Output:

まず、シミュレーター（ StatevectorEstimator ）を使ってこれを実行する。通常、これはデバッグのために推奨されることだが、我々はデバッグの実行後、すぐに実際の量子ハードウェア上での計算を行う。最先端のスーパーコンピューティング設備がなければ、もはや古典的なシミュレーションが不可能な問題が増えている。

estimator = StatevectorEstimator()
cost = cost_func_vqe(theta_list, ansatz, observable, estimator)
print(cost)

Output:

[-0.58744589]

それでは、実際の量子コンピュータ上で実行を進めていきます。構文の変更点にご注意ください。 pass_manager に関連する手順については、次の例でさらに詳しく説明します。変分アルゴリズムにおいて特に重要な手順の一つは、 IBM Quantum Compute Service セッションの利用です。セッションを開始すると、パラメータが更新されるたびに新しいキューで待機することなく、変分アルゴリズムを複数回実行できるようになります。待ち時間が長い場合や、反復処理が何度も必要な場合には、これが重要になります。なお、Open PlanのユーザーはQuantum Computeセッションを利用できません。セッションにアクセスできない場合は、一度に送信する反復回数を減らし、次回以降の実行で再利用できるよう、最新のパラメータを保存しておくことができます。反復処理の送信回数が多すぎたり、キューの待ち時間が長すぎたりすると、エラーコード 1217 が発生する場合があります。これは、ジョブの送信間隔が長すぎることを示しています。

# Estimated usage: < 1 min. Benchmarked at 7 seconds on an Eagle processor
# Load necessary packages:

from qiskit_ibm_runtime import (
    QiskitRuntimeService,
    Session,
    EstimatorOptions,
    EstimatorV2 as Estimator,
)
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

# Select the least busy backend:

service = QiskitRuntimeService()
backend = service.least_busy(
    operational=True, min_num_qubits=ansatz.num_qubits, simulator=False
)
# Or get a specific backend:
# backend = service.backend("ibm_brisbane")

# Use a pass manager to transpile the circuit and observable for the specific backend being used:

pm = generate_preset_pass_manager(backend=backend, optimization_level=1)
isa_ansatz = pm.run(ansatz)
isa_observable = observable.apply_layout(layout=isa_ansatz.layout)


# Set estimator options
estimator_options = EstimatorOptions(resilience_level=1, default_shots=10_000)

# Open a Quantum Compute session:

with Session(backend=backend) as session:
    estimator = Estimator(mode=session, options=estimator_options)
    cost = cost_func_vqe(theta_list, isa_ansatz, isa_observable, estimator)

session.close()
print(cost)

上記の2つの計算で得られた値は非常に似ていることに注意してほしい。結果を向上させるテクニックについては、さらに後述する。

非物理システムへの例示マッピング

最大切断（max-cut）問題は、グラフの頂点を2つの互いに素な集合に分割し、その2つの集合間の辺の数を最大化する組み合わせ最適化問題である。より形式的に言えば、無向グラフ $G=(V,E)$ （ $V$ を頂点の集合、 $E$ を辺の集合とする）において、最大切断問題とは、頂点を2つの互いに素な部分集合 $S$ と $T$ に分割し、一方の端点が $S$ にあり、もう一方の端点が $T$ にある辺の数を最大化するよう求める問題である。

最大カット法は、クラスタリング、ネットワーク設計、相転移など、さまざまな問題の解決に応用することができます。まずは問題グラフを作成することから始めましょう：

import rustworkx as rx
from rustworkx.visualization import mpl_draw

n = 4
G = rx.PyGraph()
G.add_nodes_from(range(n))
# The edge syntax is (start, end, weight)
edges = [(0, 1, 1.0), (0, 2, 1.0), (0, 3, 1.0), (1, 2, 1.0), (2, 3, 1.0)]
G.add_edges_from(edges)

mpl_draw(
    G, pos=rx.shell_layout(G), with_labels=True, edge_labels=str, node_color="#1192E8"
)

Output:

この問題は二項最適化問題として表現できる。各ノード $0 \leq i < n$ 、ここで $n$ はグラフのノード数（この場合 $n=4$ ）で、バイナリ変数 $x_i$ を考える。この変数は、ノード $i$ が $1$ とラベル付けするグループの1つであれば値 $1$ を持ち、 $0$ とラベル付けするもう1つのグループであれば値 $0$ を持つ。また、ノード $i$ からノード $j$ に向かう辺の重みを $w_{ij}$ （隣接行列 $w$ の要素 $(i,j)$ ）とする。グラフは無向グラフなので、 $w_{ij}=w_{ji}$ 。そうすると、我々の問題は次のコスト関数を最大化するものとして定式化できる：

\begin{aligned} C(\vec{x}) & =\sum_{i,j=0}^n w_{ij} x_i(1-x_j)\\[1mm] & = \sum_{i,j=0}^n w_{ij} x_i - \sum_{i,j=0}^n w_{ij} x_ix_j\\[1mm] & = \sum_{i,j=0}^n w_{ij} x_i - \sum_{i=0}^n \sum_{j=0}^i 2w_{ij} x_ix_j \end{aligned}

量子コンピュータでこの問題を解くために、コスト関数を観測可能なものの期待値として表現する。しかし、Qiskitがネイティブに認める観測量はパウリ作用素で構成され、その固有値は $0$ と $1$ ではなく $1$ と $-1$ です。そのため、次のように変数を変更します：

ここで $\vec{x}=(x_0,x_1,\cdots ,x_{n-1})$ 。隣接行列 $w$ を使えば、すべてのエッジの重みに簡単にアクセスできる。これを用いてコスト関数を求める：

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

ということを意味している：

\begin{array}{lcl} x_i=0 & \rightarrow & z_i=1 \\ x_i=1 & \rightarrow & z_i=-1.\end{array}

つまり、最大化したい新しいコスト関数は次のようになる：

\begin{aligned} C(\vec{z}) & = \sum_{i,j=0}^n w_{ij} \bigg(\frac{1-z_i}{2}\bigg)\bigg(1-\frac{1-z_j}{2}\bigg)\\[1mm] & = \sum_{i,j=0}^n \frac{w_{ij}}{4} - \sum_{i,j=0}^n \frac{w_{ij}}{4} z_iz_j\\[1mm] & = \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2} - \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2} z_iz_j \end{aligned}

さらに、量子コンピュータの自然な傾向として、極大値ではなく極小値（通常は最低エネルギー）を見つけるので、 $C(\vec{z})$ を最大化するのではなく、最小化することになる：

-C(\vec{z}) = \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2} z_iz_j - \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2}

さて、最小化すべきコスト関数があり、その変数は値 $-1$ と $1$ を持つことができるので、次のようにパウリ $Z$ と類推することができる：

z_i \equiv Z_i = \overbrace{I}^{n-1}\otimes ... \otimes \overbrace{Z}^{i} \otimes ... \otimes \overbrace{I}^{0}

言い換えれば、変数 $z_i$ は、量子ビット $i$ に作用する $Z$ ゲートと等価になる。さらに

Z_i|x_{n-1}\cdots x_0\rangle = z_i|x_{n-1}\cdots x_0\rangle \rightarrow \langle x_{n-1}\cdots x_0 |Z_i|x_{n-1}\cdots x_0\rangle = z_i

となると、これから考えるのは観測値である：

\hat{H} = \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2} Z_iZ_j

これに独立項を後から加えることになる：

\texttt{offset} = - \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2}

この演算子は、エッジで結ばれたノード上のZ演算子を持つ項の線形結合である（0番目の量子ビットが最も右にあることを思い出してほしい）： $IIZZ + IZIZ + IZZI + ZIIZ + ZZII$ .いったん演算子が構築されれば、QAOAアルゴリズムのansatzは、Qiskit回路ライブラリの QAOAAnsatz 回路を使うことで簡単に構築できます。

from qiskit.circuit.library import QAOAAnsatz
from qiskit.quantum_info import SparsePauliOp

hamiltonian = SparsePauliOp.from_list(
    [("IIZZ", 1), ("IZIZ", 1), ("IZZI", 1), ("ZIIZ", 1), ("ZZII", 1)]
)


ansatz = QAOAAnsatz(hamiltonian, reps=2)
# Draw
ansatz.decompose(reps=3).draw("mpl")

Output:

# Sum the weights, and divide by 2

offset = -sum(edge[2] for edge in edges) / 2
print(f"""Offset: {offset}""")

Output:

Offset: -2.5

IBM Quantum 推定器は、ハミルトニアンとパラメータ化されたアンザッツを直接受け取り、必要なエネルギーを返すため、QAOAインスタンスのコスト関数は非常に単純です：

def cost_func(params, ansatz, hamiltonian, estimator):
    """Return estimate of energy from estimator

    Parameters:
        params (ndarray): Array of ansatz parameters
        ansatz (QuantumCircuit): Parameterized ansatz circuit
        hamiltonian (SparsePauliOp): Operator representation of Hamiltonian
        estimator (Estimator): Estimator primitive instance

    Returns:
        float: Energy estimate
    """
    pub = (ansatz, hamiltonian, params)
    cost = estimator.run([pub]).result()[0].data.evs
    #    cost = estimator.run(ansatz, hamiltonian, parameter_values=params).result().values[0]
    return cost

import numpy as np

x0 = 2 * np.pi * np.random.rand(ansatz.num_parameters)

estimator = StatevectorEstimator()
cost = cost_func_vqe(x0, ansatz, hamiltonian, estimator)
print(cost)

Output:

1.473098768180865

# Estimated usage: < 1 min, benchmarked at 6 seconds on ibm_osaka, 5-23-24
# Load some necessary packages:

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime import Session, EstimatorV2 as Estimator

# Select the least busy backend:

backend = service.least_busy(
    operational=True, min_num_qubits=ansatz.num_qubits, simulator=False
)

# Or get a specific backend:
# backend = service.backend("ibm_brisbane")

# Use a pass manager to transpile the circuit and observable for the specific backend being used:

pm = generate_preset_pass_manager(backend=backend, optimization_level=1)
isa_ansatz = pm.run(ansatz)
isa_hamiltonian = hamiltonian.apply_layout(layout=isa_ansatz.layout)

# Set estimator options
estimator_options = EstimatorOptions(resilience_level=1, default_shots=10_000)

# Open a Quantum Compute session:

with Session(backend=backend) as session:
    estimator = Estimator(mode=session, options=estimator_options)
    cost = cost_func_vqe(x0, isa_ansatz, isa_hamiltonian, estimator)

# Close session after done
session.close()
print(cost)

Output:

1.1120776913677988

アプリケーションでこの例を再び取り上げ、オプティマイザを活用して探索空間を反復する方法を探ります。一般的に言えば、これには以下が含まれる：

オプティマイザを活用して最適なパラメータを見つける
最適なパラメータをアサッツに結合して固有値を求める
固有値を問題定義に置き換える

測定戦略：速度対精度

前述したように、私たちはノイズの多い量子コンピュータをブラックボックスのオラクルとして使っています。ノイズは取得した値を非決定論的なものにし、ランダムな揺らぎをもたらします。これは、量子の有用性を徐々に探求し、量子の優位性へと前進していく中で、私たちが取り組まなければならない一般的な問題である：

シミュレーション・コストが回路の複雑さによってどのように変化するかを示すグラフ。古典的なコンピュータを使えば、指数関数的に成長する。量子エラーの緩和を使えば、それが有利になるクロスオーバーがあるはずだ。量子エラー訂正は、シミュレーション・コストをリニアに増加させることが可能であり、確実にアドバンテージにつながる。

IBM Quantum プリミティブのエラー抑制およびエラー緩和オプションを活用することで、ノイズに対処し、現在の量子コンピュータの有用性を最大限に引き出すことができます。

エラー抑制

エラー抑制とは、エラーを最小限に抑えるために、コンパイル時に回路を最適化および変換する手法を指します。これは基本的なエラー処理の手法であり、通常、全体の実行時間に基本的な実装に対して、新しい技術によってもたらされる余分なコスト。程度の処理オーバーヘッドが生じます。このオーバーヘッドには、量子ハードウェア上で実行するために回路をトランスパイルする以下の処理が含まれます：

量子システムで利用可能なネイティブゲートを使って回路を表現する
仮想量子ビットを物理量子ビットにマッピングする
接続要件に基づくSWAPの追加
1Q と 2Q ゲートの最適化
アイドル量子ビットに動的デカップリングを追加してデコヒーレンスの影響を防ぐ。

プリミティブでは、 optimization_level オプションを設定し、高度なトランスパイルオプションを選択することで、エラー抑制技術を使用することができます。後のコースでは、結果を向上させるためのさまざまな回路構成方法について掘り下げていきますが、ほとんどの場合、 optimization_level=3 を設定することをお勧めします。

ここでは、単純な理想的動作をする回路例を見て、トランスパイルプロセスにおける最適化の増加の価値を可視化する。

from qiskit.circuit import Parameter, QuantumCircuit
from qiskit.quantum_info import SparsePauliOp

theta = Parameter("theta")

qc = QuantumCircuit(2)
qc.x(1)
qc.h(0)
qc.cp(theta, 0, 1)
qc.h(0)
observables = SparsePauliOp.from_list([("ZZ", 1)])

qc.draw("mpl")

Output:

上記の回路は、 $[0,2\pi]$ のような適切な区間に位相を挿入すれば、与えられた観測値の正弦波の期待値を得ることができる。

## Setup phases
import numpy as np

phases = np.linspace(0, 2 * np.pi, 50)

# phases need to be expressed as a list of lists in order to work
individual_phases = [[phase] for phase in phases]

シミュレーターを使って、最適化されたトランスパレーションの有用性を示すことができる。エラー緩和の有用性を実証するために、実際のハードウェアを使用することについては後述する。 QiskitRuntimeService、実際のバックエンド（この場合はibm_brisbane）を取得し、 AerSimulator、そのバックエンド（ノイズの挙動を含む）をシミュレートする。

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_aer import AerSimulator

# get a real backend from the runtime service
service = QiskitRuntimeService()
backend = service.backend("ibm_brisbane")

# generate a simulator that mimics the real quantum system with the latest calibration results
backend_sim = AerSimulator.from_backend(backend)

これで、パスマネージャーを使用して、回路をバックエンドの「命令セットアーキテクチャ（ISA）」にトランスパイルすることができます。これは、 IBM Quantum Compute Serviceにおける新しい要件です。バックエンドに送信されるすべての回路は、そのバックエンドのターゲットの制約に準拠している必要があります。つまり、バックエンドのISA（そのデバイスが理解・実行できる命令セット）に基づいて記述されている必要があります。これらの目標制約は、デバイスのネイティブ基底ゲート、量子ビット間の接続性、および（該当する場合）パルスやその他の命令のタイミング仕様といった要因によって定義されます。

今回のケースでは、最適化レベルを0に設定した場合と、3に設定した場合の2回行う。毎回、Estimatorプリミティブを使用して、異なる位相の値における観測値の期待値を推定する。

# Import estimator and specify that we are using the simulated backend:

from qiskit_ibm_runtime import EstimatorV2 as Estimator

estimator = Estimator(mode=backend_sim)

circuit = qc

# Use a pass manager to transpile the circuit and observable for the backend being simulated.
# Start with no optimization:

from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

pm = generate_preset_pass_manager(backend=backend_sim, optimization_level=0)
isa_circuit = pm.run(circuit)
isa_observables = observables.apply_layout(layout=isa_circuit.layout)

noisy_exp_values = []
pub = (isa_circuit, isa_observables, [individual_phases])
cost = estimator.run([pub]).result()[0].data.evs
noisy_exp_values = cost[0]

# Repeat above steps, but now with optimization = 3:

exp_values_with_opt_es = []
pm = generate_preset_pass_manager(backend=backend_sim, optimization_level=3)
isa_circuit = pm.run(circuit)
isa_observables = observables.apply_layout(layout=isa_circuit.layout)

pub = (isa_circuit, isa_observables, [individual_phases])
cost = estimator.run([pub]).result()[0].data.evs
exp_values_with_opt_es = cost[0]

最後に結果をプロットしてみると、最適化なしでも計算の精度はかなり良かったが、最適化をレベル3に上げることで確実に向上していることがわかる。より深く、より複雑な回路では、最適化レベル0と3の差はより大きくなる可能性が高いことに注意されたい。これは、おもちゃの模型として使われる非常に単純な回路である。

import matplotlib.pyplot as plt

plt.plot(phases, noisy_exp_values, "o", label="opt=0")
plt.plot(phases, exp_values_with_opt_es, "o", label="opt=3")
plt.plot(phases, 2 * np.sin(phases / 2) ** 2 - 1, label="ideal")
plt.ylabel("Expectation")
plt.legend()
plt.show()

Output:

エラー軽減

エラー緩和とは、実行時にデバイスのノイズをモデル化することで、ユーザーが回路エラーを低減できるようにする手法を指します。通常、これに伴い、モデルの学習に関連する量子事前処理のオーバーヘッドと、生成されたモデルを用いて生データの誤差を軽減するための古典的後処理のオーバーヘッドが発生する。

IBM Quantum プリミティブのオプション resilience_level は、エラーに対する耐性をどの程度設定するかを指定します。レベルが高いほど、量子サンプリングのオーバーヘッドにより処理時間は長くなりますが、より正確な結果が得られます。レジリエンスレベルを使用することで、プリミティブクエリにエラー軽減を適用する際の、コストと精度のバランスを調整することができます。

誤差を軽減する技術を導入する場合、結果の通常は誤差によって引き起こされる、測定量の系統的なドリフト。、以前の、軽減されていないバイアスに比べて減少することが期待される。場合によっては、バイアスが消えることさえある。しかし、これには代償が伴う。推定量のバイアスを減らすと、統計的なばらつきが大きくなる（つまり分散が大きくなる）。これは、バイアスを減らすために必要な以上のオーバーヘッドをもたらすので、デフォルトでは行われない。以下の例のように、 options.executions.shots で回路あたりのショット数を調整することで、この動作に簡単にオプトインすることができる。

/var バイアスのトレードオフのように、分布が広がったり狭まったりすることを示す図。

このコースでは、 IBM Quantum のプリミティブが実行できるエラー軽減機能を、詳細な実装内容に踏み込むことなく概観するために、これらのエラー軽減モデルについて大まかに解説します。

ツイールド読み出しエラー消去（T-REx）

Twirled Readout Error Extinction（T-REx）は、量子測定の過程で発生するノイズを低減するために、パウリ・ツワーリングとして知られる技術を使用する。このテクニックは、特定のノイズの形を想定していないため、非常に一般的で効果的である。

全体のワークフロー

ランダムなビット反転でゼロ状態のデータを取得（測定前のパウリX）
ランダム化されたビット反転（測定前のパウリX）により、希望する（ノイズの多い）状態のデータを取得する
各データセットについて特殊関数を計算し、除算する。

これは options.resilience_level = 1 で設定できる。

ゼロノイズ外挿

ゼロノイズ外挿（ZNE）は、まず目的の量子状態を準備する回路のノイズを増幅し、いくつかの異なるノイズレベルの測定値を得て、それらの測定値を用いてノイズのない結果を推測することで機能する。

全体のワークフロー

複数のノイズファクターに対して回路ノイズを増幅
すべてのノイズ増幅回路を動かす
ゼロノイズ限界まで外挿する

ZNEのステップを示す図。ノイズはさまざまな要因によって人為的に増幅される。その後、ノイズがゼロのときの値を外挿する。

これは次のように設定 options.resilience_level = 2できます。さまざまな noise_factors、 noise_amplifiers、およびを検討することで extrapolators、これをさらに最適化することは可能ですが、これは本コースの範囲外となります。ここで説明しているこれらのオプションを、ぜひ試してみてください。

各手法には、それぞれ関連するオーバーヘッドがある。必要な量子計算の回数（時間）と結果の精度とのトレードオフである：

\begin{array}{c|c|c|c} \text{Methods} & R=1 \text{, T-REx} & R=2 \text{, ZNE} \\[1mm] \hline \text{Assumptions} & \text{None} & \text{Ability to scale noise} \\[1mm] \text{Qubit overhead} & 1 & 1 \\[1mm] \text{Sampling overhead} & 2 & N_{\text{noise-factors}} \\[1mm] \text{Bias} & 0 & \mathcal{O}(\lambda^{N_{\text{noise-factors}}}) \\[1mm] \end{array}

Quantum Compute のオプションを活用したリスク軽減と抑制

IBM Quantum Compute Service でエラーの軽減および抑制機能を使用しながら、期待値を計算する方法をご紹介します。以前とまったく同じ回路と観測変数を活用できますが、今回は最適化レベルをレベル2に固定したまま、使用しているレジリエンスやエラー軽減手法を調整します。このエラー軽減処理は、最適化ループ全体を通じて複数回行われます。

エラー緩和はシミュレーターでは利用できないので、この部分は実際のハードウェアで実行する。

# Estimated usage: 8 minutes, benchmarked on an Eagle processor, 5-23-24

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime import (
    Session,
    EstimatorOptions,
    EstimatorV2 as Estimator,
)

# We select the least busy backend

# Select the least busy backend
# backend = service.least_busy(
#    operational=True, min_num_qubits=ansatz.num_qubits, simulator=False
# )

# Or use a specific backend
backend = service.backend("ibm_brisbane")

# Initialize some variables to save the results from different runs:

exp_values_with_em0_es = []
exp_values_with_em1_es = []
exp_values_with_em2_es = []

# Use a pass manager to optimize the circuit and observables for the backend chosen:

pm = generate_preset_pass_manager(backend=backend, optimization_level=2)
isa_circuit = pm.run(circuit)
isa_observables = observables.apply_layout(layout=isa_circuit.layout)

# Open a session and run with no error mitigation:

estimator_options = EstimatorOptions(resilience_level=0, default_shots=10_000)

with Session(backend=backend) as session:
    estimator = Estimator(mode=session, options=estimator_options)

    pub = (isa_circuit, isa_observables, [individual_phases])
    cost = estimator.run([pub]).result()[0].data.evs

session.close()

exp_values_with_em0_es = cost[0]

# Open a session and run with resilience = 1:

estimator_options = EstimatorOptions(resilience_level=1, default_shots=10_000)

with Session(backend=backend) as session:
    estimator = Estimator(mode=session, options=estimator_options)

    pub = (isa_circuit, isa_observables, [individual_phases])
    cost = estimator.run([pub]).result()[0].data.evs

session.close()

exp_values_with_em1_es = cost[0]

# Open a session and run with resilience = 2:

estimator_options = EstimatorOptions(resilience_level=2, default_shots=10_000)

with Session(backend=backend) as session:
    estimator = Estimator(mode=session, options=estimator_options)

    pub = (isa_circuit, isa_observables, [individual_phases])
    cost = estimator.run([pub]).result()[0].data.evs

session.close()

exp_values_with_em2_es = cost[0]

前回と同様に、3段階の誤差軽減を行った場合の期待値を位相角の関数としてプロットすることができる。非常に難しいことだが、エラーを軽減することで結果がわずかに改善されることがわかる。この効果も、より深く複雑な回路ではより顕著になる。

import matplotlib.pyplot as plt

plt.plot(phases, exp_values_with_em0_es, "o", label="unmitigated")
plt.plot(phases, exp_values_with_em1_es, "o", label="resil = 1")
plt.plot(phases, exp_values_with_em2_es, "o", label="resil = 2")
plt.plot(phases, 2 * np.sin(phases / 2) ** 2 - 1, label="ideal")
plt.ylabel("Expectation")
plt.legend()
plt.show()

Output:

サマリー

このレッスンでは、コスト関数の作成方法を学びました：

コスト関数を作成する
IBM Quantum のプリミティブを活用してノイズを低減・抑制する方法
スピードと精度を最適化する測定戦略の立て方

これがハイレベルな変分ワークロードだ：

参照状態と変分状態を準備するユニタリーを備えた量子回路を示す図。これらはコスト関数の評価に使われる。

私たちのコスト関数は、最適化ループの各反復中に実行される。次のレッスンでは、古典的なオプティマイザが、新しいパラメータを選択するために、コスト関数の評価をどのように利用するかを探ります。

import qiskit
import qiskit_ibm_runtime

print(qiskit.version.get_version_info())
print(qiskit_ibm_runtime.version.get_version_info())

Output:

1.1.0
0.23.0

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

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