Optimal Quantum Control with GRAPE
A quantum gate is not something you simply “apply.” On real hardware, a gate is the end result of a carefully shaped control pulse acting on a qubit that also has its own drift dynamics. Optimal quantum control (GRAPE: GRadient Ascent Pulse Engineering) tunes the control pulse so that the resulting gate fidelity is maximized.
In this article we synthesize a Hadamard gate on a qubit with a nonzero detuning, using GRAPE with exact analytic gradients, and compare it with a naive textbook pulse sequence.
The Problem
We consider a single qubit with the control Hamiltonian
$$
H(t) = \frac{\Delta}{2},\sigma_z + \frac{u_x(t)}{2},\sigma_x + \frac{u_y(t)}{2},\sigma_y ,
$$
where $\Delta = 0.8$ is a fixed detuning (the drift term we cannot switch off) and $u_x(t), u_y(t)$ are the two control amplitudes we are free to choose, with $|u_{x,y}| \le 2$.
The goal is to realize the Hadamard gate
$$
U_{\text{target}} = \frac{1}{\sqrt{2}}\begin{pmatrix} 1 & 1 \ 1 & -1 \end{pmatrix}
$$
in a total time $T = 6$.
Discretization
We split $[0, T]$ into $N = 60$ slices of width $\delta t = T/N$, with piecewise-constant controls $u_{x,k}, u_{y,k}$. Each slice has an exact propagator. Define
Then
$$
U_k = e^{-iH_k\delta t} = \cos\theta_k, I - i,\frac{\sin\theta_k}{\omega_k},\mathbf{h}_k\cdot\boldsymbol{\sigma}.
$$
The total gate is
$$
V = U_N U_{N-1}\cdots U_1 .
$$
Fidelity and Gradient
The (phase-insensitive) gate fidelity is
$$
F = \left|\frac{1}{2},\mathrm{Tr}!\left(U_{\text{target}}^{\dagger} V\right)\right|^2 = |g|^2 .
$$
Define the forward product $X_k = U_k\cdots U_1$ and the backward product $B_k = U_N\cdots U_{k+1}$. The derivative with respect to a control $u_{j,k}$ is
$$
\frac{\partial F}{\partial u_{j,k}} = 2,\mathrm{Re}!\left[,g^{*}\cdot\frac{1}{2},\mathrm{Tr}!\left(U_{\text{target}}^{\dagger}, B_k,\frac{\partial U_k}{\partial u_{j,k}}, X_{k-1}\right)\right].
$$
Every derivative costs only a few $2\times 2$ matrix products, so the full gradient is obtained at roughly the cost of one forward and one backward sweep. We hand it to L-BFGS-B together with box constraints on the amplitudes.
The Baseline
For comparison we use a pulse that ignores the detuning. It is a $\pi/2$ rotation about $y$ followed by a $\pi$ rotation about $x$, which is a Hadamard gate up to a global phase when $\Delta = 0$.
Source Code
1 | import numpy as np |
Code Walkthrough
Setup and constants
The Pauli matrices are stacked into a single array PAULI of shape $(3, 2, 2)$. This lets us build $\mathbf{h}\cdot\boldsymbol{\sigma}$ for all time slices with one einsum call. N, T, DELTA, and UMAX define the discretization, the gate time, the drift, and the amplitude bound. The Hadamard gate is written as $(\sigma_x+\sigma_z)/\sqrt{2}$.
slice_unitaries: the fast exponentials
Calling a generic matrix exponential 60 times per cost evaluation would be wasteful. For a $2\times2$ Hamiltonian of the form $\mathbf{h}\cdot\boldsymbol{\sigma}$, the exponential has the closed form shown earlier, so the function computes the vector $\mathbf{h}_k$, its norm $\omega_k$, and the angle $\theta_k$ for all slices at once with NumPy broadcasting. The small constant 1e-14 added to the norm prevents division by zero. The scale argument multiplies the control amplitudes, and delta replaces the detuning. We use both later to probe robustness against pulse-amplitude and detuning errors.
total_unitary and fidelity
total_unitary multiplies the slice propagators in time order, so later slices multiply from the left, which gives $V = U_N\cdots U_1$. fidelity evaluates $|\mathrm{Tr}(U_{\text{target}}^\dagger V)/2|^2$. Using the absolute value makes the measure insensitive to the global phase, which is physically irrelevant.
cost_and_grad: the heart of GRAPE
First we build the lists of forward products X and backward products B. After this, every slice’s contribution to the gradient needs only the local derivative of $U_k$ sandwiched between $B_k$ and $X_{k-1}$.
The derivative of $U_k$ with respect to $h_j$ follows from differentiating $\cos\theta_k$ and $\sin\theta_k/\omega_k$ using $\partial\omega/\partial h_j = h_j/\omega$. The result is the expression for dU. The factor 0.5 accounts for $h_j = u_j/2$. The function returns $1-F$ and the negated gradient, because the optimizer minimizes.
This analytic gradient is what makes the optimization fast. Finite differences would need $2\cdot 2N = 240$ extra cost evaluations per gradient, whereas the analytic version needs one sweep. The script also verifies the gradient numerically at the end, so you can see that it is correct.
Optimization loop
We run five random restarts. Each starts from small random amplitudes and runs L-BFGS-B with box bounds $[-2, 2]$. The callback records the infidelity at every iteration so we can plot the convergence curve for the best run. The tight tolerances (ftol, gtol) let the optimizer keep going until the infidelity is essentially at machine precision.
Baseline pulse
The naive sequence is built directly in the same discretized format. The first third of the time carries a $y$-drive of strength $(\pi/2)/(T/3)$ and the remaining two thirds carry an $x$-drive of strength $\pi/(2T/3)$. Since the rotation angle is amplitude times duration, these correspond to $R_y(\pi/2)$ and $R_x(\pi)$. This is a perfect Hadamard gate when $\Delta = 0$, but the detuning ruins it.
Robustness maps
We evaluate both pulses on a $25\times25$ grid in which the actual detuning is $\Delta + \epsilon_\Delta$ and the actual amplitudes are multiplied by a factor $s$. This costs 625 cheap evaluations per pulse, which takes a fraction of a second thanks to the closed-form propagators.
Bloch-sphere trajectories
bloch_path subdivides each slice into 6 sub-steps. Within a slice the Hamiltonian is constant, so the sub-steps use exactly the same propagator formula with $\delta t/6$. This gives a smooth curve of $(\langle\sigma_x\rangle,\langle\sigma_y\rangle,\langle\sigma_z\rangle)$ starting from $|0\rangle$.
One figure for everything
All six plots are placed in a single figure, so one plt.show() produces everything: convergence, pulse shapes, the Bloch trajectory in 3D, two 3D robustness surfaces, and the time evolution of the Bloch components.
Console Output
=== GRAPE optimisation (Hadamard gate, detuning = 0.80) === restart 1: fidelity = 1.00000000 (iterations: 12) restart 2: fidelity = 1.00000000 (iterations: 11) restart 3: fidelity = 1.00000000 (iterations: 12) restart 4: fidelity = 1.00000000 (iterations: 12) restart 5: fidelity = 1.00000000 (iterations: 12) === Result summary === Naive square pulses : F = 0.002820 GRAPE optimised : F = 1.00000000 Infidelity (GRAPE) : 5.33e-15 Gradient check : analytic = -0.02026719, numeric = -0.02026719 Fidelity at +-0.2 detuning error (GRAPE / naive): 0.81664 / 0.09848
Reading the Results
Console output
All five random restarts converge to a fidelity of essentially $1$. The infidelity of the best pulse is at the level of numerical round-off, far below $10^{-10}$. This tells us that, with $T=6$ and the given amplitude bound, the Hadamard gate is exactly reachable despite the detuning, and that the optimization landscape is benign enough for the method to find a solution from every start.
The naive pulse reaches a fidelity of only about $0.003$. The detuning is large compared with the gate duration, so the uncorrected rotations are tilted far from the ideal axes, and the gate is almost orthogonal to the target. The gradient check shows the analytic and numerical derivatives agreeing to all printed digits, which confirms the derivation.
Panel 1: convergence
The infidelity drops by many orders of magnitude within roughly a dozen iterations. The semilogarithmic scale makes the final plunge to machine precision visible. This quick convergence is the payoff of the exact gradient combined with a quasi-Newton optimizer.
Panel 2: control pulses
The GRAPE pulses are smooth, irregular-looking waveforms that stay within the $\pm 2$ bounds, while the naive pulse consists of two flat blocks. The optimized shape is not intuitive. It continuously steers the qubit while pre-compensating the rotation about $z$ that the detuning would otherwise cause.
Panel 3: Bloch-sphere trajectory (3D)
The Hadamard gate maps $|0\rangle$ to $|+\rangle$, the point $(1,0,0)$ on the sphere. The GRAPE trajectory winds around the sphere and lands on the green star. The naive trajectory starts the same way but is pushed off course by the drift, so it ends far from the target. This picture shows directly why the naive gate fails.
Panels 4 and 5: robustness surfaces (3D)
The vertical axis is the fidelity, and the horizontal axes are the detuning error and the amplitude scale. The GRAPE surface has a high plateau around the nominal point $(0,1)$ with fidelity $1$. The naive surface is low everywhere. At a detuning error of $+0.2$ the GRAPE pulse still keeps a fidelity of about $0.82$, whereas the naive pulse drops to about $0.10$.
An important point is that the GRAPE pulse was optimized only for the nominal parameters. Its fidelity therefore degrades away from the nominal point, and its robustness here is a side effect of a smooth, well-conditioned solution. When high robustness is the goal, the cost function can be averaged over an ensemble of $(\epsilon_\Delta, s)$ values, which is known as robust GRAPE.
Panel 6: Bloch components over time
The GRAPE curves show $\langle\sigma_x\rangle$ rising to $1$ at the final time while $\langle\sigma_y\rangle$ and $\langle\sigma_z\rangle$ return to $0$. The dashed naive curve for $\langle\sigma_x\rangle$ never gets close to the dotted line at $1$.
Conclusion
With about 100 lines of NumPy, we implemented GRAPE with exact analytic gradients, reaching essentially perfect fidelity for a Hadamard gate on a detuned qubit in just a few iterations. The key ingredients were the closed-form $2\times2$ exponentials, the forward and backward propagator products that turn the gradient into one cheap sweep, and a quasi-Newton optimizer with amplitude bounds. The same structure extends directly to multi-qubit systems by swapping the closed-form slice exponentials for matrix exponentials and the gradient formulas for their exact or approximate derivatives.
is the partial trace over the two-dimensional output factor, and the constraint is exactly trace preservation, $\sum_jR_j^\dagger R_j=I_{16}$.
for every pair $(j,k)$ at once. The same function scores every strategy, which keeps the comparison fair: the uncoded qubit uses $R=I$ with $A_k=K_k$, and the coded cases use $16\times2$ operators $A_k$.