jupytext | kernelspec | ||||||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
|
|
Corresponding to: Min Li (minl2@illinois.edu)
This notebook is a prototype of how to perform robust shadow estimation protocol with mitiq.
import cirq
import numpy as np
from typing import List
from mitiq.shadows.shadows import *
from mitiq.shadows.quantum_processing import *
from mitiq.shadows.classical_postprocessing import *
from mitiq.shadows.shadows_utils import *
from mitiq import MeasurementResult
from mitiq.interface.mitiq_cirq.cirq_utils import (
sample_bitstrings as cirq_sample_bitstrings,
)
# set random seed
np.random.seed(666)
There are two options: Whether to run the quantum measurement or directly use the results from the previous run.
- If True, the measurement will be run again.
- If False, the results from the previous run will be used.
import zipfile, pickle, io, requests, os
run_quantum_processing = False
run_pauli_twirling_calibration = False
file_directory = "./resources"
if not run_quantum_processing:
saved_data_name = "saved_data-rshadows"
with zipfile.ZipFile(
f"{file_directory}/rshadows-tutorial-{saved_data_name}.zip"
) as zf:
saved_data = pickle.load(zf.open(f"{saved_data_name}.pkl"))
The robust shadow estimation{cite}chen2021robust
approach based on {cite}huang2020predicting
exhibits noise resilience. The inherent randomization in the protocol simplifies the noise, transforming it into a Pauli noise channel that can be characterized relatively straightforwardly. Once the noisy channel
In this notebook, we'll use the ground state of the Ising model with periodic boundary conditions as an example to study energy and two-point correlation function estimations. We'll compare the performance of robust shadow estimation with the standard shadow protocol, taking into account the bit-flip or depolarization noise on the quantum channels.
The Hamiltonian of the Ising model is given by
\begin{equation}
H = -J\sum_{i=0}^{n-1} Z_i Z_{i+1} - g\sum_{i=1}^N X_i,
\end{equation}
We focus on the case where
# import groud state of 1-D Ising model w/ periodic boundary condition
download_ising_circuits = True
num_qubits = 8
qubits: List[cirq.Qid] = cirq.LineQubit.range(num_qubits)
if download_ising_circuits:
with open(f"{file_directory}/rshadows-tutorial-1D_Ising_g=1_{num_qubits}qubits.pkl", "rb") as file:
old_cirq_circuit = pickle.load(file)
circuit = cirq.Circuit(old_cirq_circuit.all_operations())
g = 1
# or user can import from tensorflow_quantum
else:
from tensorflow_quantum.datasets import tfi_chain
qbs = cirq.GridQubit.rect(num_qubits, 1)
circuits, labels, pauli_sums, addinfo = tfi_chain(qbs, "closed")
lattice_idx = 40 # Critical point where g == 1
g = addinfo[lattice_idx].g
circuit = circuits[lattice_idx]
qubit_map = {
cirq.GridQubit(i, 0): cirq.LineQubit(i) for i in range(num_qubits)
}
circuit = circuit.transform_qubits(qubit_map=qubit_map)
Similar with the classical shadow protocol, we define the executor to perform the computational measurement for the circuit. Here, we use add single-qubit depolarizing noise after rotation gates but before the
def cirq_executor(
circuit: cirq.Circuit,
noise_model_function=cirq.depolarize,
noise_level=(0.2,),
sampler=cirq.Simulator(),
) -> MeasurementResult:
"""
This function returns the measurement outcomes of a circuit with noisy
channel added right before measurement.
Args:
circuit: The circuit to execute.
Returns:
A one shot MeasurementResult object containing the measurement
outcomes.
"""
tmp_circuit = circuit.copy()
qubits = sorted(list(tmp_circuit.all_qubits()))
if noise_level[0] > 0:
noisy_circuit = cirq.Circuit()
operations = list(tmp_circuit)
n_ops = len(operations)
for i, op in enumerate(operations):
if i == n_ops - 1:
noisy_circuit.append(
cirq.Moment(
*noise_model_function(*noise_level).on_each(*qubits)
)
)
noisy_circuit.append(op)
tmp_circuit = noisy_circuit
# circuit.append(cirq.Moment(*noise_model_function(*noise_level).on_each(*qubits)))
executor = cirq_sample_bitstrings(
tmp_circuit,
noise_model_function=None,
noise_level=(0,),
shots=1,
sampler=sampler,
)
return executor
The PTM (Pauli Transfer Matrix) or Liouville representation provides a vector representation for all linear operators
from mitiq.utils import operator_ptm_vector_rep
operator_ptm_vector_rep(cirq.I._unitary_() / np.sqrt(2))
The classical shadow estimation involves Pauli twirling of a quantum channel represented by get_single_shot_pauli_fidelity
and mitiq.shadows.classical_postprocessing.get_pauli_fidelities
for implementation.
In the ideal noise-free scenario, Pauli fidelity is:
\begin{equation}
\hat{f}_{b}^{\mathrm{ideal}} = 3^{-|{b}|}
\end{equation}
For noisy channels, the inverse channel
from functools import partial
n_total_measurements_calibration = 20000
if run_quantum_processing:
noisy_executor = partial(cirq_executor, noise_level=(0.1,))
zero_state_shadow_output = shadow_quantum_processing(
# zero circuit of 8 qubits
circuit=cirq.Circuit(),
num_total_measurements_shadow=n_total_measurements_calibration,
executor=noisy_executor,
qubits=qubits,
)
else:
zero_state_shadow_output = saved_data["shadow_outcomes_f_plot"]
f_est_results = pauli_twirling_calibrate(
zero_state_shadow_outcomes=zero_state_shadow_output,
k_calibration=5,
locality=2,
)
# sort bitstrings (b_lists' string rep) by number of 1s
bitstrings = np.array(sorted(list(f_est_results.keys())))
# reorder f_est_results by number of '1' in bitstrings
counts = {bitstring: bitstring.count("1") for bitstring in bitstrings}
order = np.argsort(list(counts.values()))
reordered_bitstrings = bitstrings[order]
# solve for theoretical Pauli fidelities
f_theoretical = {}
bitstrings = list(f_est_results.keys())
for bitstring in bitstrings:
n_ones = bitstring.count("1")
f_val = 3 ** (-n_ones)
f_theoretical[bitstring] = f_val
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style("whitegrid")
# plot estimated vs theoretical Pauli fidelities when no errors in quantum circuit
plt.plot(
[np.abs(f_est_results[b]) for b in reordered_bitstrings],
"-*",
label="Noisy Channel",
)
plt.plot(
[f_theoretical[b] for b in reordered_bitstrings], label="Noiseless Channel"
)
plt.xlabel(r"measurement basis states $b$")
plt.xticks(
range(len(reordered_bitstrings)),
reordered_bitstrings,
rotation="vertical",
fontsize=6,
)
plt.ylabel("Pauli fidelity")
plt.legend();
The expectation value for a series of operators, denoted as ${O_\iota}{\iota\leq M}$, has a snapshot version of expectation value estimation by random Pauli measurement $\widetilde{\mathcal{M}}=\bigotimes{i}\widetilde{\mathcal{M}}{P_i}$ and Pauli-twirling calibration $\widehat{\mathcal{M}}^{-1}=\sum{b\in{0,1}^n}f_b^{-1}\bigotimes_{i}\Pi_{b_i\in b}$, which is given by
\begin{align}
\hat{o}\iota &= \langle!\langle O\iota|{\hat{\rho}}\rangle!\rangle \simeq \langle!\langle O_\iota|\widehat{\mathcal{M}}^{-1}\widetilde{\mathcal{M}}|\rho\rangle!\rangle=\sum_{b^{(1)}\in{0,1}^{n}}f_{b^{(1)}}^{-1}\left(\bigotimes_{i=1}^n \langle!\langle P_i|\Pi_{b_i^{(1)}}\widehat{\mathcal{M}}{P_i}\right)|\rho\rangle!\rangle\nonumber\
&=\sum{b^{(1)}\in{0,1}^{n}}f_{b^{(1)}}^{-1}\prod_{i=1}^n \langle!\langle P_i|\Pi_{b^{(1)}i}\bigg|U_i^{(2)\dagger}|b_i^{(2)}\rangle\langle b_i^{(2)}|U_i^{(2)}\bigg\rangle!\bigg\rangle
\end{align}
where in the last equality, ${P_i}{i\in n}$ represents Pauli operators, with
Next steps are identical to the classical protocol, where the statistical method of taking an average called "median of means" is used to achieve an acceptable failure probability of estimation. This needs \mathrm{where}\hat{o}\iota^{(j)}=N_2^{-1}\sum_{k=N_2(j-1)+1}^{N_2j}\mathrm{Tr}(O_\iota\hat{\rho}k),\qquad \forall~1\leq j\leq K_2,
\end{align}
where we have $K_2$ estimators each of which is the average of $N_2$ single-round estimators $\hat{o}i^{(j)}$, and take the median of these $K_2$ estimators as our final estimator $\hat{o}\iota(N_2,K_2)$. We can calculate the median of means of each irreducible representations with projection $\Pi_b=\bigotimes{i=1}^n\Pi_{b_i}$,
In this section, we will use the robust shadows estimation algorithm to estimate the ground state energy of the Ising model. We will use the compare_shadow_methods
function to compare the performance of the robust shadows estimation algorithm and the classical shadows estimation algorithm. The compare_shadow_methods
function takes the following parameters:
def compare_shadow_methods(
circuit,
observables,
n_measurements_calibration,
k_calibration,
n_measurement_shadow,
k_shadows,
locality,
noisy_executor,
run_quantum_processing,
shadow_measurement_result=None,
zero_state_shadow_output=None,
):
if run_quantum_processing:
zero_state_shadow_output = shadow_quantum_processing(
circuit=cirq.Circuit(),
num_total_measurements_shadow=n_measurements_calibration,
executor=noisy_executor,
qubits=qubits,
)
shadow_measurement_result = shadow_quantum_processing(
circuit,
num_total_measurements_shadow=n_measurement_shadow,
executor=noisy_executor,
)
else:
assert shadow_measurement_result is not None
assert zero_state_shadow_output is not None
file_zsso = zero_state_shadow_output[1][0]
file_k_cal = k_calibration
file_locality = locality
file_name = f"rshadows-tutorial-{file_zsso}-{file_k_cal}-{file_locality}"
if not run_pauli_twirling_calibration and os.path.exists(f"{file_directory}/{file_name}.pkl"):
# use the file
with open(f"{file_directory}/{file_name}.pkl", "rb") as file:
f_est = pickle.load(file)
else:
f_est = pauli_twirling_calibrate(
zero_state_shadow_outcomes=zero_state_shadow_output,
k_calibration=k_calibration,
locality=locality,
)
output_shadow = classical_post_processing(
shadow_outcomes=shadow_measurement_result,
observables=observables,
k_shadows=k_shadows,
)
output_shadow_cal = classical_post_processing(
shadow_outcomes=shadow_measurement_result,
calibration_results=f_est,
observables=observables,
k_shadows=k_shadows,
)
return {"standard": output_shadow, "robust": output_shadow_cal}
We use the groud state of 1-D Ising model with periodic boundary condition, with
# define obersevables lists as Ising model Hamiltonian
from mitiq import PauliString
ising_hamiltonian = [
PauliString("X", support=(i,), coeff=-g) for i in range(num_qubits)
] + [
PauliString("ZZ", support=(i, (i + 1) % num_qubits), coeff=-1)
for i in range(num_qubits)
]
Calculate the exact expectation values of the Pauli operators for the above state:
state_vector = circuit.final_state_vector()
expval_exact = []
for i, pauli_string in enumerate(ising_hamiltonian):
exp = pauli_string._pauli.expectation_from_state_vector(
state_vector, qubit_map={q: i for i, q in enumerate(qubits)}
)
expval_exact.append(exp.real)
We use bit_flip channel as an example to show how to use the robust shadow estimation (RSE) in Mitiq. The bit_flip channel is a common noise model in quantum computing. It is a Pauli channel that flips the state of a qubit with probability
noise_levels = np.linspace(0, 0.06, 4)
# if noise_model is None, then the noise model is depolarizing noise
noise_model = "bit_flip"
standard_results = []
robust_results = []
noise_model_fn = getattr(cirq, noise_model)
for noise_level in noise_levels:
noisy_executor = partial(
cirq_executor,
noise_level=(noise_level,),
noise_model_function=cirq.bit_flip,
)
experiment_name = f"{num_qubits}qubits_{noise_model}_{noise_level}"
if run_quantum_processing:
shadow_measurement_result, zero_state_shadow_output = None, None
else:
shadow_measurement_result = saved_data[experiment_name][
"shadow_outcomes"
]
zero_state_shadow_output = saved_data[experiment_name][
"zero_shadow_outcomes"
]
est_values = compare_shadow_methods(
circuit=circuit,
observables=ising_hamiltonian,
n_measurements_calibration=60000,
n_measurement_shadow=60000,
k_shadows=6,
locality=3,
noisy_executor=noisy_executor,
k_calibration=10,
run_quantum_processing=False,
shadow_measurement_result=shadow_measurement_result,
zero_state_shadow_output=zero_state_shadow_output,
)
standard_results.append(est_values["standard"])
robust_results.append(est_values["robust"])
import pandas as pd
df_energy = pd.DataFrame(
columns=["noise_level", "method", "observable", "value"]
)
for i, noise_level in enumerate(noise_levels):
est_values = {}
est_values["standard"] = list(standard_results[i].values())
est_values["robust"] = list(robust_results[i].values())
# for j in range(len(standard_est_values)):
df_energy = pd.concat(
[
df_energy,
pd.DataFrame(
{
"noise_level": noise_level,
"method": "exact",
"observable": [str(ham) for ham in ising_hamiltonian],
"value": expval_exact,
}
),
],
ignore_index=True,
)
for method in ["standard", "robust"]:
df_energy = pd.concat(
[
df_energy,
pd.DataFrame(
{
"noise_level": noise_level,
"method": method,
"observable": [str(ham) for ham in ising_hamiltonian],
"value": est_values[method],
}
),
],
ignore_index=True,
)
df_hamiltonian = df_energy.groupby(["noise_level", "method"]).sum()
df_hamiltonian = df_hamiltonian.reset_index()
noise_model = "bit_flip"
# Define a color palette
palette = {"exact": "black", "robust": "red", "standard": "green"}
plt.figure()
sns.lineplot(
data=df_hamiltonian,
x="noise_level",
y="value",
hue="method",
palette=palette, # Use the color palette defined above
markers=True,
style="method",
dashes=False,
errorbar=("ci", 95),
)
plt.title(f"Hamiltonian Estimation for {noise_model} Noise")
plt.xlabel("Noise Level")
plt.ylabel("Energy Value");
Let's estimate two point correlation fuction:
Import groud state of 1-D Ising model with periodic boundary condition
num_qubits = 16
qubits = cirq.LineQubit.range(num_qubits)
if download_ising_circuits:
with open(f"{file_directory}/rshadows-tutorial-1D_Ising_g=1_{num_qubits}qubits.pkl", "rb") as file:
old_cirq_circuit = pickle.load(file)
circuit = cirq.Circuit(old_cirq_circuit.all_operations())
g = 1
else:
qbs = cirq.GridQubit.rect(num_qubits, 1)
circuits, labels, pauli_sums, addinfo = tfi_chain(qbs, "closed")
lattice_idx = 40 # Critical point where g == 1
g = addinfo[lattice_idx].g
circuit = circuits[lattice_idx]
qubit_map = {
cirq.GridQubit(i, 0): cirq.LineQubit(i) for i in range(num_qubits)
}
circuit = circuit.transform_qubits(qubit_map=qubit_map)
Define obersevables lists as two point correlation functions between the first qubit and the rest of the qubits
two_pt_correlation = [
PauliString("ZZ", support=(0, i), coeff=-1) for i in range(1, num_qubits, 2)
]
Calculate the exact correlation function
expval_exact = []
state_vector = circuit.final_state_vector()
for i, pauli_string in enumerate(two_pt_correlation):
exp = pauli_string._pauli.expectation_from_state_vector(
state_vector, qubit_map={q: i for i, q in enumerate(qubits)}
)
expval_exact.append(exp.real)
with depolarizing noise set to
noisy_executor = partial(cirq_executor, noise_level=(0.1,))
experiment_name = f"{num_qubits}qubits_depolarize_{noise_level}"
shadow_measurement_result = saved_data[experiment_name]["shadow_outcomes"]
zero_state_shadow_output = saved_data[experiment_name]["zero_shadow_outcomes"]
est_values = compare_shadow_methods(
circuit=circuit,
observables=two_pt_correlation,
n_measurements_calibration=50000,
n_measurement_shadow=50000,
k_shadows=5,
locality=2,
noisy_executor=noisy_executor,
k_calibration=5,
run_quantum_processing=False,
shadow_measurement_result=shadow_measurement_result,
zero_state_shadow_output=zero_state_shadow_output,
)
df_corr = pd.DataFrame(
columns=["method", "qubit_index", "observable", "value"]
)
qubit_idxes = [max(corr.support()) for corr in two_pt_correlation]
# for j in range(len(standard_est_values)):
for method in ["standard", "robust"]:
df_corr = pd.concat(
[
df_corr,
pd.DataFrame(
{
"method": method,
"qubit_index": qubit_idxes,
"observable": [str(corr) for corr in two_pt_correlation],
"value": list(est_values[method].values()),
}
),
],
ignore_index=True,
)
df_corr = pd.concat(
[
df_corr,
pd.DataFrame(
{
"method": "exact",
"qubit_index": qubit_idxes,
"observable": [str(corr) for corr in two_pt_correlation],
"value": expval_exact,
}
),
],
ignore_index=True,
)
# Define a color palette
palette = {"exact": "black", "robust": "red", "standard": "green"}
plt.figure()
sns.lineplot(
data=df_corr,
x="qubit_index",
y="value",
hue="method",
palette=palette, # Use the color palette defined above
markers=True,
style="method",
dashes=False,
errorbar=("ci", 95),
)
plt.title("Correlation Function Estimation w/ 0.3 Depolarization Noise")
plt.xlabel(r"Correlation Function $\langle Z_0Z_i \rangle$")
plt.ylabel("Correlation");