Designing noise-robust single-qubit gates for IBM Qiskit
Increasing robustness against dephasing and control noise using Boulder Opal pulses
Boulder Opal enables you to design numerically optimized controls which can improve the performance of quantum computing hardware. Specifically, robust control solutions are able to reduce sensitivity to dephasing and/or control noise, including slowly varying noise sources that arise as hardware drifts over time.
In this application note we present an entire workflow - from gate optimization through to experimental validation on IBM Quantum Hardware. We will cover:
- Designing optimized dephasing-noise-robust controls using graphs
- Incorporating band-limits on the optimization to produce smooth pulses
- Validating performance using filter functions and quasi-static-noise-susceptibility scans
- Formatting for Qiskit and experimental execution on IBM Quantum hardware.
A detailed discussion of optimized-pulse performance on IBM hardware can be found in our technical manuscript, “Error-robust quantum logic optimization using a cloud quantum computer interface”, published as Phys. Rev. Applied 15, 064054 (2021).
Some cells in this notebook require an account with IBM-Q to execute correctly. If you want to run them, please go to the IBM-Q experience to set up an account.
Imports and initialization
import time
import warnings
from pathlib import Path
# Disable warning from qiskit
warnings.filterwarnings("ignore")
# Choose to run experiments or to use saved data
use_saved_data = True
use_IBM = False
if use_saved_data == False:
timestr = time.strftime("%Y%m%d-%H%M%S")
print("Time label for data saved throughout this experiment:" + timestr)
data_folder = Path("resources/")
import jsonpickle
import matplotlib.pyplot as plt
# General imports
import numpy as np
from scipy import interpolate
from qctrlvisualizer import QCTRL_STYLE_COLORS, get_qctrl_style, plot_controls
# Q-CTRL imports
from qctrl import Qctrl
# Parameters
giga = 1e9
mega = 1e6
nano = 1e-9
plt.style.use(get_qctrl_style())
colors = {
"Q-CTRL": QCTRL_STYLE_COLORS[0],
"Square": QCTRL_STYLE_COLORS[1],
"IBM default X-gate": QCTRL_STYLE_COLORS[2],
}
bloch_colors = {
"x": QCTRL_STYLE_COLORS[0],
"y": QCTRL_STYLE_COLORS[1],
"z": QCTRL_STYLE_COLORS[2],
}
bloch_markers = {"x": "x", "y": "s", "z": "o"}
bloch_lines = {"x": "--", "y": "-.", "z": "-"}
bloch_basis = ["x", "y", "z"]
# Auxiliary functions
def get_detuning_scan(detunings, control):
graph = qctrl.create_graph()
shift_I = graph.pwc(*qctrl.utils.pwc_pairs_to_arrays(control["I"]))
shift_Q = graph.pwc(*qctrl.utils.pwc_pairs_to_arrays(control["Q"]))
duration = np.sum(shift_I.durations)
noise = graph.constant_pwc(
constant=detunings, duration=duration, batch_dimension_count=1
)
hamiltonian = (
shift_I * graph.pauli_matrix("X") / 2
+ shift_Q * graph.pauli_matrix("Y") / 2
+ noise * graph.pauli_matrix("Z")
)
target = np.array([[1, 0], [0, 0]])
unitaries = graph.time_evolution_operators_pwc(
hamiltonian=hamiltonian, sample_times=np.array([duration])
)
infidelities = graph.unitary_infidelity(
unitaries[:, 0], target, name="infidelities"
)
result = qctrl.functions.calculate_graph(
graph=graph, output_node_names=["infidelities"]
)
return 1 - result.output["infidelities"]["value"]
def get_amplitude_scan(amplitudes, control):
graph = qctrl.create_graph()
drive_I = graph.pwc(*qctrl.utils.pwc_pairs_to_arrays(control["I"]))
drive_Q = graph.pwc(*qctrl.utils.pwc_pairs_to_arrays(control["Q"]))
drive = drive_I + 1j * drive_Q
duration = np.sum(drive.durations)
noise = graph.constant_pwc(
constant=amplitudes, duration=duration, batch_dimension_count=1
)
target = np.array([[1, 0], [0, 0]])
hamiltonian = graph.hermitian_part((1 + noise) * drive * graph.pauli_matrix("P"))
unitaries = graph.time_evolution_operators_pwc(
hamiltonian=hamiltonian, sample_times=np.array([duration])
)
infidelities = graph.unitary_infidelity(
unitaries[:, 0], target, name="infidelities"
)
result = qctrl.functions.calculate_graph(
graph=graph, output_node_names=["infidelities"]
)
return 1 - result.output["infidelities"]["value"]
def simulation_coherent(control, time_samples):
graph = qctrl.create_graph()
drive_I = graph.pwc(*qctrl.utils.pwc_pairs_to_arrays(control["I"]))
drive_Q = graph.pwc(*qctrl.utils.pwc_pairs_to_arrays(control["Q"]))
drive = drive_I + 1j * drive_Q
duration = np.sum(drive.durations)
hamiltonian = graph.hermitian_part(drive * graph.pauli_matrix("P"))
initial_state_vector = np.array([1.0, 0.0])
sample_times = np.linspace(0, duration, time_samples)
unitaries = graph.time_evolution_operators_pwc(
hamiltonian=hamiltonian, sample_times=sample_times
)
states = unitaries @ initial_state_vector[:, None]
for name in ["X", "Y", "Z"]:
graph.real(
graph.expectation_value(states[..., 0], graph.pauli_matrix(name)), name=name
)
result = qctrl.functions.calculate_graph(
graph=graph, output_node_names=["X", "Y", "Z"]
)
return {k.lower(): v["value"] for k, v in result.output.items()}, sample_times
def save_var(file_name, var):
# Save a single variable to a file using jsonpickle
f = open(file_name, "w+")
to_write = jsonpickle.encode(var)
f.write(to_write)
f.close()
def load_var(file_name):
# Return a variable from a json file
f = open(file_name, "r+")
encoded = f.read()
decoded = jsonpickle.decode(encoded)
f.close()
return decoded
# Starting a session with the API
qctrl = Qctrl()
# IBM-Q imports
if use_IBM == True:
import qiskit.pulse as pulse
from qiskit import IBMQ
from qiskit.compiler import assemble
from qiskit.pulse import (
Acquire,
AcquireChannel,
DriveChannel,
MeasureChannel,
MemorySlot,
Play,
)
from qiskit.pulse.library import Waveform
from qiskit.tools.jupyter import *
from qiskit.tools.monitor import job_monitor
# IBM credentials and backend selection
IBMQ.enable_account("your IBM token")
provider = IBMQ.get_provider(
hub="your hub", group="your group", project="your project"
)
backend = provider.get_backend("ibmq_armonk")
backend_defaults = backend.defaults()
backend_config = backend.configuration()
# Backend properties
qubit = 0
qubit_freq_est = backend_defaults.qubit_freq_est[qubit]
dt = backend_config.dt
print(f"Qubit: {qubit}")
print(f"Hardware sampling time: {dt/nano} ns")
print(f"Qubit frequency estimate: {qubit_freq_est/giga} GHz")
# Drive and measurement channels
drive_chan = DriveChannel(qubit)
meas_chan = MeasureChannel(qubit)
inst_sched_map = backend_defaults.instruction_schedule_map
measure_schedule = inst_sched_map.get(
"measure", qubits=backend_config.meas_map[qubit]
)
else:
# Use default dt for armonk backend
dt = 2 / 9 * nano
# Use 0 qubit for armonk
qubit = 0
# Use last frequency update
qubit_freq_est = 4974444325.742604
qubit_frequency_updated = qubit_freq_est
print("IBM offline parameters")
print(f"Qubit: {qubit}")
print(f"Hardware sampling time: {dt/nano} ns")
print(f"Qubit frequency estimate: {qubit_freq_est/giga} GHz")
IBM offline parameters
Qubit: 0
Hardware sampling time: 0.2222222222222222 ns
Qubit frequency estimate: 4.9744443257426045 GHz
Creation and verification of dephasing-robust Boulder Opal pulses
We use Boulder Opal to perform optimizations to achieve an X-gate operation on the qubit. The total Hamiltonian of the driven quantum system is:
\[H(t) = \left(1 + \beta_\gamma (t)\right)H_c(t) + \eta(t) \sigma_z,\]where $\beta_\gamma(t)$ is a fractional time-dependent amplitude fluctuation process, $\eta(t)$ is a small stochastic slowly-varying dephasing noise process, and $H_c(t)$ is the control term given by
\begin{align} H_c(t) = & \frac{1}{2}\left(\gamma^*(t) \sigma_- + \gamma(t) \sigma_+ \right) \\ = & \frac{1}{2}\left(I(t) \sigma_x + Q(t) \sigma_y \right). \end{align}
Here, $\gamma(t) = I(t) + i Q(t)$ is the time-dependent complex-valued control pulse waveform and $\sigma_k$, $k=x,y,z$, are the Pauli matrices.
Creating dephasing-robust pulses
We begin by defining basic operators, constants, and optimization constraints. The latter are defined with the hardware backend characteristics in mind. The duration of each segment of the piecewise constant control pulse, for example, is required to be an integer multiple of the backend sampling time, dt
. This guarantees that Boulder Opal pulses can be properly implemented given the time resolution of the hardware.
sigma_z = np.array([[1.0, 0.0], [0.0, -1.0]], dtype=complex)
sigma_x = np.array([[0.0, 1.0], [1.0, 0.0]], dtype=complex)
sigma_y = np.array([[0.0, -1.0j], [1.0j, 0.0]], dtype=complex)
sigma_p = np.array([[0.0, 0.0], [1.0, 0.0]], dtype=complex)
X_gate = np.array([[0.0, 1.0], [1.0, 0.0]], dtype=complex)
X90_gate = np.array([[1.0, -1j], [-1j, 1.0]], dtype=complex) / np.sqrt(2)
number_of_segments = {}
number_of_optimization_variables = {}
cutoff_frequency = {}
segment_scale = {}
duration_int = {}
duration = {}
scheme_names = ["Square", "Q-CTRL"]
rabi_rotation = np.pi
number_of_runs = 20
omega_max = 2 * np.pi * 8.5e6
I_max = omega_max / np.sqrt(2)
Q_max = omega_max / np.sqrt(2)
# Dephasing-robust pulse parameters
number_of_optimization_variables["Q-CTRL"] = 64
number_of_segments["Q-CTRL"] = 256
segment_scale["Q-CTRL"] = 5
cutoff_frequency["Q-CTRL"] = omega_max * 2
duration["Q-CTRL"] = number_of_segments["Q-CTRL"] * segment_scale["Q-CTRL"] * dt
In the following cell, we create a band-limited robust pulse using a sinc filter to suppress high frequencies in the pulse waveform. The optimization setup involves defining two piecewise constant signals for the $I$ and $Q$ pulses, convolving them with the sinc filter of specified cutoff frequency, and finally discretizing the output signal into the desired number of segments (256 in this case). You can find more details on how to set up optimizations using filters in this user guide.
Note that the execution cells in this notebook contain an option (set at the initialization cell) to run all the commands from scratch or to use previously saved data.
if use_saved_data == False:
robust_dephasing_controls = {}
graph = qctrl.create_graph()
# Create I & Q variables
I_values = graph.optimization_variable(
count=number_of_optimization_variables["Q-CTRL"],
lower_bound=-I_max,
upper_bound=I_max,
)
Q_values = graph.optimization_variable(
count=number_of_optimization_variables["Q-CTRL"],
lower_bound=-Q_max,
upper_bound=Q_max,
)
# Anchor ends to zero with amplitude rise/fall envelope
time_points = np.linspace(
-1.0, 1.0, number_of_optimization_variables["Q-CTRL"] + 2
)[1:-1]
envelope_function = 1.0 - np.abs(time_points)
I_values = I_values * envelope_function
Q_values = Q_values * envelope_function
# Create I & Q signals
I_signal = graph.pwc_signal(values=I_values, duration=duration["Q-CTRL"])
Q_signal = graph.pwc_signal(values=Q_values, duration=duration["Q-CTRL"])
# Apply the sinc filter and re-discretize signals
I_signal = graph.utils.filter_and_resample_pwc(
pwc=I_signal,
cutoff_frequency=cutoff_frequency["Q-CTRL"],
segment_count=number_of_segments["Q-CTRL"],
name="I",
)
Q_signal = graph.utils.filter_and_resample_pwc(
pwc=Q_signal,
cutoff_frequency=cutoff_frequency["Q-CTRL"],
segment_count=number_of_segments["Q-CTRL"],
name="Q",
)
# Create Hamiltonian control terms
I_term = I_signal * graph.pauli_matrix("X") / 2.0
Q_term = Q_signal * graph.pauli_matrix("Y") / 2.0
control_hamiltonian = I_term + Q_term
# Create dephasing noise term
dephasing = graph.constant_pwc_operator(
duration=duration["Q-CTRL"],
operator=graph.pauli_matrix("Z") / 2.0 / duration["Q-CTRL"],
)
# Create infidelity
infidelity = graph.infidelity_pwc(
hamiltonian=control_hamiltonian,
target=graph.target(X_gate),
noise_operators=[dephasing],
name="infidelity",
)
# Calculate optimization
result = qctrl.functions.calculate_optimization(
graph=graph, cost_node_name="infidelity", output_node_names=["I", "Q"]
)
print(f"Cost: {result.cost}")
robust_dephasing_controls["Q-CTRL"] = result.output
scheme_names.append("Q-CTRL")
filename = data_folder / ("robust_dephasing_controls_" + timestr)
save_var(filename, robust_dephasing_controls)
else:
filename = data_folder / "robust_dephasing_controls_demo"
robust_dephasing_controls = load_var(filename)
Here is a ‘primitive’ benchmark pulse: a square pulse with the same duration as the optimized controls to serve as a reference. This duration matches that of the default IBM-Q U3 gate in this backend.
number_of_segments["Square"] = 1
segment_scale["Square"] = 1280
duration["Square"] = segment_scale["Square"] * dt
square_pulse_value = np.array([rabi_rotation / duration["Square"]])
square_sequence = {
"I": qctrl.utils.pwc_arrays_to_pairs(duration["Square"], square_pulse_value),
"Q": qctrl.utils.pwc_arrays_to_pairs(duration["Square"], np.zeros([1])),
}
scheme_names.append("Square")
Now, robust and primitive controls are combined into a single dictionary and plotted.
gate_schemes = {**robust_dephasing_controls, **{"Square": square_sequence}}
for scheme_name, gate in gate_schemes.items():
plot_controls(plt.figure(), gate)
plt.suptitle(scheme_name)
plt.show()
The control plots represent the $I(t)$ and $Q(t)$ components of the control pulses for: Boulder Opal optimizations (top) and the primitive control (bottom). Note that the robust pulse obtained is mostly symmetrical. If desired, time-symmetrical pulses can be enforced within the optimization.
Simulated time evolution of the driven qubit
Another useful way to characterize the pulses is to use the Boulder Opal simulation tool to obtain the time evolution of the system. This way one can observe the effect of arbitrary pulses on the qubit dynamics and check if the pulses are performing the desired operation. The next two cells run the simulations, extract the results, and plot them.
# Boulder Opal simulations
simulated_bloch = {}
gate_times = {}
for scheme_name, control in gate_schemes.items():
simulated_bloch[scheme_name], gate_times[scheme_name] = simulation_coherent(
control, 100
)
Your task calculate_graph (action_id="1141444") has completed.
Your task calculate_graph (action_id="1141445") has completed.
# Plot results
fig, axs = plt.subplots(1, len(gate_schemes.keys()), figsize=(15, 5))
fig.suptitle(f"Time evolution under control pulses for qubit {qubit}", y=1.1)
for idx, scheme_name in enumerate(gate_schemes.keys()):
ax = axs[idx]
for meas_basis in bloch_basis:
ax.plot(
gate_times[scheme_name] / nano,
simulated_bloch[scheme_name][meas_basis],
ls=bloch_lines[meas_basis],
label=f"{meas_basis}: simulation",
)
ax.set_title(scheme_name)
ax.set_xlabel("Time (ns)")
axs[0].set_ylabel("Bloch components")
hs, ls = axs[0].get_legend_handles_labels()
fig.legend(handles=hs, labels=ls, loc="center", bbox_to_anchor=(0.5, 1.0), ncol=3)
plt.show()
The simulated Bloch vector components $(x, y, z)$ of the state dynamics for the duration of the robust and primitive pulses.
Robustness characterization with filter functions and quasi-static scan
Filter functions provide a useful way to determine the sensitivity of pulses to different time varying noise mechanisms. The next cells generate filter functions for the Boulder Opal robust pulses and the primitive square pulse under dephasing noise.
sample_count = 4096
frequencies = np.logspace(-2, np.log10(omega_max), 2000)
filter_functions = {}
for scheme_name, control in gate_schemes.items():
durations, I_values, _ = qctrl.utils.pwc_pairs_to_arrays(control["I"])
durations, Q_values, _ = qctrl.utils.pwc_pairs_to_arrays(control["Q"])
shift_I = qctrl.types.filter_function.Shift(
control=[
qctrl.types.RealSegmentInput(duration=duration, value=value)
for duration, value in zip(durations, I_values)
],
operator=sigma_x / 2,
)
shift_Q = qctrl.types.filter_function.Shift(
control=[
qctrl.types.RealSegmentInput(duration=duration, value=value)
for duration, value in zip(durations, Q_values)
],
operator=sigma_y / 2,
)
dephasing_noise = qctrl.types.filter_function.Drift(noise=True, operator=sigma_z)
filter_function = qctrl.functions.calculate_filter_function(
duration=np.sum(durations),
frequencies=frequencies,
sample_count=sample_count,
shifts=[shift_I, shift_Q],
drifts=[dephasing_noise],
)
filter_functions[scheme_name] = filter_function.samples
Your task calculate_filter_function (action_id="1141446") has started. You can use the `qctrl.get_result` method to retrieve previous results.
Your task calculate_filter_function (action_id="1141446") has completed.
Your task calculate_filter_function (action_id="1141447") has started. You can use the `qctrl.get_result` method to retrieve previous results.
Your task calculate_filter_function (action_id="1141447") has completed.
We next verify the pulse robustness against quasi-static noise by comparing the optimized $\pi$-pulses with the primitive pulse while detuning the drive frequency from the qubit frequency: $\Delta = \omega_\text{drive}-\omega_\text{qubit}$. This is tantamount to an effective dephasing term in the Hamiltonian, the strength and effect of which may be measured by scanning the detuning $\Delta$. The scan is calculated in the cell below.
gate_infidelity = {}
scan_array = np.linspace(-np.pi * mega, np.pi * mega, 21)
gate_infidelity["dephasing"] = {
scheme_name: get_detuning_scan(scan_array, control)
for scheme_name, control in gate_schemes.items()
}
Your task calculate_graph (action_id="1141448") has completed.
Your task calculate_graph (action_id="1141449") has completed.
Here’s a plot of the filter functions and the detuning scan.
fig, axs = plt.subplots(1, 2, figsize=(15, 5))
fig.suptitle("Robustness verification", y=1.1)
ax = axs[0]
ax.set_xlabel("Frequency (Hz)")
ax.set_ylabel("Inverse power")
ax.set_title("Filter functions")
ax.set_xlim(np.min(frequencies), np.max(frequencies))
for scheme_name, _data in filter_functions.items():
inverse_powers = np.array([_d.inverse_power for _d in _data])
inverse_power_precisions = np.array([_d.inverse_power_uncertainty for _d in _data])
ax.loglog(
frequencies, inverse_powers, "-", label=scheme_name, color=colors[scheme_name]
)
y_upper = inverse_powers + inverse_power_precisions
y_lower = inverse_powers - inverse_power_precisions
ax.fill_between(
frequencies,
y_lower,
y_upper,
hatch="||",
alpha=0.35,
facecolor="none",
edgecolor=colors[scheme_name],
linewidth=0.0,
)
ax = axs[1]
ax.set_title("Robustness scan")
ax.set_ylabel("Ground state population")
ax.set_xlabel("Detuning (MHz)")
for name, infidelity in gate_infidelity["dephasing"].items():
ax.plot(scan_array * 1e-6 / np.pi, infidelity, label=name, color=colors[name])
hs, ls = axs[0].get_legend_handles_labels()
fig.legend(handles=hs, labels=ls, loc="center", bbox_to_anchor=(0.5, 1.0), ncol=3)
plt.show()
Left plot: Filter functions for the control pulses under dephasing noise. The hallmark of robustness is the significant drop of the filter function values at low frequencies. This feature is clearly demonstrated for the optimized Boulder Opal pulses. Right plot: Theoretical ground state population as a function of the detuning. The Boulder Opal solution shows flat behavior for a large range of detunings as compared to the primitive pulse, indicating robustness against dephasing.
Experimental time evolution of the driven qubit
We now compare the implementation and performance of the pulses discussed previously. A good test to check that the Boulder Opal pulses are being implemented properly in the IBM hardware is to compare the theoretical with the experimental time evolution of the qubit state. The next cell executes the different steps to obtain such a comparison:
- Calibrate pulses for the specific IBM hardware,
- Set up the schedules and run them on IBM-Q,
- Extract and plot the resulting dynamics.
Note that here we use Rabi rate calibration data obtained previously for the ibmq_armonk
device. To learn how to calibrate devices see the How to automate calibration of control hardware user guide.
# Calibrate pulses
rabi_calibration_exp_I = load_var(data_folder / "rabi_calibration_exp_I_demo")
pulse_amp_array = np.linspace(0.1, 1, 10)
pulse_amp_array = np.concatenate((-pulse_amp_array[::-1], pulse_amp_array))
f_amp_to_rabi = interpolate.interp1d(pulse_amp_array, rabi_calibration_exp_I)
amplitude_interpolated_list = np.linspace(-1, 1, 100)
rabi_interpolated_exp_I = f_amp_to_rabi(amplitude_interpolated_list)
f_rabi_to_amp = interpolate.interp1d(
rabi_interpolated_exp_I, amplitude_interpolated_list
)
waveform = {}
for scheme_name, control in gate_schemes.items():
I_values = qctrl.utils.pwc_pairs_to_arrays(control["I"])[1] / (2 * np.pi)
Q_values = qctrl.utils.pwc_pairs_to_arrays(control["Q"])[1] / (2 * np.pi)
A_I_values = np.repeat(f_rabi_to_amp(I_values), segment_scale[scheme_name])
A_Q_values = np.repeat(f_rabi_to_amp(Q_values), segment_scale[scheme_name])
waveform[scheme_name] = A_I_values + 1j * A_Q_values
# Create IBM schedules
ibm_evolution_times = {}
times_int = {}
pulse_evolution_program = {}
time_min = 64
time_step = 64
for scheme_name in gate_schemes.keys():
time_max = int(segment_scale[scheme_name] * number_of_segments[scheme_name])
times_int[scheme_name] = np.arange(time_min, time_max + time_step, time_step)
ibm_evolution_times[scheme_name] = times_int[scheme_name] * dt
if use_IBM == True:
# Refresh backend
backend.properties(refresh=True)
qubit_frequency_updated = backend.properties().qubit_property(qubit, "frequency")[0]
for scheme_name in gate_schemes.keys():
evolution_schedules = []
for meas_basis in bloch_basis:
for time_idx in times_int[scheme_name]:
schedule = pulse.Schedule(
name="Basis_%s_duration_%d" % (meas_basis, time_idx)
)
schedule += Play(Waveform(waveform[scheme_name][:time_idx]), drive_chan)
if meas_basis == "x":
schedule += inst_sched_map.get("u2", [0], P0=0.0, P1=np.pi)
schedule += measure_schedule << schedule.duration
if meas_basis == "y":
schedule += inst_sched_map.get("u2", [0], P0=0.0, P1=np.pi / 2)
schedule += measure_schedule << schedule.duration
if meas_basis == "z":
schedule += measure_schedule << schedule.duration
evolution_schedules.append(schedule)
pulse_evolution_program[scheme_name] = assemble(
evolution_schedules,
backend=backend,
meas_level=2,
meas_return="single",
shots=num_shots_per_point,
schedule_los=[{drive_chan: qubit_frequency_updated}]
* len(times_int[scheme_name])
* 3,
)
if use_saved_data == False:
# Run the jobs
evolution_exp_results = {}
for scheme_name in gate_schemes.keys():
job = backend.run(pulse_evolution_program[scheme_name])
job_monitor(job)
evolution_exp_results[scheme_name] = job.result(timeout=120)
# Extract results
evolution_results_ibm = {}
for scheme_name in gate_schemes.keys():
evolution_basis = {}
for meas_basis in bloch_basis:
evolution_exp_data = np.zeros(len(times_int[scheme_name]))
for idx, time_idx in enumerate(times_int[scheme_name]):
counts = evolution_exp_results[scheme_name].get_counts(
"Basis_%s_duration_%d" % (meas_basis, time_idx)
)
excited_pop = 0
for bits, count in counts.items():
excited_pop += count if bits[::-1][qubit] == "1" else 0
evolution_exp_data[idx] = excited_pop / num_shots_per_point
evolution_basis[meas_basis] = evolution_exp_data
evolution_results_ibm[scheme_name] = evolution_basis
filename = data_folder / ("bloch_vectors_dephasing_" + timestr)
save_var(filename, evolution_results_ibm)
else:
filename = data_folder / "bloch_vectors_dephasing_demo"
evolution_results_ibm = load_var(filename)
else:
filename = data_folder / "bloch_vectors_dephasing_demo"
evolution_results_ibm = load_var(filename)
# Plot results
fig, axs = plt.subplots(1, len(gate_schemes.keys()), figsize=(15, 5))
fig.set_figheight(5)
fig.set_figwidth(16)
fig.suptitle(f"Time evolution under control pulses for qubit {qubit}", y=1.15)
for idx, scheme_name in enumerate(gate_schemes.keys()):
ax = axs[idx]
for meas_basis in bloch_basis:
ax.plot(
gate_times[scheme_name] / nano,
simulated_bloch[scheme_name][meas_basis],
ls=bloch_lines[meas_basis],
color=bloch_colors[meas_basis],
label=f"{meas_basis}: Q-CTRL simulation",
)
ax.plot(
ibm_evolution_times[scheme_name] / nano,
1 - 2 * evolution_results_ibm[scheme_name][meas_basis],
bloch_markers[meas_basis],
color=bloch_colors[meas_basis],
label=f"{meas_basis}: IBM experiments",
)
ax.set_title(scheme_name)
ax.set_xlabel("Time (ns)")
axs[0].set_ylabel("Bloch components")
hs, ls = axs[0].get_legend_handles_labels()
fig.legend(handles=hs, labels=ls, loc="center", bbox_to_anchor=(0.5, 1.02), ncol=3)
plt.show()
Comparison between the Bloch vector components (x, y, z) for the IBM experiments (symbols) and numerical simulation (lines). The smooth and band-limited Boulder Opal pulse waveform successfully implements the desired gate on the IBM hardware.
Experimental robustness verification with a quasi-static scan
To demonstrate the pulse robustness, we perform the experimental detuning scan, similar to what was done for the theoretical controls. As an extra comparison case, added to the scan here is the default X-gate for this IBM backend. The steps to achieve this are combined in the next execution cell, as follows:
- Define detuning sweep interval,
- Create schedules and run jobs on IBM-Q,
- Extract results and produce robustness plot.
scheme_names.append("IBM default X-gate")
# Define detuning sweep interval
detuning_array = np.linspace(
qubit_frequency_updated - 1e6, qubit_frequency_updated + 1e6, 21
)
number_of_runs = 20
num_shots_per_point = 512
# Create schedules and run jobs
if use_IBM == True and use_saved_data == False:
drive_frequencies = [{drive_chan: freq} for freq in detuning_array]
pi_experiment_results = {}
for run in range(number_of_runs):
for scheme_name in scheme_names:
# Default IBM X-gate
if scheme_name == "IBM default X-gate":
pi_schedule = pulse.Schedule(name="IBM default X-gate")
pi_schedule |= inst_sched_map.get(
"u3", [qubit], P0=np.pi, P1=0.0, P2=np.pi
)
else:
pi_schedule = pulse.Schedule(name=scheme_name)
pi_schedule += Play(Waveform(waveform[scheme_name]), drive_chan)
pi_schedule += measure_schedule << pi_schedule.duration
pi_experiment_ibm = assemble(
pi_schedule,
backend=backend,
meas_level=2,
meas_return="single",
shots=num_shots_per_point,
schedule_los=drive_frequencies,
)
job = backend.run(pi_experiment_ibm)
print(job.job_id())
job_monitor(job)
pi_experiment_results[scheme_name + str(run)] = job.result(timeout=120)
# Extract results
detuning_sweep_results = {}
for scheme_name in scheme_names:
qctrl_sweep_values = np.zeros((number_of_runs, len(detuning_array)))
for run in range(number_of_runs):
i = 0
for result in pi_experiment_results[scheme_name + str(run)].results:
counts = result.data.counts["0x1"]
excited_pop = counts / num_shots_per_point
qctrl_sweep_values[run][i] = 1 - excited_pop
i = i + 1
detuning_sweep_results[scheme_name] = qctrl_sweep_values
filename = data_folder / ("detuning_sweep_results_" + timestr)
save_var(filename, detuning_sweep_results)
else:
filename = data_folder / "detuning_sweep_results_demo"
detuning_sweep_results = load_var(filename)
# Robustness plot
fig, ax = plt.subplots(figsize=(10, 5))
for scheme_name, probability in detuning_sweep_results.items():
ax.errorbar(
detuning_array / mega - qubit_frequency_updated / mega,
np.mean(probability, axis=0),
yerr=np.std(probability, axis=0)
/ np.sqrt(len(detuning_sweep_results[scheme_name])),
fmt="-o",
markersize=5,
capsize=3,
label=f"{scheme_name} run on IBM-Q",
color=colors[scheme_name],
)
if scheme_name != "IBM default X-gate":
plt.plot(
detuning_array / mega - qubit_frequency_updated / mega,
gate_infidelity["dephasing"][scheme_name] + 0.075,
linestyle="dashed",
label=f"{scheme_name} simulation",
color=colors[scheme_name],
)
fig.suptitle("Detuning scan")
plt.xlabel("Detuning (MHz)")
plt.ylabel("Ground state population")
plt.legend()
plt.show()
Experimental (symbols) and simulated (lines) ground state populations as a function of the detuning. The Boulder Opal control solution is mostly flat across the detuning scan, showing its robustness against dephasing noise as compared to the primitive square pulse as well as the default IBM X-gate. The theoretical curves have been shifted up by $0.075$ to account for measurement errors in the experimental data.
Other noise sources: Boulder Opal pulse robustness against control error
Besides robustness to dephasing, we next use Boulder Opal to design solutions that are robust against errors in the control amplitudes. The workflow is exactly as in the dephasing case:
- Initialize Boulder Opal optimization parameters,
- Run optimization and extract results,
- Characterize robustness using filter functions,
- Calibrate and implement pulses,
- Verify robustness with an amplitude scan.
Creating pulses robust to control amplitude errors
The next cell contains all Boulder Opal commands to generate pulses for the control scenario defined previously. The only difference to the dephasing case is that the noise_operators
are now time-varying terms of the control Hamiltonian.
number_of_segments = 256
number_of_optimization_variables = 64
segment_scale = 10
duration_int = number_of_segments * segment_scale
duration = duration_int * dt
cutoff_frequency = omega_max * 2.0
scheme_names = []
if use_saved_data == False:
robust_amplitude_controls = {}
graph = qctrl.create_graph()
# Create I & Q variables
I_values = graph.optimization_variable(
count=number_of_optimization_variables, lower_bound=-I_max, upper_bound=I_max
)
Q_values = graph.optimization_variable(
count=number_of_optimization_variables, lower_bound=-Q_max, upper_bound=Q_max
)
# Anchor ends to zero with amplitude rise/fall envelope
time_points = np.linspace(-1.0, 1.0, number_of_optimization_variables + 2)[1:-1]
envelope_function = 1.0 - np.abs(time_points)
I_values = I_values * envelope_function
Q_values = Q_values * envelope_function
# Create I & Q signals
I_signal = graph.pwc_signal(values=I_values, duration=duration)
Q_signal = graph.pwc_signal(values=Q_values, duration=duration)
# Apply the sinc filter and re-discretize signals
I_signal = graph.utils.filter_and_resample_pwc(
pwc=I_signal,
cutoff_frequency=cutoff_frequency,
segment_count=number_of_segments,
name="I",
)
Q_signal = graph.utils.filter_and_resample_pwc(
pwc=Q_signal,
cutoff_frequency=cutoff_frequency,
segment_count=number_of_segments,
name="Q",
)
# Create Hamiltonian control terms
I_term = I_signal * graph.pauli_matrix("X") / 2.0
Q_term = Q_signal * graph.pauli_matrix("Y") / 2.0
control_hamiltonian = I_term + Q_term
# Create dephasing noise term
amplitude = control_hamiltonian
# Create infidelity
infidelity = graph.infidelity_pwc(
hamiltonian=control_hamiltonian,
target=graph.target(X_gate),
noise_operators=[amplitude],
name="infidelity",
)
# Calculate optimization
result = qctrl.functions.calculate_optimization(
graph=graph, cost_node_name="infidelity", output_node_names=["I", "Q"]
)
print(f"Cost: {result.cost}")
robust_amplitude_controls["Q-CTRL"] = result.output
filename = data_folder / ("robust_amplitude_controls_" + timestr)
save_var(filename, robust_amplitude_controls)
else:
filename = data_folder / "robust_amplitude_controls_demo"
robust_amplitude_controls = load_var(filename)
scheme_names.append("Q-CTRL")
# Create square pulse
square_pulse_value = np.array([rabi_rotation / duration])
square_sequence = {
"I": qctrl.utils.pwc_arrays_to_pairs(duration, square_pulse_value),
"Q": qctrl.utils.pwc_arrays_to_pairs(duration, np.zeros(1)),
}
scheme_names.append("Square")
# Combine all pulses into a single dictionary
gate_schemes = {**robust_amplitude_controls, **{"Square": square_sequence}}
# Add IBM default X-gate
scheme_names.append("IBM default X-gate")
for scheme_name, gate in gate_schemes.items():
plot_controls(plt.figure(), gate)
plt.suptitle(scheme_name)
plt.show()
$I(t)$ and $Q(t)$ components of the control pulses for the optimized (top) and primitive (bottom) controls. Note that the duration of the pulses has increased. This is to enable testing of the robustness of the pulses later, by scanning over a large range of control amplitude errors. Therefore using longer pulses (and lower Rabi rates) will help to always stay within the maximum Rabi rate bounds given by the backend.
Time evolution of the driven qubit: simulations and experiment
As for the amplitude robust pulses, we calibrate the control-error-robust pulses using the interpolation obtained previously and run them on the IBM machine to compare the theoretical and experimental time evolution of the qubit state under no added noise. The cell below combines the following steps:
- Pulse calibration,
- Set up IBM schedules,
- Run jobs on IBM-Q and extract results,
- Run Boulder Opal simulations,
- Plot results.
# Pulse calibration
waveform = {}
for scheme_name, control in gate_schemes.items():
I_values = qctrl.utils.pwc_pairs_to_arrays(control["I"])[1] / (2 * np.pi)
Q_values = qctrl.utils.pwc_pairs_to_arrays(control["Q"])[1] / (2 * np.pi)
if scheme_name == "Square":
A_I_values = np.repeat(
f_rabi_to_amp(I_values), segment_scale * number_of_segments
)
A_Q_values = np.repeat(
f_rabi_to_amp(Q_values), segment_scale * number_of_segments
)
else:
A_I_values = np.repeat(f_rabi_to_amp(I_values), segment_scale)
A_Q_values = np.repeat(f_rabi_to_amp(Q_values), segment_scale)
waveform[scheme_name] = A_I_values + 1j * A_Q_values
# Create IBM schedules
time_min = 64
time_max = duration_int
time_step = 128
num_shots_per_point = 2048
times_int = np.arange(time_min, time_max + time_step, time_step)
ibm_evolution_times = times_int * dt
num_time_points = len(times_int)
if use_IBM == True:
for scheme_name in gate_schemes.keys():
evolution_schedules = []
for meas_basis in bloch_basis:
for time_idx in times_int:
schedule = pulse.Schedule(
name="Basis_%s_duration_%d" % (meas_basis, time_idx)
)
schedule += Play(Waveform(waveform[scheme_name][:time_idx]), drive_chan)
if meas_basis == "x":
schedule += inst_sched_map.get("u2", [0], P0=0.0, P1=np.pi)
schedule += measure_schedule << schedule.duration
if meas_basis == "y":
schedule += inst_sched_map.get("u2", [0], P0=0.0, P1=np.pi / 2)
schedule += measure_schedule << schedule.duration
if meas_basis == "z":
schedule += measure_schedule << schedule.duration
evolution_schedules.append(schedule)
pulse_evolution_program[scheme_name] = assemble(
evolution_schedules,
backend=backend,
meas_level=2,
meas_return="single",
shots=num_shots_per_point,
schedule_los=[{drive_chan: qubit_frequency_updated}] * len(times_int) * 3,
)
if use_saved_data == False:
evolution_exp_results = {}
for scheme_name in gate_schemes.keys():
job = backend.run(pulse_evolution_program[scheme_name])
job_monitor(job)
evolution_exp_results[scheme_name] = job.result(timeout=120)
# Extract results
evolution_results_ibm = {}
for scheme_name in gate_schemes.keys():
evolution_basis = {}
for meas_basis in bloch_basis:
evolution_exp_data = np.zeros(len(times_int))
for idx, time_idx in enumerate(times_int):
counts = evolution_exp_results[scheme_name].get_counts(
"Basis_%s_duration_%d" % (meas_basis, time_idx)
)
excited_pop = 0
for bits, count in counts.items():
excited_pop += count if bits[::-1][qubit] == "1" else 0
evolution_exp_data[idx] = excited_pop / num_shots_per_point
evolution_basis[meas_basis] = evolution_exp_data
evolution_results_ibm[scheme_name] = evolution_basis
filename = data_folder / ("bloch_vectors_amplitude_" + timestr)
save_var(filename, evolution_results_ibm)
else:
filename = data_folder / "bloch_vectors_amplitude_demo"
evolution_results_ibm = load_var(filename)
else:
filename = data_folder / "bloch_vectors_amplitude_demo"
evolution_results_ibm = load_var(filename)
# Boulder Opal simulations
simulated_bloch = {}
gate_times = {}
for scheme_name, control in gate_schemes.items():
simulated_bloch[scheme_name], gate_times[scheme_name] = simulation_coherent(
control, 100
)
Your task calculate_graph (action_id="1141450") has completed.
Your task calculate_graph (action_id="1141451") has completed.
# Plot results
fig, axs = plt.subplots(1, len(gate_schemes.keys()), figsize=(15, 5))
fig.suptitle(f"Time evolution under control pulses for qubit {qubit}", y=1.15)
for idx, scheme_name in enumerate(gate_schemes.keys()):
ax = axs[idx]
for meas_basis in bloch_basis:
ax.plot(
gate_times[scheme_name] / nano,
simulated_bloch[scheme_name][meas_basis],
ls=bloch_lines[meas_basis],
color=bloch_colors[meas_basis],
label=f"{meas_basis}: Q-CTRL simulation",
)
ax.plot(
ibm_evolution_times / nano,
1 - 2 * evolution_results_ibm[scheme_name][meas_basis],
bloch_markers[meas_basis],
color=bloch_colors[meas_basis],
label=f"{meas_basis}: IBM experiments",
)
ax.set_title(scheme_name)
ax.set_xlabel("Time (ns)")
axs[0].set_ylabel("Bloch components")
hs, ls = axs[0].get_legend_handles_labels()
fig.legend(handles=hs, labels=ls, loc="center", bbox_to_anchor=(0.5, 1.02), ncol=3)
plt.show()
Time evolution of the Bloch vector components (x, y, z) for simulations (lines) and experimental runs (symbols) under the control-robust pulses.
Experimental robustness verification with a quasi-static scan and filter functions
We again assess the sensitivity of the pulses to the noise using filter functions and quasi-static noise-susceptibility scans. The only difference here is the change in the noise process, which is achieved by replacing the dephasing term by noise on the controls. The filter functions calculated in the cell below use noise in the $I$ and $Q$ components of the controls.
sample_count = 4096
frequencies = np.logspace(-2, np.log10(omega_max), 2000)
filter_functions = {}
for scheme_name, control in gate_schemes.items():
durations, I_values, _ = qctrl.utils.pwc_pairs_to_arrays(control["I"])
durations, Q_values, _ = qctrl.utils.pwc_pairs_to_arrays(control["Q"])
phasors = I_values + Q_values * 1j
drive = qctrl.types.filter_function.Drive(
control=[
qctrl.types.ComplexSegmentInput(duration=duration, value=value)
for duration, value in zip(durations, phasors)
],
operator=sigma_p / 2,
noise=True,
)
filter_function = qctrl.functions.calculate_filter_function(
duration=np.sum(durations),
frequencies=frequencies,
sample_count=sample_count,
drives=[drive],
)
filter_functions[scheme_name] = filter_function.samples
Your task calculate_filter_function (action_id="1141452") has started. You can use the `qctrl.get_result` method to retrieve previous results.
Your task calculate_filter_function (action_id="1141452") has completed.
Your task calculate_filter_function (action_id="1141453") has started. You can use the `qctrl.get_result` method to retrieve previous results.
Your task calculate_filter_function (action_id="1141453") has completed.
We also demonstrate robustness by performing a quasi-static control-error scan, both numerically and experimentally. The steps follow the same structure as in the dephasing case and are combined in the next cell:
- Define control error range,
- Create schedules and run them on IBM-Q,
- Create Boulder Opal theoretical control error scan.
scheme_names = ["Q-CTRL", "Square", "IBM default X-gate"]
# Define control amplitude error range
amp_error = 0.25
amp_scalings = np.linspace(1 - amp_error, amp_error + 1, 21)
# Define schedules and run jobs
if use_IBM == True and use_saved_data == False:
IBM_default_X_gate_schedule = inst_sched_map.get(
"u3", [qubit], P0=np.pi, P1=0.0, P2=np.pi
)
pi_experiment_results = {}
for run in range(number_of_runs):
for scheme_name in scheme_names:
pi_schedule = []
for amp_scaling in amp_scalings:
if scheme_name == "IBM default X-gate":
this_schedule = pulse.Schedule(name=scheme_name)
for idx, inst in enumerate(
IBM_default_X_gate_schedule.instructions
):
if inst[1].name is None:
this_schedule |= pulse.Schedule(inst)
else: # if "X90" in inst[1].name:
Ival = np.real(inst[1].instructions[0][1].pulse.samples)
Qval = np.imag(inst[1].instructions[0][1].pulse.samples)
phasors = Ival + 1j * Qval
this_schedule |= (
Play(
Waveform(phasors * amp_scaling), inst[1].channels[0]
)
<< this_schedule.duration
)
this_schedule += measure_schedule << this_schedule.duration
pi_schedule.append(this_schedule)
else:
this_schedule = pulse.Schedule(name=scheme_name)
this_schedule += (
Play(Waveform(waveform[scheme_name] * amp_scaling), drive_chan)
<< 1
)
this_schedule += measure_schedule << duration_int
pi_schedule.append(this_schedule)
num_shots_per_point = 512
pi_experiment_ibm = assemble(
pi_schedule,
backend=backend,
meas_level=2,
meas_return="single",
shots=num_shots_per_point,
schedule_los=[{drive_chan: qubit_frequency_updated}]
* len(amp_scalings),
)
job = backend.run(pi_experiment_ibm)
print(job.job_id())
job_monitor(job)
pi_experiment_results[scheme_name + str(run)] = job.result(timeout=120)
# Extracting results
amplitude_sweep_results = {}
for scheme_name in scheme_names:
qctrl_sweep_values = np.zeros((number_of_runs, len(amp_scalings)))
for run in range(number_of_runs):
i = 0
for result in pi_experiment_results[scheme_name + str(run)].results:
counts = result.data.counts["0x1"]
excited_pop = counts / num_shots_per_point
qctrl_sweep_values[run][i] = 1 - excited_pop
i = i + 1
amplitude_sweep_results[scheme_name] = qctrl_sweep_values
filename = data_folder / ("amplitude_sweep_results_" + timestr)
save_var(filename, amplitude_sweep_results)
else:
filename = data_folder / "amplitude_sweep_results_demo"
amplitude_sweep_results = load_var(filename)
# Boulder Opal control amplitude scan
gate_infidelity = {}
amplitude_array = np.linspace(-0.25, 0.25, 21)
gate_infidelity["amplitude"] = {
scheme_name: get_amplitude_scan(amplitude_array, control)
for scheme_name, control in gate_schemes.items()
}
Your task calculate_graph (action_id="1141454") has completed.
Your task calculate_graph (action_id="1141455") has completed.
Next, we plot the numerical filter functions and the control error scan.
fig, axs = plt.subplots(1, 2, figsize=(15, 5))
fig.suptitle("Robustness verification", y=1)
ax = axs[0]
ax.set_xlabel("Frequency (Hz)")
ax.set_ylabel("Inverse power")
ax.set_title("Filter functions")
ax.set_xlim(np.min(frequencies), np.max(frequencies))
for scheme_name, _data in filter_functions.items():
inverse_powers = np.array([_d.inverse_power for _d in _data])
inverse_power_precisions = np.array([_d.inverse_power_uncertainty for _d in _data])
ax.loglog(
frequencies, inverse_powers, "-", label=scheme_name, color=colors[scheme_name]
)
y_upper = inverse_powers + inverse_power_precisions
y_lower = inverse_powers - inverse_power_precisions
ax.fill_between(
frequencies,
y_lower,
y_upper,
hatch="||",
alpha=0.35,
facecolor="none",
edgecolor=colors[scheme_name],
linewidth=0.0,
)
ax.legend()
ax = axs[1]
ax.set_title("Robustness scan")
ax.set_ylabel("Ground state population")
ax.set_xlabel("Control amplitude error")
for scheme_name, probability in amplitude_sweep_results.items():
ax.errorbar(
amp_scalings - 1,
np.mean(probability, axis=0),
yerr=np.std(probability, axis=0)
/ np.sqrt(len(amplitude_sweep_results[scheme_name])),
fmt="-o",
capsize=5,
label=f"{scheme_name} run on IBM-Q",
color=colors[scheme_name],
)
if scheme_name != "IBM default X-gate":
ax.plot(
amp_scalings - 1,
gate_infidelity["amplitude"][scheme_name] + 0.071,
linestyle="dashed",
label=f"{scheme_name} simulation",
color=colors[scheme_name],
)
ax.legend()
plt.show()
Left plot: Filter functions for the control pulses under control noise in the $I$ component. Again, robustness is clearly demonstrated for the optimized Boulder Opal pulses. Right plot: Experimental (symbols) and simulated (lines) ground state populations as a function of the error in the control amplitude (note that the control amplitudes vary between -1 and 1). The Boulder Opal control solution is robust, with a flat response for a large range of control errors as compared to the primitive pulse. The theoretical curves have been shifted up by $0.071$ to account for measurement errors in the experimental data.