Pulse Smoothing with GEOPE¤
This notebook demonstrates how to use GEOPE (Geodesic Pulse Engineering) to synthesise piecewise-constant quantum control pulses and then smooth them using two different strategies:
- Geometric smoothing — iteratively projects parameter updates into the null space of the Jacobian so that fidelity is preserved while neighbouring pulse segments are brought closer together.
- Gaussian filtering — a simple signal-processing baseline that convolves the pulse with a Gaussian kernel, then re-optimises to recover lost fidelity.
We target the controlled-Z (CZ) gate on two qubits under an Ising-type Hamiltonian with time-dependent transverse-field controls.
Setup¤
Load the required libraries. GEOPE is built on JAX for automatic differentiation and JIT compilation. We also use SciPy for signal processing and Matplotlib (with the scienceplots style) for publication-quality figures.
import sys
sys.path.append("../")
from geope.geope import Geope
from geope.gecko import Gecko
from geope.parameters import Parameters
from geope.utils import History
from geope.utils import construct_full_pauli_basis, construct_restricted_pauli_basis, qft_unitary
import os
import jax.numpy as jnp
import jax
import numpy as np
import scipy.linalg as spla
import scipy.optimize as spo
import scipy.signal as sps
import re
import matplotlib.pyplot as plt
# Prefer the publication 'science' style when scienceplots is installed;
# otherwise fall back to matplotlib's default so the notebook runs anywhere.
try:
import scienceplots # noqa: F401
plt.style.use('science')
except Exception:
pass
fontsize_title=24
fontsize_label=26
fontsize_legend=18
fontsize_ticks=22
CZ with \(h_1(t)\) control¤
With Hamiltonian model: \(H = \sigma^z_1 \sigma^z_2 + h_1(t)\sigma^x_1\)
Solve with GEOPE and perform repeated Geometric smoothing steps (fewer steps would be sufficient)¤
We construct the CZ gate as a 2-qubit target unitary and set up the Lie algebra bases:
- Full basis — the complete 15-element Pauli basis for 2 qubits.
- Drift basis — the fixed \(\sigma^z_1 \sigma^z_2\) interaction that is always on.
- Projected (control) basis — only \(\sigma^x_1\), the single controllable field \(h_1(t)\).
Two Geope optimisers are created: one for geometric smoothing and one for the Gaussian filter comparison. The optimisation loop runs GEOPE for up to 6 rounds, doubling the number of piecewise-constant segments and applying geometric smoothing after each round until the pulse-difference cost drops below the tolerance.
nqubits = 2
target = np.array([[1,0,0,0],[0,1.j,0,0],[0,0,1.j,0],[0,0,0,1]])
gates = 4
basis = construct_full_pauli_basis(nqubits)
max_step_size = 0.1
drift = construct_restricted_pauli_basis(nqubits, ['zz'])
restriction = construct_restricted_pauli_basis(nqubits, {1: ['x']})
drift_parameters = [1.]
params = Parameters(basis=basis,
projected_basis=restriction,
drift_basis=drift,
drift_values=drift_parameters,
target=target,
piecewise_steps=gates,
init_spread=0.1)
opt_geope = Geope(params,
verbose=True,
history=History())
c=0
success = False
steps = []
while c < 6 and not success:
opt_geope.optimize(max_steps=100, max_step_size=max_step_size, gram_schmidt_step_size=1.2*max_step_size)
if c==0:
# Warm-start a second optimiser from the current solution for the
# Gaussian-filter comparison later on.
params_gf = Parameters(basis=basis,
projected_basis=restriction,
drift_basis=drift,
drift_values=drift_parameters,
target=target,
piecewise_steps=gates,
init_values=np.array(opt_geope.params.parameters))
opt_geope_gf = Geope(params_gf,
verbose=True,
history=History())
steps.append(opt_geope.history.steps[-1])
success, _ = Gecko(opt_geope.params).smooth(piecewise_steps_multiplier=2, smoothing_rate=0.001, max_smoothing_steps=500, diff_tol=0.00001)
steps.append(opt_geope.history.steps[-1])
c+=1
opt_geope.optimize(max_steps=100, max_step_size=max_step_size, gram_schmidt_step_size=1.2*max_step_size)
print("Finished optimising.")
[54/100] [Fidelity = 0.9999999896440438] A solution!
[500/500] [Fidelity = 0.9999999864423904] Smoothing : cost = 0.8998307133739202 (aim = 1e-05)
[500/500] [Fidelity = 0.9999999153903602] Smoothing : cost = 0.09889156642061901 (aim = 1e-05)
[500/500] [Fidelity = 0.9999990456541067] Smoothing : cost = 0.014138945780945067 (aim = 1e-05)
[1/100] [Fidelity = 0.9999999997512626] A solution!
[500/500] [Fidelity = 0.9999845886293534] Smoothing : cost = 0.0019036364470704757 (aim = 1e-05)
[1/100] [Fidelity = 0.9999999989137855] A solution!
[500/500] [Fidelity = 0.9997573649359299] Smoothing : cost = 0.00024352050174859845 (aim = 1e-05)
[1/100] [Fidelity = 0.9999999986650027] A solution!
[500/500] [Fidelity = 0.9994460171911798] Smoothing : cost = 1.888789989691024e-05 (aim = 1e-05)
[1/100] [Fidelity = 0.9999999849140624] A solution!
Finished optimising.
Gaussian filter¤
As an alternative to geometric smoothing, we apply a Gaussian filter (from scipy.ndimage) to the piecewise-constant pulse. The steps are:
- Subdivide the pulse into many fine segments (factor of 64).
- Pad each control channel with zeros, convolve with a Gaussian kernel (\(\sigma = 8\)), and trim.
- Replace the parameters with the smoothed version and re-evaluate the fidelity before and after filtering.
- Re-optimise with GEOPE to recover any fidelity lost from the filtering.
import scipy.ndimage as spnimag
factor = 64
# Subdivide to a fine grid (no smoothing iterations — just resample the pulse).
Gecko(opt_geope_gf.params).smooth(piecewise_steps_multiplier=factor, smoothing_rate=0.01, max_smoothing_steps=0, diff_tol=0.01)
# Gaussian-filter each controllable strength as an alternative smoothing.
before_params = np.array(opt_geope_gf.params.parameters)
sigma = 8
new_params = []
for i, strength in enumerate(opt_geope_gf.params.parameters.T):
if opt_geope_gf.params.basis.labels[i] in opt_geope_gf.params.projected_basis.labels:
strength = np.concatenate((np.zeros(5000), strength, np.zeros(5000)))
smoothed_strength = spnimag.gaussian_filter1d(strength, sigma=sigma)
new_params.append(smoothed_strength[5000:-5000])
else:
new_params.append(strength)
after_params = np.stack(new_params, axis=1)
# Fidelity before/after the Gaussian filter.
free_params = jnp.array([p[opt_geope_gf.params.proj_drift_indices] for p in before_params]).astype(np.complex128)
fid = opt_geope_gf.params.fid_U_fn(opt_geope_gf.params.compute_U_fn(free_params))
print(f"Before: {fid}")
free_params = jnp.array([p[opt_geope_gf.params.proj_drift_indices] for p in after_params]).astype(np.complex128)
fid = opt_geope_gf.params.fid_U_fn(opt_geope_gf.params.compute_U_fn(free_params))
print(f"After: {fid}")
# Adopt the Gaussian-smoothed pulse as the live state, then re-optimise.
opt_geope_gf.params.parameters = after_params
opt_geope_gf.optimize(max_steps=100, max_step_size=max_step_size, gram_schmidt_step_size=1.2*max_step_size)
[0/0] [Fidelity = 0] Smoothing : cost = inf (aim = 0.01)
Before: 0.9999999896440391
After: 0.9959626911167784
[2/100] [Fidelity = 0.9999999990273412] A solution!
<geope.parameters.Parameters at 0x12a640cd0>
Plot results¤
We compare the two smoothing strategies side-by-side. The left panel shows the original piecewise pulse (black) overlaid with the Gaussian-filtered result (blue). The right panel shows the same original pulse overlaid with the geometrically smoothed result (green). Both methods reduce high-frequency discontinuities, but geometric smoothing preserves fidelity by construction.
alpha=0.05
y_limits = [-1.5,8.5]
geometric_smoothing_colour='green'
gf_smoothing_colour='blue'
orig_scaling = opt_geope.history.parameters[steps[0]].shape[0]
smoothed_scaling = opt_geope.params.parameters.shape[0]
gf_scaling = opt_geope_gf.params.parameters.shape[0]
orig_gates = [float(i)/orig_scaling for i in range(-1, opt_geope.history.parameters[steps[0]].shape[0]+1)]
smoothed_gates = [float(i)/smoothed_scaling for i in range(-1, opt_geope.params.parameters.shape[0]+1)]
gf_gates = [float(i)/gf_scaling for i in range(-1, opt_geope_gf.params.parameters.shape[0]+1)]
fig, axs = plt.subplots(1, 2, figsize=(10, 2.5), sharex=True)
for i, strength in enumerate(opt_geope.history.parameters[steps[0]].T):
strength = [s*orig_scaling for s in strength]
if opt_geope.params.basis.labels[i] in opt_geope.params.projected_basis.labels:
axs[0].step(orig_gates,
np.append(np.array([0]), np.append(strength,np.array([0]))),
where='post',
label=opt_geope.params.basis.plot_labels[i],
color='black')
axs[0].fill_between(orig_gates,
np.append(np.array([0]), np.append(strength,np.array([0]))),
step='post',
color='black',
alpha=0.05)
for i, strength in enumerate(opt_geope_gf.params.parameters.T):
strength = [s*gf_scaling for s in strength]
if opt_geope_gf.params.basis.labels[i] in opt_geope_gf.params.projected_basis.labels:
axs[0].step(gf_gates,
np.append(np.array([0]), np.append(strength,np.array([0]))),
where='post',
label=opt_geope_gf.params.basis.plot_labels[i],
color=gf_smoothing_colour)
axs[0].fill_between(gf_gates,
np.append(np.array([0]), np.append(strength,np.array([0]))),
step='post',
color=gf_smoothing_colour,
alpha=alpha)
for i, strength in enumerate(opt_geope.history.parameters[steps[0]].T):
strength = [s*orig_scaling for s in strength]
if opt_geope.params.basis.labels[i] in opt_geope.params.projected_basis.labels:
axs[1].step(orig_gates,
np.append(np.array([0]), np.append(strength,np.array([0]))),
where='post',
label=opt_geope.params.basis.plot_labels[i],
color='black')
axs[1].fill_between(orig_gates,
np.append(np.array([0]), np.append(strength,np.array([0]))),
step='post',
color='black',
alpha=0.05)
for i, strength in enumerate(opt_geope.params.parameters.T):
strength = [s*smoothed_scaling for s in strength]
if opt_geope.params.basis.labels[i] in opt_geope.params.projected_basis.labels:
axs[1].step(smoothed_gates,
np.append(np.array([0]), np.append(strength,np.array([0]))),
where='post',
label=opt_geope.params.basis.plot_labels[i],
color=geometric_smoothing_colour)
axs[1].fill_between(smoothed_gates,
np.append(np.array([0]), np.append(strength,np.array([0]))),
step='post',
color=geometric_smoothing_colour,
alpha=alpha)
for ax in axs:
# ax.legend(loc='upper right')
ax.grid(False)
ax.set_xlim(-0.01,smoothed_gates[-1]+0.008)
# ax.set_ylim(y_limits[0],y_limits[1])
axs[0].set_title(f"(a) Gaussian filter")
axs[1].set_title(f"(b) Geometric smoothing")
axs[0].set_xlabel("Unit time, $t$")
axs[1].set_xlabel("Unit time, $t$")
axs[0].set_ylabel(f"$h_1(t)$")
axs[1].set_ylabel(f"$h_1(t)$")
plt.tight_layout()
plt.show()
CZ with \(h_1(t)\) and \(h_2(t)\) control¤
With Hamiltonian model: \(H = \sigma^z_1 \sigma^z_2 + h_1(t)\sigma^x_1 + h_2(t)\sigma^x_2\)
Solve with GEOPE and then repeated geometric smoothing¤
We now repeat the experiment with two independent controls, \(h_1(t)\sigma^x_1\) and \(h_2(t)\sigma^x_2\), while keeping the same \(\sigma^z_1\sigma^z_2\) drift. Having more control degrees of freedom generally makes the optimisation landscape more favourable for convergence.
nqubits = 2
target = np.array([[1,0,0,0],[0,1.j,0,0],[0,0,1.j,0],[0,0,0,1]])
gates = 4
basis = construct_full_pauli_basis(nqubits)
max_step_size = 0.1
drift = construct_restricted_pauli_basis(nqubits, ['zz'])
restriction = construct_restricted_pauli_basis(nqubits, {1: ['x'], 2: ['x']})
drift_parameters = [1.]
params = Parameters(basis=basis,
projected_basis=restriction,
drift_basis=drift,
drift_values=drift_parameters,
target=target,
piecewise_steps=gates,
init_spread=0.01)
opt = Geope(params,
verbose=True,
history=History())
c=0
success = False
steps_2 = []
while c < 6 and not success:
opt.optimize(max_steps=100, max_step_size=max_step_size, gram_schmidt_step_size=1.2*max_step_size)
success, _ = Gecko(opt.params).smooth(piecewise_steps_multiplier=2, smoothing_rate=0.001, max_smoothing_steps=500, diff_tol=0.00001)
steps_2.append(opt.history.steps[-1])
c+=1
opt.optimize(max_steps=100, max_step_size=max_step_size, gram_schmidt_step_size=1.2*max_step_size)
print("Finished optimising.")
[50/100] [Fidelity = 0.9999999999989022] A solution!
[500/500] [Fidelity = 0.9999999998647647] Smoothing : cost = 0.7182509435122655 (aim = 1e-05)
[500/500] [Fidelity = 0.9999999917004302] Smoothing : cost = 0.07466545252254528 (aim = 1e-05)
[500/500] [Fidelity = 0.9999995415604149] Smoothing : cost = 0.009589850630987084 (aim = 1e-05)
[1/100] [Fidelity = 0.9999999999856406] A solution!
[500/500] [Fidelity = 0.9999897419316898] Smoothing : cost = 0.0012013125811606627 (aim = 1e-05)
[1/100] [Fidelity = 0.9999999996706552] A solution!
[500/500] [Fidelity = 0.9997881084148444] Smoothing : cost = 0.00014533717482892892 (aim = 1e-05)
[1/100] [Fidelity = 0.9999999901373964] A solution!
[245/500] [Fidelity = 0.9992115104384376] Smoothing : cost = 1.7404440013375502e-05 (aim = 1e-05)
Plot¤
Each subplot shows one control channel. The black step function is the initial (unsmoothed) pulse; the green step function is the geometrically smoothed result. Both panels share the same time axis.
alpha = 0.05
y_limits = [-8.5,3]
init_step = steps_2[0]
smoothed_step = -1
init_scaling = opt.history.parameters[init_step].shape[0]
smoothed_scaling = opt.history.parameters[smoothed_step].shape[0]
gates = [float(i)/init_scaling for i in range(-1,opt.history.parameters[init_step].shape[0]+1)]
smoothed_gates = [float(i)/smoothed_scaling for i in range(-1,opt.history.parameters[smoothed_step].shape[0]+1)]
fig, axs = plt.subplots(1, 2, figsize=(10, 2.5), sharex=True)
a=0
for i, strength in enumerate(opt.history.parameters[init_step].T):
strength = [s*init_scaling for s in strength]
if opt.params.basis.labels[i] in opt.params.projected_basis.labels:
axs[a].step(gates,
np.append(np.array([0]), np.append(strength,np.array([0]))),
where='post',
label=opt.params.basis.plot_labels[i],
color="black")
axs[a].fill_between(gates,
np.append(np.array([0]), np.append(strength,np.array([0]))),
step='post',
color="black",
alpha=alpha)
a+=1
a=0
for i, strength in enumerate(opt.history.parameters[smoothed_step].T):
strength = [s*smoothed_scaling for s in strength]
if opt.params.basis.labels[i] in opt.params.projected_basis.labels:
axs[a].step(smoothed_gates,
np.append(np.array([0]), np.append(strength,np.array([0]))),
where='post',
label=opt.params.basis.plot_labels[i],
color="green")
axs[a].fill_between(smoothed_gates,
np.append(np.array([0]), np.append(strength,np.array([0]))),
step='post',
color="green",
alpha=alpha)
a+=1
for ax in axs:
# ax.legend(loc='upper right')
ax.grid(False)
ax.set_xlim(-0.01,gates[-1]+0.008)
# ax.set_ylim(y_limits[0],y_limits[1])
axs[0].set_xlabel("Unit time, $t$")
axs[1].set_xlabel("Unit time, $t$")
axs[0].set_ylabel(f"$h_{{1}}(t)$")
axs[1].set_ylabel(f"$h_{{2}}(t)$")
plt.tight_layout()
plt.show()