Extrapolace k nulovému šumu s hradlovým foldingem a vlastní extrapolací
Extrapolace k nulovému šumu (ZNE) je technika pro zmírňování kvantového šumu při výpočtech střední hodnoty pozorovatelných. Funguje tak, že spustí kvantový obvod při rostoucích úrovních šumu, vezme střední hodnoty vzhledem k těmto úrovním šumu a extrapoluje, aby odhadla střední hodnotu při nulovém šumu. Fáze zesilování šumu i extrapolace lze implementovat různými způsoby. Tento průvodce zesiluje šum pomocí „digitálního hradlového foldingu“, což znamená, že dvouqubitová hradla se nahradí ekvivalentními posloupnostmi hradla a jeho inverze. Například nahrazení unitárního operátoru za by přineslo faktor zesílení šumu 3. Zesilování můžeš implementovat také pomocí pravděpodobnostního zesilování chyb (PEA), které vyžaduje další krok naučení modelu šumu a je popsáno v průvodci pravděpodobnostní zesilování chyb. Tento průvodce extrapoluje pomocí několika vestavěných strategií fitování (lineární a exponenciální) a také ukazuje, jak implementovat vlastní fitovací funkci.
Pamatuj, že ačkoli ZNE může být účinná při zmírňování kvantového šumu za nižší cenu vzorkování než jiné metody, nemá žádné meze na zkreslení (bias).
Pracovní postup je následující:
- Zadej obvod a pozorovatelnou
- Vypočítej přesné střední hodnoty pomocí simulátoru stavového vektoru v Qiskitu
- Vypočítej zašuměné, nezmírněné střední hodnoty
- Připrav program ZNE s digitálním hradlovým foldingem a vzorkuj obvody s různými úrovněmi šumu
- Extrapoluj k limitě nulového šumu pomocí vestavěných i vlastních fitovacích funkcí
- Porovnej zmírněné střední hodnoty
Zadání obvodu a pozorovatelné
Vypočítej magnetizaci zprůměrovanou přes uzly trotterizovaného Isingova modelu s příčným polem (TFIM) na osmi qubitech v prstencové topologii. Projdeš několik hodnot síly příčného pole a vykreslíš magnetizaci jako funkci síly pole.
# Added by doQumentation — installs the packages this notebook needs if they are missing
import importlib.util
_needed = {"IPython": "IPython", "matplotlib": "matplotlib", "numpy": "numpy", "qiskit": "qiskit", "qiskit_aer": "qiskit-aer", "qiskit_ibm_runtime": "qiskit-ibm-runtime", "qiskit_mitigation": "qiskit-mitigation", "scipy": "scipy"}
_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 Parameter, QuantumCircuit
from qiskit.quantum_info import SparsePauliOp
warnings.filterwarnings("ignore")
def trotter_circuit(num_qubits, theta, layers, num_steps):
"""Trotterized transverse-field Ising evolution: RX on every site, then a CZ layer per coupling layer."""
circuit = QuantumCircuit(num_qubits)
for _ in range(num_steps):
circuit.rx(theta, range(num_qubits))
for layer in layers:
circuit.barrier()
for edge in layer:
circuit.cz(*edge)
circuit.sdg(edge)
circuit.barrier()
return circuit
num_qubits = 8
trotter_steps = 15
ring_coupling = [(i, (i + 1) % num_qubits) for i in range(num_qubits)]
layers = [
ring_coupling[0::2],
ring_coupling[1::2],
] # two disjoint layers of CZ gates
theta = Parameter("theta")
circuit = trotter_circuit(num_qubits, theta, layers, trotter_steps)
# Parameter values: theta in [0, 0.25 pi] - values for the transverse-field sweep
data_points = 6
parameter_values = np.linspace(0, 0.25 * np.pi, data_points).reshape(
(data_points, 1)
)
# Define observables - Z on each site
observables = [
SparsePauliOp("I" * (num_qubits - q - 1) + "Z" + "I" * q)
for q in range(num_qubits)
]
print("Visualization of 1 Trotter step")
trotter_circuit(num_qubits, theta, layers, 1).draw("mpl", scale=0.7)
Visualization of 1 Trotter step
Výpočet přesných středních hodnot pomocí simulátoru stavového vektoru v Qiskitu
Použij třídu Statevector z Qiskitu k výpočtu přesné střední hodnoty na každém uzlu pro každou hodnotu .
import matplotlib.pyplot as plt
from qiskit.quantum_info import Statevector
exact_evs = np.array(
[
[
Statevector(circuit.assign_parameters([value]))
.expectation_value(obs)
.real
for value in parameter_values.ravel()
]
for obs in observables
]
) # shape (observables, parameter values)
x_vals = parameter_values.ravel() / np.pi
fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(x_vals, exact_evs.mean(axis=0), "ko-", label="Exact")
ax.set_xlabel(r"$\theta/\pi$")
ax.set_ylabel(r"$\langle Z \rangle_{\text{avg}}$")
ax.legend()
plt.tight_layout()
plt.show()
Výpočet zašuměných, nezmírněných středních hodnot
QuantumProgram obsahuje kvantový obvod a samplex ve formátu vhodném pro primitivu Executor. MitigationTask.prepare sestaví QuantumProgram z obvodu, pozorovatelných a hodnot parametrů a obvod ohraničí do boxů podle zadaných voleb boxingu. Nastavíš volbu boxingu measure_annotations na "change_basis", která otočí měření do požadovaných bází, a přidáš add_tags: "unique_box", aby simulátor mohl namapovat model šumu na boxy obvodu. Program je ve výchozím nastavení Pauliho twirlován, což na hardwaru mění koherentní šum na Pauliho šum, který zmírňovací metody předpokládají. num_randomizations určuje, kolik twirlovaných instancí každého obvodu se vzorkuje, a shots_per_randomization určuje, kolik shotů každá instance dostane, takže jejich součin je rozpočet shotů na obvod.
Uživatel s přístupem k QPU by program spustil pomocí Executor sestaveného z backendu QiskitRuntimeService. Zde předáš Executor simulátor AerSimulator, který program spustí lokálně, a pomocí jeho simulátorové volby layer_noise_model vložíš náhodný Pauli-Lindbladův model šumu na každou vrstvu provázávacích hradel, vygenerovaný pomocí random_pauli_lindblad_noise z qiskit-mitigation s chybou na vrstvené hradlo (EPLG) 0.003, typickou pro současné QPU Heron.
from qiskit.circuit import CircuitInstruction
from qiskit.quantum_info import PauliLindbladMap, pauli_basis
from qiskit_aer import AerSimulator
from qiskit_ibm_runtime import Executor
from qiskit_mitigation import MitigationTask, find_combined_unique_layers
def generate_random_pauli_lindblad_noise(
layers,
num_qubits: int,
eplg: float = 3e-3,
seed: int | None = None,
) -> list[tuple[CircuitInstruction, 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 contain_measurement: # do not add measurement noise
layer_noise = PauliLindbladMap.identity(num_qubits)
else:
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))
return layer_noise_model
num_randomizations = 50
shots_per_randomization = 64
boxing_options = {
"measure_annotations": "change_basis",
"add_tags": "unique_box",
}
# find the unique layers to generate the layered noise for
box_layers = find_combined_unique_layers(
[circuit], custom_boxing_options=boxing_options
)
layer_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(),
options={"simulator": {"layer_noise_model": layer_noise_model}},
)
task = MitigationTask()
program = task.prepare(
circuit=circuit,
observables=observables,
parameters=parameter_values,
num_randomizations=num_randomizations,
shots_per_randomization=shots_per_randomization,
custom_boxing_options=boxing_options,
)
program._semantic_role = "simulation" # flag the program for the simulator
unmitigated_res = task.postprocess(noisy_executor.run(program).result())
fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(x_vals, exact_evs.mean(axis=0), "k-", label="Exact")
ax.errorbar(
x_vals,
np.mean(unmitigated_res.data.evs, axis=0),
np.mean(unmitigated_res.data.stds, axis=0),
fmt="o-",
label="Unmitigated",
)
ax.set_xlabel(r"$\theta/\pi$")
ax.set_ylabel(r"$\langle Z \rangle_{\text{avg}}$")
ax.legend()
plt.tight_layout()
plt.show()
Příprava programu ZNE s digitálním hradlovým foldingem a vzorkování obvodů s různými úrovněmi šumu
Připrav program ZNE obsahující obvody s rostoucími úrovněmi šumu. Každý faktor šumu představuje násobek provázávacích hradel, které obvod obsahuje (faktor šumu 3.0 obsahuje třikrát více provázávacích hradel než původní obvod). Obvod každého faktoru šumu je twirlován se stejným počtem randomizací a shotů jako nezmírněný program.
from qiskit_mitigation import GateFolding
task_gf = GateFolding()
program_zne = task_gf.prepare(
circuit=circuit,
observables=observables,
parameters=parameter_values,
num_randomizations=num_randomizations,
shots_per_randomization=shots_per_randomization,
custom_boxing_options=boxing_options,
noise_factors=[1, 3, 5], # ZNE option
)
Vykresli první Trotterův krok twirlované šablony s faktorem šumu 1 a s faktorem šumu 3, abys zobrazil zesílení. Pro faktor šumu 3.0 je každá vrstva obvodu nahrazena vrstvou, její inverzí a znovu vrstvou.
import IPython.display
def first_trotter_step(template, noise_factor):
"""The template up to the closing barrier of its first Trotter step (two layers, each folded ``noise_factor`` times)."""
toy = template.copy_empty_like()
for instruction in template.data:
toy.append(instruction)
if (
instruction.operation.name == "barrier"
and instruction.operation.label.startswith(
f"R{2 * noise_factor - 1}@"
)
):
break
return toy
IPython.display.display(
first_trotter_step(program_zne.items[0].circuit, 1).draw(
"mpl", scale=0.7, plot_barriers=False
)
)
first_trotter_step(program_zne.items[1].circuit, 3).draw(
"mpl", fold=-1, plot_barriers=False
)
Nakonec spusť program na zašuměném simulátoru. Jeho výsledky vstupují do extrapolace v dalším kroku.
program_zne._semantic_role = (
"simulation" # flag the program for the simulator
)
job_zne_results = noisy_executor.run(program_zne).result()
Extrapolace k limitě nulového šumu pomocí vestavěných a vlastních fitovacích funkcí
Nyní, když máš data pro každou úroveň šumu, můžeš do dat nafitovat funkci a extrapolovat zpět k hodnotě při nulovém šumu, abys odhadl střední hodnotu bezšumové pozorovatelné. Existuje mnoho způsobů, jak data fitovat, a správná volba funkce závisí na povaze vzorkovaných dat. GateFolding.postprocess nabízí několik vestavěných extrapolačních funkcí, včetně linear, exponential a double_exponential. Metoda postprocess přijímá seznam extrapolačních metod, vyhodnocuje je v pořadí a vrací první hodnotu s fyzikálně platnými hodnotami a konečným rozptylem. To znamená, že pořadí, v jakém extrapolátory zadáš, má velký vliv na to, který extrapolátor je vybrán. Zadání fallback jako extrapolátoru znamená, že pokud nebyly nalezeny žádné platné odpovědi, bude vybrána a vrácena střední hodnota nejmenšího faktoru šumu. Pokud fallback zadán není a žádný extrapolátor nedá fyzikálně platnou hodnotu, vrátí se výsledek s nejmenším std. Při Pauliho šumu pozorovatelná klesá exponenciálně s faktorem šumu, takže s 30 vrstvami (15 Trotterových kroků, 2 vrstvy na krok) provázávacích hradel je přirozenou volbou exponenciální extrapolátor; můžeš také extrapolovat lineárně, což je běžná volba, když je pokles malý.
Pro použití vlastní fitovací funkce předej extrapolátor "custom" spolu s n-ticí custom_fit obsahující fitovací funkci, seznam s počátečním odhadem každého z jejích parametrů a meze. Ty se předají do scipy.optimize.curve_fit. Postprocessing lze opakovat s různými extrapolátory bez opětovného spuštění programu. Zde použiješ vlastní fitovací funkci druhé odmocniny. Pro vizualizaci extrapolace vykresli magnetizaci zprůměrovanou přes uzly vypočtenou pro každý faktor šumu pro , tři fity a průměr přes uzly extrapolovaných středních hodnot v bodě nulového šumu. Tři datové body vypadají téměř lineárně, a přesto lineární extrapolace skončí viditelně pod přesnou hodnotou. Jak se předpovídalo, exponenciální fit odhadne bezšumovou střední hodnotu nejúčinněji.
from scipy.optimize import curve_fit
zne_mitigated_exp = task_gf.postprocess(
job_zne_results, extrapolator=["exponential", "linear"]
)
zne_mitigated_lin = task_gf.postprocess(
job_zne_results, extrapolator=["linear"]
)
def sqrt_fit(x, a, b):
return a * np.sqrt(x) + b
custom_fit = (sqrt_fit, [1.0, 1.0], (-np.inf, np.inf))
zne_mitigated_custom = task_gf.postprocess(
job_zne_results,
extrapolator=["custom", "exponential", "linear"],
custom_fit=custom_fit,
)
# The result holds the expectation value at each noise factor and the extrapolated value, but
# not the fit parameters, so we refit the same models here to draw the curves.
noise_factors = np.array([1, 3, 5])
angle = 3 # theta = 0.15 pi, the angle with the largest unmitigated error
amplifications = zne_mitigated_exp.data.evs_noise_factors[:, angle, :].mean(
axis=0
)
x_fit = np.linspace(0, 5, 51)
fits = {
"exponential": (
lambda x, a, b: a * np.exp(b * x),
zne_mitigated_exp,
"black",
),
"linear": (lambda x, a, b: a * x + b, zne_mitigated_lin, "blue"),
"square-root": (sqrt_fit, zne_mitigated_custom, "green"),
}
plt.plot(
noise_factors, amplifications, "o", color="red", label="amplifications"
)
for name, (fit, res, color) in fits.items():
params, _ = curve_fit(fit, noise_factors, amplifications, p0=[1.0, -0.1])
plt.plot(
x_fit, fit(x_fit, *params), "--", color=color, label=f"{name} fit"
)
plt.plot(
[0],
res.data.evs[:, angle].mean(),
"*",
color=color,
markersize=15,
label=f"{name} extrapolation",
)
plt.plot(
[0],
exact_evs[:, angle].mean(),
"_",
color="red",
markersize=15,
markeredgewidth=3,
label="exact",
)
plt.xlabel("noise factor")
plt.ylabel(r"$\langle Z \rangle_{\text{avg}}$")
plt.legend()
plt.show()
Porovnání zmírněných středních hodnot
Vykresli magnetizaci zprůměrovanou přes uzly pro každou metodu napříč celým rozsahem Isingových rotačních úhlů. Exponenciální extrapolace obnoví přesnou magnetizaci v rámci svých chybových úseček, zatímco lineární extrapolace odstraní velkou část zkreslení, ale soustavně podstřeluje, protože skutečný pokles je exponenciální. Fit druhou odmocninou silně přestřeluje.
results = {
"Unmitigated": unmitigated_res,
"ZNE exponential": zne_mitigated_exp,
"ZNE linear": zne_mitigated_lin,
"ZNE square-root": zne_mitigated_custom,
}
fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(x_vals, exact_evs.mean(axis=0), "k-", label="Exact")
for label, res in results.items():
ax.errorbar(
x_vals,
np.mean(res.data.evs, axis=0),
np.mean(res.data.stds, axis=0),
fmt="o-",
label=label,
)
ax.set_xlabel(r"$\theta/\pi$")
ax.set_ylabel(r"$\langle Z \rangle_{\text{avg}}$")
ax.legend()
ax.set_title(
f"Error-Mitigated Magnetization [{num_qubits}-qubit, {trotter_steps} Trotter steps]"
)
plt.tight_layout()
plt.show()
