TREXを用いたPEC

確率的エラーキャンセレーション(PEC)は、観測量の期待値計算において量子ノイズを軽減する手法です。期待値の不偏推定を計算できます。この手法は、学習したノイズモデルの逆を注入したノイズのある回路のアンサンブルを実行することで機能します。ノイズモデルは、Pauliツワーリング(ランダム化コンパイルとしても知られています)を用いてPauliチャネルに近似されます。このようなノイズモデルは効率的に学習でき、PECを用いて軽減できます(詳細はProbabilistic error cancellation with sparse Pauli-Lindblad models on noisy quantum processorsを参照してください)。

Twirled Readout Error eXtinction (TREX)は、測定エラーを軽減する手法で、Pauli観測量の期待値の推定に役立ちます。これは測定ツワーリング(測定の前にランダムに X ゲートを注入し、測定結果を反転させること)に依拠しており、読み出しエラーの遷移行列を対角化して、逆行列を求めやすくします。再スケーリング項は、ゼロ状態で初期化されたツワール化済みのキャリブレーション回路を測定することで、対角ノイズチャネルから学習されます。

このガイドでは、PECを使ってゲートノイズを軽減し、それをTREXと組み合わせて、観測量の期待値を推定する際の測定ノイズも軽減する方法を説明します。

ワークフローは次のとおりです。

1. 回路と観測量を指定します。
2. Qiskitのステートベクトルシミュレーターを使って、正確な期待値を計算します。
3. ノイズのある未軽減の期待値を計算します。
4. qiskit-noise-learning を使ってエンタングリングレイヤーのノイズを学習します。
5. PEC + TREXプログラムを準備し、TREXキャリブレーション回路を追加します。
6. 軽減された期待値を比較します。

回路と観測量を指定する

8量子ビット上の線形な横磁場イジング模型のトロッター化された時間発展のミラー回路を準備します。回路をミラーリングして回路のユニタリを打ち消し、恒等回路を作ります。これにより、各観測量の理想的な期待値を簡単に計算できます。さらに、ミラーリングは回路の深さを2倍にし、その結果ノイズの影響も2倍になります。

角度 $\theta= \frac{\pi}{2}$ を使って、回路内のすべてのゲートをクリフォード群に属するものに保ち、stabilizer シミュレーターによるより高速なシミュレーションを可能にします。

# 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", "qiskit_mitigation": "qiskit-mitigation", "qiskit_noise_learning": "qiskit-noise-learning", "samplomatic": "samplomatic"}
_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")

import warnings

import numpy as np
from qiskit.circuit import QuantumCircuit
from qiskit.quantum_info import SparsePauliOp

warnings.filterwarnings("ignore")

def construct_ising_circuit(
    num_qubits: int,
    num_trotter_steps: int,
    rx_angle: float,
    barrier: bool = True,
) -> QuantumCircuit:
    circuit = QuantumCircuit(num_qubits)

    for _step in range(num_trotter_steps):
        circuit.rx(rx_angle, range(num_qubits))
        if barrier:
            circuit.barrier()
        for first_qubit in (1, 2):
            for idx in range(first_qubit, num_qubits, 2):
                # equivalent to Rzz(-pi/2):
                circuit.sdg([idx - 1, idx])
                circuit.cz(idx - 1, idx)
        if barrier:
            circuit.barrier()

    return circuit

def mirror_circuit(circuit: QuantumCircuit) -> QuantumCircuit:
    mirror_circ = circuit.copy_empty_like()
    mirror_circ.compose(circuit, inplace=True)
    mirror_circ.barrier()
    mirror_circ.compose(circuit.inverse(), inplace=True)
    mirror_circ.measure_active()
    return mirror_circ

# define observables
num_qubits = 8
observables = [
    SparsePauliOp("I" * ((num_qubits - 1) - q) + "Z" + "I" * q)
    for q in range(num_qubits)
]

# Instantiate circuit
num_trotter_steps = 5
rx_angle = np.pi / 2

