コスト関数

このレッスンでは、コスト関数の評価方法を学びます。

• まず、Qiskit Runtime primitives について学びます。
• コスト関数 $C(\vec\theta)$ を定義します。これはオプティマイザーが最小化（または最大化）すべき問題の目標を定義する、問題固有の関数です。
• Qiskit Runtime primitives を使った測定戦略を定義し、速度と精度のトレードオフを最適化します。

 

Primitives

古典・量子を問わず、あらゆる物理系はさまざまな状態を取り得ます。たとえば、道路上の車は質量・位置・速度・加速度といった物理量によってその状態が特徴付けられます。同様に、量子系も異なる配置や状態を持ちますが、測定や状態の時間発展の扱い方が古典系とは異なります。これにより、量子力学に固有の重ね合わせやエンタングルメントといった特性が生じます。車の状態を速度や加速度などの物理量で記述できるように、量子系の状態も数学的オブジェクトであるオブザーバブルを用いて記述できます。

量子力学では、状態は規格化された複素数の列ベクトル、すなわちケット（$|\psi\rangle$）で表され、オブザーバブルはケットに作用するエルミート線形演算子（$\hat{H}=\hat{H}^{\dagger}$）で表されます。オブザーバブルの固有ベクトル（$|\lambda\rangle$）は固有状態と呼ばれます。オブザーバブルをその固有状態のひとつ（$|\lambda\rangle$）に対して測定すると、対応する固有値（$\lambda$）が読み出しとして得られます。

量子系をどのように測定するか、そして何を測定できるかについては、Qiskit が提供する 2 つの primitive が役立ちます。

• Sampler：量子状態 $|\psi\rangle$ が与えられた場合、この primitive は各計算基底状態が得られる確率を求めます。
• Estimator：量子オブザーバブル $\hat{H}$ と状態 $|\psi\rangle$ が与えられた場合、この primitive は $\hat{H}$ の期待値を計算します。

Sampler primitive

Sampler primitive は、状態 $|\psi\rangle$ を準備する量子 Circuit が与えられたとき、計算基底の各状態 $|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$ は Qubit の数、$k$ は任意の出力バイナリ文字列 $\{0,1\}^n$ の整数表現（すなわち 2 進数の整数）です。

Qiskit Runtime の Sampler は、量子デバイス上で Circuit を複数回実行し、各実行で測定を行い、得られたビット列から確率分布を復元します。実行回数（ショット数）が多いほど結果は正確になりますが、より多くの時間と量子リソースが必要です。

ただし、可能な出力の数は Qubit 数 $n$ に対して指数的に増加する（すなわち $2^n$）ため、密な確率分布を捉えるにはショット数も指数的に増やす必要があります。そのため、Sampler が効率的に機能するのは疎な確率分布の場合に限られます。ここで、目標状態 $|\psi\rangle$ は計算基底状態の線形結合として表現可能であり、その項数が Qubit 数に対して高々多項式的に増加することが条件です。

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

また Sampler は、Circuit の一部分（全状態の部分集合）から確率を取得するように設定することもできます。

Estimator primitive

Estimator primitive は、量子状態 $|\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$$

しかし、オブザーバブルの期待値を計算することは常に可能とは限りません。固有基底が未知の場合が多いからです。Qiskit Runtime の Estimator は、複雑な代数的処理を用いて、実際の量子デバイス上でオブザーバブルを固有基底が既知の他のオブザーバブルの組み合わせに分解することで期待値を推定します。

簡単に言えば、Estimator は測定方法が不明なオブザーバブルを、パウリ演算子と呼ばれる、より単純で測定可能なオブザーバブルに分解します。任意の演算子は $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$ は Qubit 数、$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 は元の Circuit から各オブザーバブル $\hat{P}_k$ に対して新しい Circuit $V_k|\psi\rangle$ を導出し、パウリオブザーバブルを計算基底で効果的に対角化して測定します。パウリオブザーバブルは $V_k$ が事前に分かっているため容易に測定できますが、一般のオブザーバブルではそうではありません。

