TREX를 사용한 PEC

확률적 오류 상쇄(PEC)는 관측량 기댓값 계산에서 양자 노이즈를 완화하는 기법입니다. 이를 통해 기댓값의 편향 없는 추정치를 계산할 수 있습니다. 이 기법은 학습된 노이즈 모델의 역을 주입한 노이즈 Circuit들의 앙상블을 실행하는 방식으로 동작합니다. 이후 Pauli twirling(randomized compiling이라고도 함)을 사용해 노이즈 모델을 Pauli 채널로 근사합니다. 이러한 노이즈 모델은 효율적으로 학습할 수 있으며 PEC를 사용해 완화할 수 있습니다(자세한 내용은 Probabilistic error cancellation with sparse Pauli-Lindblad models on noisy quantum processors 참조).

Twirled Readout Error eXtinction(TREX)는 측정 오류를 완화하는 기법으로, Pauli 관측량의 기댓값을 추정할 때 유용합니다. 이 기법은 측정 twirling(측정 전에 X Gate를 무작위로 주입하고 측정된 결과를 뒤집는 것)에 의존하며, 이는 판독 오류 전이 행렬을 대각화하여 쉽게 역변환할 수 있게 만듭니다. 영 상태로 초기화된 twirling된 보정 Circuit을 측정하여 대각 노이즈 채널로부터 재조정 항을 학습합니다.

이 가이드는 PEC를 사용해 Gate 노이즈를 완화하고, 관측량의 기댓값을 추정할 때 TREX와 결합해 측정 노이즈도 함께 완화하는 방법을 보여줍니다.

워크플로는 다음과 같습니다.

1. Circuit과 관측량을 지정합니다.
2. Qiskit의 상태 벡터 시뮬레이터를 사용해 정확한 기댓값을 계산합니다.
3. 노이즈가 있고 완화되지 않은 기댓값을 계산합니다.
4. qiskit-noise-learning을 사용해 얽힘 레이어의 노이즈를 학습합니다.
5. PEC + TREX 프로그램을 준비하고 TREX 보정 Circuit을 추가합니다.
6. 완화된 기댓값을 비교합니다.

Circuit과 관측량 지정하기

Qubit 8개에서 선형 가로장 Ising 모델의 Trotter화된 시간 발전에 대한 미러 Circuit을 준비합니다. Circuit을 미러링하면 Circuit 유니터리가 상쇄되어 항등 Circuit이 만들어집니다. 이를 통해 각 관측량의 이상적인 기댓값을 쉽게 계산할 수 있습니다. 또한 미러링은 Circuit 깊이를 두 배로 늘리고, 그에 따라 노이즈의 영향도 두 배가 됩니다.

각도 $\theta= \frac{\pi}{2}$를 사용해 Circuit의 모든 Gate를 Clifford 군에 속하게 하면, 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은 Executor primitive에 적합한 형식으로 양자 Circuit과 samplex를 담고 있습니다. MitigationTask.prepare는 Circuit, 관측량, 매개변수 값으로 QuantumProgram을 만들고, 주어진 박싱 옵션으로 Circuit을 박스로 묶습니다. measure_annotations 박싱 옵션을 "change_basis"로 설정하면 측정이 필요한 기저로 회전되며, 시뮬레이터가 노이즈 모델을 Circuit의 박스에 매핑할 수 있도록 add_tags: "unique_box"를 추가합니다. 프로그램은 기본적으로 Pauli twirling되며, 하드웨어에서는 이것이 결맞음 노이즈를 완화 방법이 가정하는 Pauli 노이즈로 바꿉니다. num_randomizations는 각 Circuit에서 twirling된 인스턴스를 몇 개 샘플링할지 정하고, shots_per_randomization은 각 인스턴스가 받는 shot 수를 정하므로, 둘의 곱이 Circuit당 shot 한도가 됩니다.