circuit = construct_ising_circuit(
    num_qubits, num_trotter_steps, rx_angle, barrier=False
)
mirrored_circuit = mirror_circuit(circuit)
print("Visualization of 1 mirrored Trotter step")
mirror_circuit(
    construct_ising_circuit(num_qubits, 1, rx_angle, barrier=False)
).draw("mpl", scale=0.7)

Visualization of 1 mirrored Trotter step

ノイズのある未軽減の期待値を計算する

QuantumProgramには、量子回路とsamplexが、Executorプリミティブに適した形式で含まれます。MitigationTask.prepare は、回路、観測量、パラメーター値から QuantumProgram を構築し、指定されたボクシングオプションで回路をボックス化します。measure_annotations ボクシングオプションを "change_basis" に設定して測定を必要な基底に回転させ、シミュレーターがノイズモデルを回路のボックスにマッピングできるように add_tags: "unique_box" を追加します。プログラムはデフォルトでPauliツワーリングされ、ハードウェア上ではコヒーレントノイズが、軽減手法が想定するPauliノイズに変換されます。num_randomizations は各回路のツワール化されたインスタンスをいくつサンプリングするかを設定し、shots_per_randomization は各インスタンスが受け取るショット数を設定します。したがって、その積が回路あたりのショット上限になります。

QPUにアクセスできる場合は、QiskitRuntimeService のバックエンドから構築した Executor でプログラムを実行できます。この例では AerSimulator を使ってプログラムをローカルで実行します。そのシミュレーターオプション layer_noise_model は、エンタングリングゲートのすべてのレイヤーにランダムなPauli-Lindbladノイズモデルを注入します。このモデルは qiskit-mitigation の random_pauli_lindblad_noise で生成し、レイヤーゲートあたりのエラー(EPLG)を現在のHeron QPUに典型的な 0.003 としています。

import matplotlib.pyplot as plt
from qiskit.circuit import CircuitInstruction
from qiskit.quantum_info import PauliLindbladMap, pauli_basis
from qiskit.transpiler import PassManager
from qiskit_aer import AerSimulator
from qiskit_ibm_runtime import Executor
from qiskit_mitigation import MitigationTask, find_combined_unique_layers
from qiskit_noise_learning.aer_executor.insert_noise_pass import (
    InsertNoisePass,
)
from samplomatic import Tag
from samplomatic.utils import get_annotation

def generate_random_pauli_lindblad_noise(
    layers,
    num_qubits: int,
    eplg: float = 3e-3,
    readout_noise_scale: float = 3e-2,
    seed: int | None = None,
) -> tuple[
    list[tuple[CircuitInstruction, PauliLindbladMap]], PauliLindbladMap
]:
    """Generate random Pauli-Lindblad noise over the full Pauli basis."""
    layer_noise_model = []
    for layer in layers:
        layer_instructions = layer.operation.params[0]
        contain_measurement = any(
            inst.name.startswith("measure") for inst in layer_instructions
        )
        if not contain_measurement:
            edges = [
                (inst.qubits[0]._index, inst.qubits[1]._index)
                for inst in layer_instructions
            ]  # extract the edges of the layer

            basis_paulis = [p for p in pauli_basis(2) if np.sum(p.x + p.z)]
            rng = np.random.default_rng(seed=seed)
            num_generators = len(basis_paulis)
            rates = rng.random(len(basis_paulis)) * eplg / num_generators

            layer_noise = PauliLindbladMap.from_sparse_list(
                [
                    (pauli.to_label(), edge, rate)
                    for pauli, rate in zip(basis_paulis, rates, strict=True)
                    for edge in edges
                ],
                num_qubits=num_qubits,
            )
            layer_noise_model.append((layer, layer_noise))
    # generate a readout noise model
    rng = np.random.default_rng(seed=seed)
    rates = rng.random(num_qubits) * readout_noise_scale
    readout_noise_model = PauliLindbladMap.from_sparse_list(
        [("X", [q], rate) for q, rate in enumerate(rates)],
        num_qubits=num_qubits,
    )
    return layer_noise_model, readout_noise_model