各 $\hat{P}_{k}$ について、Estimator は対応する Circuit を量子デバイス上で複数回実行し、計算基底で出力状態を測定し、各可能な出力 $j$ が得られる確率 $p_{kj}$ を計算します。次に各出力 $j$ に対応する $P_k$ の固有値 $\lambda_{kj}$ を求め、$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$ 個のパウリの期待値を計算することは（指数的に増加するため）非現実的であり、Estimator が効率的に機能するのは $w_k$ の大部分がゼロの場合（すなわち密な分解ではなく疎なパウリ分解）に限られます。形式的には、この計算が効率的に解けるためには、非ゼロ項の数が Qubit 数 $n$ に対して高々多項式的に増加することが必要です：$\hat{H} = \sum^{\text{Poly}(n)}_k w_k \hat{P}_k.$

Sampler で説明したように、確率のサンプリングも効率的である必要があるという暗黙の仮定があることに読者は気づくでしょう。すなわち

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

期待値計算のガイド付き例

単一 Qubit の状態 $|+\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$ は可換でない（すなわち、固有基底が共有されていない）ため、同時に測定することはできません。そのため、補助 Circuit が必要です。

# Added by doQumentation — required packages for this notebook
!pip install -q matplotlib numpy qiskit qiskit-aer qiskit-ibm-runtime rustworkx

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

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

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

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

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

数学的厳密性（任意）

$|\psi\rangle$ を $\hat{H}$ の固有状態の基底で展開すると $|\psi\rangle = \sum_\lambda a_\lambda |\lambda\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}$ がエルミートであるため、ユニタリ変換 $V$ が存在し、$\hat{H}=V^\dagger \Lambda 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}$$

系が状態 $|\phi\rangle = V |\psi\rangle$ にある場合に $| j\rangle$ を測定する確率が $p_j = |\langle j|\phi \rangle|^2$ であることから、上記の期待値は次のように表せます。

$$\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$ が既知です。これが Runtime の 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 つの可能性が生じます。

• Ansatz が探索空間上の解状態を定義していない場合：この場合、オプティマイザーは解を見つけられないため、探索空間をより正確に表現できる別の Ansatz を試す必要があります。
• オプティマイザーが有効な解を見つけられない場合：最適化はグローバルに定義される場合とローカルに定義される場合があります。これについては後のセクションで詳しく説明します。

いずれにせよ、古典的な最適化ループを実行しながら、コスト関数の評価を量子コンピューターに頼ることになります。この観点から見ると、オプティマイザーがコスト関数を評価するたびにブラックボックス量子オラクルを呼び出す純粋に古典的な取り組みと考えることもできます。

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

まずシミュレーターである StatevectorEstimator を使って実行します。デバッグには通常これが推奨されますが、デバッグ実行の直後に実際の量子ハードウェアでの計算も行います。関心のある問題の多くは、最先端のスーパーコンピューター施設なしには古典的にシミュレートできなくなってきています。

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

[-0.58744589]

次に、実際の量子コンピューターで実行します。構文の変更に注意してください。pass_manager に関連するステップは次の例でさらに詳しく説明します。変分アルゴリズムにおいて特に重要なステップのひとつが、Qiskit Runtime セッションの使用です。セッションを開始することで、パラメーターが更新されるたびに新しいキューで待機することなく、変分アルゴリズムの複数回のイテレーションを実行できます。これは、キュー待ち時間が長い場合や多くのイテレーションが必要な場合に重要です。IBM Quantum® Network のパートナーのみが Runtime セッションを使用できます。セッションへのアクセス権がない場合は、一度に送信するイテレーション数を減らし、最新のパラメーターを保存して将来の実行に使用することができます。イテレーションを送信しすぎたり、キュー待ち時間が長すぎる場合は、ジョブ送信間の長い遅延を示すエラーコード 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 Runtime 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つの集合の間の辺の数を最大化する組み合わせ最適化問題です。より正式には、$V$ を頂点の集合、$E$ を辺の集合とする無向グラフ $G=(V,E)$ が与えられたとき、Max-Cut 問題は頂点を2つの互いに素な部分集合 $S$ と $T$ に分割し、一方の端点が $S$、もう一方が $T$ にある辺の数を最大化することを求めます。

Max-Cut はクラスタリング、ネットワーク設計、相転移など、さまざまな問題を解くために応用できます。まず問題グラフを作成してみましょう：

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