QPU에 접근할 수 있다면 QiskitRuntimeService의 Backend로 만든 Executor로 프로그램을 실행할 수 있습니다. 이 예제에서는 AerSimulator를 사용해 프로그램을 로컬에서 실행합니다. 시뮬레이터 옵션 layer_noise_model은 모든 얽힘 Gate 레이어에 무작위 Pauli-Lindblad 노이즈 모델을 주입합니다. 이 노이즈 모델은 qiskit-mitigation의 random_pauli_lindblad_noise로 생성하며, 레이어드 Gate당 오류(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는 Circuit에 있는 각 얽힘 Gate 레이어의 twirling된 노이즈 모델을 학습하는 것으로 시작합니다. 마르코프 노이즈 모델을 가정하면 Circuit의 깊이는 특정 얽힘 레이어의 노이즈 채널에 영향을 주지 않는다고 봅니다. 따라서 고유한 얽힘 Gate 레이어마다 노이즈 모델만 학습하면 됩니다. 이 가이드에서 사용하는 Circuit은 두 개의 고유한 레이어로 나눌 수 있습니다. 노이즈 모델을 학습하려는 경우, 고유한 레이어는 독립적으로 학습해야 하므로 고유한 얽힘 레이어의 집합이 최소가 되도록 Circuit을 설계하는 것이 중요합니다. Circuit에 고유한 레이어가 너무 많으면 많은 레이어를 학습하는 샘플링 비용 때문에 노이즈 학습이 불가능해질 수 있습니다.

노이즈 학습에는 qiskit-noise-learning을 사용합니다. 학습기는 학습용 Circuit에서 각 레이어를 _fragment_depths_번 반복하고, 각 Circuit을 _num_randomizations_번 twirling하며, 측정된 충실도의 감쇠를 Pauli-Lindblad 모델에 맞춥니다. 이 매개변수들은 학습된 모델의 정확도에 영향을 줄 수 있습니다. 학습된 노이즈 모델이 정확할수록 노이즈 상쇄도 더 정확해집니다. 다음 코드는 두 레이어를 하나의 맵으로 합쳐 학습된 노이즈 모델을 그립니다. 각 Qubit은 두 레이어에 걸쳐 합산한 가중치 1의 X, Y, Z 비율을 파이 차트로 나타내고, 각 엣지는 학습된 가중치 2의 비율을 보여줍니다. 예를 들어 3e-4라는 비율은 레이어를 10,000번 사용할 때 해당 오류가 대략 세 번 발생한다는 뜻입니다.

판독 노이즈 완화는 TREX를 사용해 후처리 단계에서 수행됩니다. 이는 학습된 노이즈 모델에 해당하는 반노이즈를 주입하기 위해 프로그램 준비 단계에서 필요한 Gate 노이즈 모델과 대조됩니다. 따라서 판독 노이즈는 Gate 노이즈를 학습할 때 함께 실행하여 학습할 수도 있고, 완화된 프로그램을 실행할 때 학습할 수도 있습니다. 판독 노이즈를 Gate 노이즈와 함께 학습하면 자기 일관적인 노이즈 모델을 학습하는 데 도움이 될 수 있습니다(예를 들어 Disambiguating Pauli noise in quantum computers 참조). 노이즈 모델에서 Gate 노이즈와 판독 노이즈가 분리되어 있는 경우에는, 노이즈를 학습한 시점과 실행한 시점 사이의 노이즈 드리프트를 피하기 위해 완화된 프로그램을 실행할 때 판독 노이즈를 학습하는 것을 권장합니다.

이 가이드에서는 TREX 보정 Circuit을 사용해 완화된 프로그램을 실행하는 동안 판독 노이즈를 학습합니다.

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 완화에 필요한 모든 속성(예: 측정 twirling)이 PEC로 완화되는 항목에 추가됩니다. 주입되는 노이즈는 무작위화마다 한 번씩 샘플링되므로 PEC의 통계적 오차는 무작위화당 shot 수가 아니라 무작위화 횟수로 조절됩니다. 따라서 무작위화는 많이, 각각의 shot은 적게 사용하여 전체 shot 수를 완화하지 않은 프로그램과 같게 유지합니다.

프로그램을 실행하는 동안 판독 노이즈를 학습하기 위해 TREX 보정 Circuit을 프로그램에 추가합니다. TREX 보정 Circuit은 매우 얕아서 이에 대한 추가 무작위화가 비교적 저렴하며, 판독 노이즈를 더 정확하게 학습할 수 있습니다.

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)

shot 오버헤드를 살펴보고, 프로그램을 실행한 뒤 후처리하기

PEC는 편향 없는 기댓값 추정치를 반환하지만 큰 shot 오버헤드가 따릅니다. 필요한 shot 오버헤드는 학습된 노이즈 모델이 Circuit에 미치는 영향을 반영하는 Circuit 감마 인자를 계산하여 근사할 수 있습니다. 기본적으로 PEC 작업의 무작위화 횟수는 편향 없는 추정에 필요한 대략적인 shot 오버헤드(~$\gamma^2$)만큼 조정됩니다. 이 자동 조정을 피하려면 매개변수 scale_randomizations_by_gamma를 False로 설정하면 됩니다. shot 오버헤드는 max_sampling_overhead 매개변수로 상한이 정해집니다.

이 값은 노이즈 비율과 Circuit 깊이에 따라 지수적으로 증가하므로, 실행하기 전에 PEC가 편향 없는 추정에 도달하는 데 필요한 shot 수를 검토하여 감당할 수 있는 수준인지 판단해야 합니다.

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