def add_readout_noise(program, readout_noise_model):
    """Injects readout noise to the circuits in a quantum program."""
    readout_noise_map = {}
    for item in program.items:
        readout_tag = item.circuit[-1].label[
            -4:
        ]  # extract the readout layer tag from the circuit
        readout_noise_map[readout_tag] = readout_noise_model
        item.circuit = PassManager(
            [InsertNoisePass(noise_dict=readout_noise_map, noise_after=False)]
        ).run(item.circuit)

num_randomizations = 300
shots_per_randomization = 64

boxing_options = {"add_tags": "unique_box"}

# find the unique layers to generate the layered noise for
box_layers = find_combined_unique_layers(
    [mirrored_circuit], custom_boxing_options=boxing_options
)
layer_noise_model, readout_noise_model = generate_random_pauli_lindblad_noise(
    box_layers, num_qubits, seed=42
)

# a local noisy simulator in place of a QPU backend
noisy_executor = Executor(
    AerSimulator(method="stabilizer"),
    options={"simulator": {"layer_noise_model": layer_noise_model}},
)

task = MitigationTask()
program = task.prepare(
    circuit=mirrored_circuit,
    observables=observables,
    parameters=None,
    num_randomizations=num_randomizations,
    shots_per_randomization=shots_per_randomization,
    custom_boxing_options=boxing_options,
)
# For the simulation - we inject readout noise to the circuits of the program
add_readout_noise(program, readout_noise_model)

program._semantic_role = "simulation"  # flag the program for the simulator
unmitigated_res = task.postprocess(noisy_executor.run(program).result())

x_vals = np.array(range(num_qubits))
fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(x_vals, [1] * len(x_vals), "k-", alpha=0.5, label="Exact")
ax.errorbar(
    x_vals,
    unmitigated_res.data.evs,
    unmitigated_res.data.stds,
    fmt="o",
    label="Unmitigated",
)

x_tick_labels = [r"$Z_{ {i} }$" for i in range(num_qubits)]
ax.set_xlabel(r"$\theta/\pi$")
ax.set_ylabel("Expectation value")
ax.legend()
ax.set_xticks(x_vals, x_tick_labels)
plt.tight_layout()
plt.show()

エンタングリングレイヤーのノイズを学習する

PECは、回路内の各エンタングリングゲートのレイヤーについて、ツワール化されたノイズモデルを学習することから始まります。マルコフ型のノイズモデルを仮定する場合、回路の深さは特定のエンタングリングレイヤーのノイズチャネルに影響しないと想定されます。したがって、一意なエンタングリングゲートのレイヤーごとにノイズモデルを学習するだけで済みます。このガイドで使用する回路は、2つの一意なレイヤーに分割できます。一意なレイヤーは個別に学習する必要があるため、ノイズモデルを学習したい場合は、一意なエンタングリングレイヤーの集合が最小になるように回路を設計することが重要です。回路内の一意なレイヤーが多すぎると、多数のレイヤーを学習するためのサンプリングコストのせいで、ノイズ学習が実行不可能になる場合があります。

ノイズの学習には qiskit-noise-learning を使用します。学習器は学習用回路内で各レイヤーを fragment_depths 回繰り返し、各回路を num_randomizations 回ツワール化し、測定された忠実度の減衰をPauli-Lindbladモデルにフィットさせます。これらのパラメーターは、学習したモデルの精度に影響を与えることがあります。学習したノイズモデルが正確であるほど、ノイズキャンセレーションもより正確になります。 次のコードは、2つのレイヤーを1つのマップにまとめた学習済みノイズモデルを描画します。各量子ビットは、2つのレイヤーにわたって合計した重み1のX、Y、Zレートの円グラフで表され、各エッジは学習した重み2のレートを示します。たとえば、3e-4 というレートは、そのレイヤーを10,000回使用するうちおよそ3回その誤りが発生することを意味します。