この問題は2値最適化問題として表現できます。グラフのノード数を $n$（この場合は $n=4$）として、各ノード $0 \leq i < n$ に対して2値変数 $x_i$ を考えます。この変数は、ノード $i$ がグループ $1$ と呼ぶ集合に属する場合は $1$、グループ $0$ と呼ぶ集合に属する場合は $0$ の値をとります。また、$w_{ij}$（隣接行列 $w$ の要素 $(i,j)$）をノード $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$ は Qubit $i$ に作用する $Z$ Gate と等価です。さらに：

$$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番目の Qubit が最も右に位置することに注意してください）：$IIZZ + IZIZ + IZZI + ZIIZ + ZZII$。演算子が構築されたら、Qiskit Circuit ライブラリの QAOAAnsatz Circuit を使って QAOA アルゴリズム用のアンザッツを簡単に構築できます。

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

# Sum the weights, and divide by 2

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

Offset: -2.5

Runtime Estimator はハミルトニアンとパラメーター化されたアンザッツを直接受け取り、必要なエネルギーを返します。そのため、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)

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

1.1120776913677988

このサンプルは「Applications」のセクションで改めて取り上げ、オプティマイザーを活用して探索空間を反復する方法を学びます。一般的に、このプロセスには次のステップが含まれます：

• オプティマイザーを活用して最適なパラメーターを見つける
• 最適パラメーターをアンザッツにバインドして固有値を求める
• 固有値を問題の定義に変換する

測定戦略：速度と精度のトレードオフ

前述のとおり、私たちはノイズのある量子コンピューターをブラックボックスオラクルとして使用しています。ノイズにより取得される値が非決定的になり、ランダムな揺らぎが生じます。この揺らぎは、特定のオプティマイザーが提案された解に収束することを妨げたり、完全に防いだりする可能性があります。これは、量子ユーティリティの探索を段階的に進め、量子優位性に向けて前進する上で対処すべき一般的な課題です：

Qiskit Runtime Primitive のエラー抑制およびエラー軽減オプションを使用することで、ノイズに対処し、現在の量子コンピューターの有用性を最大化できます。

エラー抑制

エラー抑制とは、エラーを最小化するためにコンパイル時に Circuit を最適化・変換する技術を指します。これは基本的なエラー処理技術であり、通常は全体の実行時間に対して古典的な前処理のオーバーヘッドが生じます。このオーバーヘッドには、量子ハードウェア上で実行するための Circuit の Transpile が含まれます：

• 量子システムで利用可能なネイティブゲートを使って Circuit を表現する
• 仮想 Qubit を物理 Qubit にマッピングする
• 接続性の要件に基づいて SWAP を追加する
• 1Q ゲートおよび 2Q ゲートを最適化する
• アイドル状態の Qubit にダイナミカルデカップリングを追加してデコヒーレンスの影響を防ぐ

Primitive では optimization_level オプションの設定と高度なトランスパイルオプションの選択により、エラー抑制技術を利用できます。後続のコースでは、結果を改善するためのさまざまな Circuit 構築方法を詳しく説明しますが、ほとんどの場合は optimization_level=3 に設定することを推奨します。

ここでは、単純な理想的な動作を持つサンプル Circuit を用いて、トランスパイルプロセスにおける最適化レベルを上げることの価値を可視化してみましょう。

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

上記の Circuit は、$[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 を用いて実際の Backend（この例では ibm_brisbane）を取得し、AerSimulator でその Backend のノイズ挙動を含めてシミュレーションします。

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)

次に、パスマネージャーを使って Circuit を Backend の「命令セットアーキテクチャ（ISA）」にトランスパイルします。これは Qiskit Runtime の新しい要件です。Backend に送信されるすべての Circuit は、Backend のターゲットの制約に準拠していなければなりません。つまり、Backend の ISA（デバイスが理解・実行できる命令セット）に基づいて記述される必要があります。これらのターゲット制約は、デバイスのネイティブな基底ゲート、Qubit の接続性、そして必要に応じてパルスやその他の命令タイミング仕様といった要素によって決まります。

