コスト関数

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

• まず、IBM Quantumプリミティブについて学びます。

• コスト関数 $C(\vec\theta)$ を定義します。これはオプティマイザーが最小化（または最大化）すべき問題の目標を定義する、問題固有の関数です。

• IBM Quantum プリミティブを使った測定戦略を定義し、速度と精度のトレードオフを最適化します。

 

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 進数の整数）です。

IBM Quantum 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プリミティブは、量子状態$|\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は、測定方法がわからない任意のオブザーバブルを、パウリ演算子と呼ばれる、より単純で測定可能なオブザーバブルに分解します。 任意の演算子は、$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}$を計算します。次に、各出力$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$がゼロである場合（すなわち密ではなく疎なパウリ分解の場合）にのみ効率的になり得ます。形式的には、この計算が効率的に解けるためには、非ゼロ項の数がキュビット数$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 — installs the packages this notebook needs if they are missing
import importlib.util

_needed = {"matplotlib": "matplotlib", "numpy": "numpy", "qiskit": "qiskit", "qiskit_aer": "qiskit-aer", "qiskit_ibm_runtime": "qiskit-ibm-runtime", "rustworkx": "rustworkx"}
_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")

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}$$

Since we do not know the eigenvalues or eigenstates of the target observable $\hat{H}$, first we need to consider its diagonalization. Given that $\hat{H}$ is Hermitian, there exists a unitary transformation $V$ such that $\hat{H}=V^\dagger \Lambda V$, where $\Lambda$ is the diagonal eigenvalue matrix, so $\langle j | \Lambda | k \rangle = 0$ if $j\neq k$, and $\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$ が既知です。これが 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 つの可能性が生じます。

• 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に関する手順は、次の例でさらに詳しく説明します。変分アルゴリズムにおいて特に重要な手順の1つは、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 つの計算で得られた値は非常に近いことに注意してください。結果を改善するための手法については、以下でさらに詳しく説明します。

非物理系へのマッピングの例

最大カット（マックスカット）問題は、グラフの頂点を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"
)

この問題は二値最適化問題として表現できます。グラフの各ノード$0 \leq i < n$（$n$はグラフのノード数であり、この場合$n=4$です）に対して、二値変数$x_i$を考えます。この変数は、ノード$i$がグループ$1$と呼ぶ方に属する場合は値$1$を、グループ$0$と呼ぶもう一方に属する場合は値$0$をとります。また、隣接行列$w$の要素$(i,j)$である$w_{ij}$を、ノード$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$。演算子が構築されたら、Qiskit回路ライブラリのQAOAAnsatz回路を使用することで、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

IBM Quantum 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 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)

1.1120776913677988

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

• オプティマイザーを活用して最適なパラメーターを見つける

• 最適パラメーターをアンザッツにバインドして固有値を求める

• 固有値を問題の定義に変換する

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

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

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

エラー抑制

エラー抑制とは、エラーを最小化するためにコンパイル時に 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を使用して実際のバックエンド（この場合は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)

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

この場合、これを2回行うことに注意してください。1回はoptimization_level = 0で、もう1回はそれを3に設定して行います。それぞれの場合において、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 に上げることで確実に向上したことがわかります。より深く複雑な 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 のエラーを低減できる技術を指します。一般的に、モデルの学習に関連する量子的な前処理オーバーヘッドと、生成したモデルを用いて生の結果のエラーを緩和するための古典的な後処理オーバーヘッドが発生します。

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

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

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

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}$$

緩和・抑制のためのQuantum Computeのオプションの使い方

以下に、IBM Quantum Compute Service でエラー緩和と抑制を使用しながら期待値を計算する方法を示します。先ほどとまったく同じ 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()

まとめ

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

• コスト関数の作成方法

• IBM Quantumプリミティブを活用してノイズを緩和・抑制する方法

• 速度と精度を最適化するための測定戦略の定義方法

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

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

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