Ekstrapolacja do zerowego szumu z zagęszczaniem bramek i niestandardową ekstrapolacją
Ekstrapolacja do zerowego szumu (ZNE) to technika łagodzenia szumu kwantowego w obliczeniach wartości oczekiwanych obserwabli. Polega na wykonaniu obwodu kwantowego przy rosnących skalach szumu, obliczeniu wartości oczekiwanych względem tych skal szumu i ekstrapolacji w celu oszacowania wartości oczekiwanej przy zerowym szumie. Zarówno etap wzmacniania szumu, jak i etap ekstrapolacji można zaimplementować na różne sposoby. Ten przewodnik wzmacnia szum za pomocą „cyfrowego zagęszczania bramek” (digital gate folding), co oznacza, że bramki dwukubitowe są zastępowane równoważnymi ciągami składającymi się z bramki i jej odwrotności. Na przykład zastąpienie unitarnej przez daje współczynnik wzmocnienia szumu równy 3. Wzmacnianie można też zaimplementować za pomocą probabilistycznego wzmacniania błędów (PEA), które wymaga dodatkowego kroku uczenia modelu szumu i jest omówione w przewodniku o probabilistycznym wzmacnianiu błędów. Ten przewodnik ekstrapoluje przy użyciu kilku wbudowanych strategii dopasowania (liniowej i wykładniczej), a także pokazuje, jak zaimplementować własną funkcję dopasowania.
Pamiętaj, że choć ZNE może skutecznie łagodzić szum kwantowy przy niższym koszcie próbkowania niż inne metody, nie daje żadnych ograniczeń na obciążenie (bias).
Przepływ pracy wygląda następująco:
- Określ obwód i obserwablę
- Oblicz dokładne wartości oczekiwane za pomocą symulatora wektora stanu Qiskit
- Oblicz zaszumione, niełagodzone wartości oczekiwane
- Przygotuj program ZNE z cyfrowym zagęszczaniem bramek i próbkuj obwody o różnych skalach szumu
- Wyekstrapoluj do granicy zerowego szumu za pomocą wbudowanych i niestandardowych funkcji dopasowania
- Porównaj złagodzone wartości oczekiwane
Określ obwód i obserwablę
Oblicz średnią po węzłach magnetyzację trotteryzowanego modelu Isinga z polem poprzecznym (TFIM) na ośmiu kubitach w topologii pierścienia. Przeprowadzisz przemiatanie po kilku wartościach natężenia pola poprzecznego i narysujesz magnetyzację jako funkcję natężenia pola.
# 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
Oblicz dokładne wartości oczekiwane za pomocą symulatora wektora stanu Qiskit
Użyj klasy Statevector z Qiskit, aby obliczyć dokładną wartość oczekiwaną w każdym węźle dla każdej wartości .
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()
Oblicz zaszumione, niełagodzone wartości oczekiwane
QuantumProgram zawiera obwód kwantowy i samplex w formacie odpowiednim dla prymitywu Executor. MitigationTask.prepare buduje QuantumProgram z obwodu, obserwabli i wartości parametrów oraz otacza obwód pudełkami zgodnie z podanymi opcjami boxingu. Ustawiasz opcję boxingu measure_annotations na "change_basis", która obraca pomiary do wymaganych baz, i dodajesz add_tags: "unique_box", aby symulator mógł odwzorować model szumu na pudełka obwodu. Program jest domyślnie poddawany Pauli twirlingowi, co na sprzęcie zamienia szum koherentny w szum Pauliego, którego zakładają metody łagodzenia. num_randomizations określa, ile zatwirlowanych instancji każdego obwodu jest próbkowanych, a shots_per_randomization określa, ile strzałów otrzymuje każda instancja, więc ich iloczyn to budżet strzałów na obwód.
Użytkownik z dostępem do QPU wykonałby program za pomocą Executor zbudowanego z backendu QiskitRuntimeService. Tutaj przekażesz do Executor symulator AerSimulator, który uruchamia program lokalnie, i użyjesz jego opcji symulatora layer_noise_model, aby wstrzyknąć losowy model szumu Pauli-Lindblada na każdą warstwę bramek splątujących, wygenerowany za pomocą random_pauli_lindblad_noise z qiskit-mitigation z błędem na warstwową bramkę (EPLG) równym 0.003, typowym dla obecnych 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()
Przygotuj program ZNE z cyfrowym zagęszczaniem bramek i próbkuj obwody o różnych skalach szumu
Przygotuj program ZNE zawierający obwody o rosnących skalach szumu. Każdy współczynnik szumu oznacza krotność bramek splątujących w obwodzie (współczynnik szumu 3.0 oznacza trzy razy więcej bramek splątujących niż w oryginalnym obwodzie). Obwód dla każdego współczynnika szumu jest twirlowany z taką samą liczbą randomizacji i strzałów jak program niełagodzony.
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
)
Narysuj pierwszy krok Trottera zatwirlowanego szablonu dla współczynnika szumu 1 i dla współczynnika szumu 3, aby zobaczyć wzmocnienie. Dla współczynnika szumu 3.0 każda warstwa obwodu jest zastępowana przez warstwę, jej odwrotność i ponownie warstwę.
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
)
Na koniec uruchom program na zaszumionym symulatorze. Jego wyniki zasilają ekstrapolację w następnym kroku.
program_zne._semantic_role = (
"simulation" # flag the program for the simulator
)
job_zne_results = noisy_executor.run(program_zne).result()
Wyekstrapoluj do granicy zerowego szumu za pomocą wbudowanych i niestandardowych funkcji dopasowania
Skoro masz dane dla każdej skali szumu, możesz dopasować funkcję do danych i wyekstrapolować z powrotem do wartości przy zerowym szumie, aby oszacować wartość oczekiwaną obserwabli bez szumu. Istnieje wiele sposobów dopasowywania danych, a właściwy wybór funkcji zależy od natury próbkowanych danych. GateFolding.postprocess udostępnia kilka wbudowanych funkcji ekstrapolacji, w tym linear, exponential i double_exponential. Metoda postprocess przyjmuje listę metod ekstrapolacji, ocenia je po kolei i zwraca pierwszą wartość o fizycznie poprawnych wartościach i skończonej wariancji. Oznacza to, że kolejność, w jakiej podasz ekstrapolatory, ma duży wpływ na to, który z nich zostanie wybrany. Podanie fallback jako ekstrapolatora oznacza, że jeśli nie znaleziono żadnych poprawnych odpowiedzi, zostanie wybrana i zwrócona wartość oczekiwana dla najmniejszego współczynnika szumu. Jeśli fallback nie zostanie podany i żaden ekstrapolator nie da fizycznie poprawnej wartości, zwracany jest wynik o najmniejszym std. Przy szumie Pauliego obserwabla zanika wykładniczo ze współczynnikiem szumu, więc przy 30 warstwach (15 kroków Trottera, 2 warstwy na krok) bramek splątujących naturalnym wyborem jest ekstrapolator wykładniczy; możesz też ekstrapolować liniowo, co jest częstym wyborem, gdy zanik jest niewielki.
Aby użyć niestandardowej funkcji dopasowania, przekaż ekstrapolator "custom" razem z krotką custom_fit zawierającą funkcję dopasowania, listę z początkowym oszacowaniem każdego z jej parametrów oraz ograniczenia. Są one przekazywane do scipy.optimize.curve_fit. Przetwarzanie końcowe można powtarzać z różnymi ekstrapolatorami bez ponownego wykonywania programu. Tutaj używasz niestandardowej funkcji dopasowania typu pierwiastek kwadratowy. Aby zwizualizować ekstrapolację, narysuj średnią po węzłach magnetyzację obliczoną dla każdego współczynnika szumu dla , trzy dopasowania oraz średnią po węzłach wyekstrapolowanych wartości oczekiwanych w punkcie zerowego szumu. Trzy punkty danych wyglądają niemal liniowo, a mimo to ekstrapolacja liniowa ląduje wyraźnie poniżej dokładnej wartości. Zgodnie z przewidywaniami dopasowanie wykładnicze najskuteczniej szacuje wartość oczekiwaną bez szumu.
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()
Porównaj złagodzone wartości oczekiwane
Narysuj średnią po węzłach magnetyzację dla każdej metody w całym przemiataniu po kątach obrotu Isinga. Ekstrapolacja wykładnicza odtwarza dokładną magnetyzację w granicach swoich słupków błędów, podczas gdy ekstrapolacja liniowa usuwa dużą część obciążenia, ale konsekwentnie ją zaniża, ponieważ prawdziwy zanik jest wykładniczy. Dopasowanie pierwiastkowe znacznie ją zawyża.
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()