読み出しノイズの軽減は、TREXを使って後処理で行われます。これはゲートノイズモデルとは対照的です。ゲートノイズモデルは、学習したノイズモデルに対応する反ノイズを注入するために、プログラムの準備時に必要となります。 そのため、読み出しノイズは、ゲートノイズを学習するときに実行して学習することも、軽減されたプログラムを実行するときに学習することもできます。読み出しノイズをゲートノイズと一緒に学習すると、自己無撞着なノイズモデルを学習するうえで有益な場合があります(たとえばDisambiguating Pauli noise in quantum computersを参照してください)。ゲートノイズと読み出しノイズがノイズモデル内で分離されている場合は、ノイズが学習されたときと実行されたときの間のノイズのドリフトを避けるため、軽減されたプログラムを実行するときに読み出しノイズを学習することをお勧めします。

このガイドでは、TREXキャリブレーション回路を使って、軽減されたプログラムを実行しながら読み出しノイズを学習します。

from qiskit.providers.fake_provider import GenericBackendV2
from qiskit.quantum_info import PauliList
from qiskit_ibm_runtime.results import LayerError, PauliLindbladError
from qiskit_ibm_runtime.visualization.embeddings import Embedding
from qiskit_mitigation import PEC
from qiskit_noise_learning.aer_executor import AerExecutor
from qiskit_noise_learning.noise_learner import LearningOptions, NoiseLearner
from samplomatic import InjectNoise

line_coupling = [(i, i + 1) for i in range(num_qubits - 1)]
backend = GenericBackendV2(
    num_qubits=num_qubits,
    basis_gates=["id", "rz", "sx", "x", "cz"],
    coupling_map=line_coupling,
)

task_pec = PEC()
# first learn only the layers PEC will inject noise into, i.e. not the measurement layer
pec_layers = [
    layer
    for layer in task_pec.find_unique_layers(mirrored_circuit, boxing_options)
    if get_annotation(layer.operation, InjectNoise)
]
layers = task_pec.find_unique_layers(mirrored_circuit, boxing_options)

learning_options = LearningOptions(
    fragment_depths=[2, 4, 8, 16],
    num_randomizations=250,
    shots_per_randomizations=16,
)

# qiskit-noise-learning is not compatible with the Executor in local mode,
# and requires the simulator present in the package for running simulation
noise_dict = {
    get_annotation(layer.operation, Tag).ref: noise_layer
    for layer, noise_layer in layer_noise_model
}
noise_learner_simulator = AerExecutor(
    qasm_simulator=AerSimulator(method="stabilizer"), noise_dict=noise_dict
)
learner = NoiseLearner(
    backend, options=learning_options, executor=noise_learner_simulator
)
learned_noise = (
    learner.run(pec_layers).result().to_dict()
)  # keyed by InjectNoise reference

# Draw the learned model on the line, laid out as a line with the qubits in order.
# Composing the two layers' maps sums their rates, giving one map over all the edges.
embedding = Embedding(
    [(0, 0), (0, 2), (0, 4), (2, 4), (4, 4), (4, 2), (4, 0), (2, 0)],
    line_coupling,
)
first, second = learned_noise.values()
learned_model = first.compose(second).simplify()
error = PauliLindbladError(
    PauliList(
        learned_model.get_qubit_sparse_pauli_list_copy().to_pauli_list()
    ),
    learned_model.rates,
)
fig = LayerError(circuit, list(range(num_qubits)), error).draw_map(
    embedding, radius=0.4, height=460, width=620
)
fig.update_traces(
    marker={
        "cmin": 0,
        "cmax": max(learned_model.rates),
        "colorbar": {"tickformat": ".1e", "title": "rate"},
    },
    selector={"mode": "markers"},
)
fig.show()

PECとTREXを組み合わせたプログラムを準備する

ここでは、TREX タスクと PEC タスクの両方を準備し、TREX タスクを PEC タスクの変数として追加します。これにより、TREXによる軽減に必要なすべてのプロパティ(測定ツワーリングなど)が、PECで軽減される項目に確実に追加されます。注入されるノイズはランダム化ごとに1回サンプリングされるため、PECの統計誤差はランダム化あたりのショット数ではなくランダム化の回数で制御されます。そのため、それぞれ少数のショットで多数のランダム化を使用し、総ショット数を未軽減のプログラムと同じに保ちます。