なお、今回のケースではこれを2回行います。1回目は optimization_level = 0、2回目は 3 に設定します。それぞれで Estimator Primitive を使い、さまざまな位相の値における観測量の期待値を推定します。

# 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 に上げることで確実に向上したことがわかります。より深く複雑な Circuit では、最適化レベル 0 と 3 の差はより顕著になる可能性があります。ここで使用したのは、toy モデルとして用いたごくシンプルな Circuit です。

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

エラー緩和

エラー緩和とは、実行時にデバイスのノイズをモデル化することで、Circuit のエラーを低減できる技術を指します。一般的に、モデルの学習に関連する量子的な前処理オーバーヘッドと、生成したモデルを用いて生の結果のエラーを緩和するための古典的な後処理オーバーヘッドが発生します。

Qiskit Runtime プリミティブの resilience_level オプションは、エラーに対してどの程度の耐性を持たせるかを指定します。レベルが高いほど、量子サンプリングのオーバーヘッドによる処理時間の増加と引き換えに、より正確な結果が得られます。レジリエンスレベルを使用することで、プリミティブクエリにエラー緩和を適用する際のコストと精度のトレードオフを調整できます。

エラー緩和技術を実装する場合、緩和前のバイアスと比較して、結果のバイアスが低減されることが期待されます。場合によっては、バイアスが完全に消えることもあります。ただし、これにはコストが伴います。推定量のバイアスを低減するにつれて、統計的なばらつき（すなわち分散）が増加します。これはサンプリングプロセスにおける Circuit あたりのショット数をさらに増やすことで対処できます。ただし、バイアスを低減するために必要な量を超えたオーバーヘッドが発生するため、デフォルトでは行われません。以下の例に示すように、options.executions.shots で Circuit あたりのショット数を調整することで、この動作を簡単に有効にできます。

このコースでは、Qiskit Runtime プリミティブが実行できるエラー緩和について、完全な実装の詳細を必要とせずに概要レベルで説明するために、これらのエラー緩和モデルを取り上げます。

Twirled readout error extinction（T-REx）

Twirled readout error extinction（T-REx）は、Pauli twirling と呼ばれる技術を使用して、量子測定プロセス中に発生するノイズを低減します。この技術はノイズの特定の形式を前提としないため、非常に汎用的かつ効果的です。

全体的なワークフロー：

1. ゼロ状態のデータをランダムなビットフリップ（測定前の Pauli X）付きで取得する
2. 目的の（ノイズのある）状態のデータをランダムなビットフリップ（測定前の Pauli X）付きで取得する
3. 各データセットに対して特殊関数を計算し、除算する

 

これは options.resilience_level = 1 で設定でき、以下の例で示します。

ゼロノイズ外挿法

ゼロノイズ外挿法（ZNE）は、目的の量子状態を準備する Circuit のノイズを最初に増幅し、複数の異なるノイズレベルで測定値を取得したうえで、それらの測定値を用いてノイズのない結果を推定するという手法です。

全体的なワークフロー：

1. 複数のノイズ係数に対して Circuit ノイズを増幅する
2. ノイズを増幅したすべての Circuit を実行する
3. ゼロノイズの限界へ外挿する

 

これは 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}$$

Qiskit Runtime の緩和・抑制オプションの使い方

以下に、Qiskit Runtime でエラー緩和と抑制を使用しながら期待値を計算する方法を示します。先ほどとまったく同じ Circuit とオブザーバブルを利用できますが、今回は最適化レベルをレベル 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 つのエラー緩和レベルについて、位相角の関数として期待値をプロットできます。注意深く見ると、エラー緩和によって結果がわずかに改善されていることが分かります。この効果は、より深く複雑な Circuit においてはるかに顕著に現れます。

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

まとめ

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

• コスト関数の作成方法
• Qiskit Runtime プリミティブを活用してノイズを緩和・抑制する方法
• 速度と精度を最適化するための測定戦略の定義方法

変分ワークロードの全体像は以下のとおりです。

コスト関数は最適化ループの反復ごとに実行されます。次のレッスンでは、古典的なオプティマイザーがコスト関数の評価を使用して新しいパラメーターを選択する方法について説明します。

import qiskit
import qiskit_ibm_runtime

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

1.1.0
0.23.0