プログラムの実行中に読み出しノイズを学習するために、TREXキャリブレーション回路がプログラムに追加されます。TREXキャリブレーション回路は非常に浅いため、そのための追加のランダム化は比較的低コストで行え、読み出しノイズをより正確に学習できます。

from qiskit_mitigation import TREX

trex = TREX()

program_pec = task_pec.prepare(
    circuit=mirrored_circuit,
    observables=observables,
    parameters=None,
    num_randomizations=1200,
    shots_per_randomization=16,
    custom_boxing_options=boxing_options,
    trex=trex,
    noise_maps=learned_noise,  # PEC option
    noise_gain=0,  # PEC option
    max_sampling_overhead=100,  # PEC option
)
program_pec = trex.prepare(
    num_randomizations=1500,
    quantum_program=program_pec,
    custom_boxing_options=boxing_options,
)

# For the simulation - we inject readout noise to the circuits of the program
add_readout_noise(program_pec, readout_noise_model)

ショットオーバーヘッドを確認し、プログラムを実行して後処理する

PECは不偏な期待値の推定を返しますが、その代償として大きなショットオーバーヘッドが生じます。必要なショットオーバーヘッドは、学習したノイズモデルが回路に与える影響を反映する回路ガンマ係数を計算することで近似できます。デフォルトでは、PECタスクのランダム化の回数は、不偏推定を得るために必要なショットオーバーヘッドのおよその量(~$\gamma^2$)でスケールされます。この自動スケーリングを避けるには、パラメーター scale_randomizations_by_gamma を False に設定できます。ショットオーバーヘッドは max_sampling_overhead パラメーターで上限が設定されます。

PECが不偏推定に到達するために必要なショット数は、ノイズレートと回路の深さに対して指数的に増加するため、実行前にその数を確認し、この量が許容できるかどうかを判断する必要があります。

print(f"The circuit calculated gamma is: {task_pec.gamma}")
print(f"The circuit shot overhead is ~{task_pec.gamma**2}")

The circuit calculated gamma is: 1.2774846896089953
The circuit shot overhead is ~1.631967132185391

program_pec._semantic_role = (
    "simulation"  # flag the program for the simulator
)
job_pec_results = noisy_executor.run(program_pec).result()

# post-processing
pec_res = task_pec.postprocess(job_pec_results)
# In order to post process the results without TREX, we call `postprocess` again, but instead of the default options we
# pass only the PEC item result and the `measure_noise_data` as `None`
pec_without_trex_res = task_pec.postprocess(
    job_pec_results[task_pec._program_item_index], measure_noise_data=None
)

軽減された期待値を比較する

最後に、各観測量の期待値をプロットします。PEC + TREXは誤差範囲内で正確な期待値を回復しますが、未軽減の結果は大きく下方にずれています。

x = np.array(range(num_qubits))
x_tick_labels = [r"$Z_{ {i} }$" for i in range(num_qubits)]

fig, ax = plt.subplots(figsize=(8, 5))
ax.errorbar(
    x,
    unmitigated_res.data.evs,
    unmitigated_res.data.stds,
    fmt="o",
    label="Twirl only",
)
ax.errorbar(
    x,
    pec_without_trex_res.data.evs,
    pec_without_trex_res.data.stds,
    fmt="o",
    label="PEC only",
)
ax.errorbar(
    x,
    pec_res.data.evs,
    pec_res.data.stds,
    fmt="o",
    label="PEC + TREX",
)

ax.plot(x, [1] * len(x), "--", color="black", alpha=0.5, label="Exact")

tick_labels = [f"$Z_{ {i} }$" for i in range(num_qubits)]
ax.set_xticks(x, tick_labels)

ax.set_ylabel("Expectation value")
ax.legend()
ax.set_title("PEC+TREX Mitigated Expectation Values")
plt.tight_layout()
plt.show()