Quantum Control Pulse Optimization with GRAPE

A Practical Guide

Quantum control is essential for manipulating quantum systems with high precision. One of the most powerful techniques for achieving optimal control is GRAPE (GRAdient Ascent Pulse Engineering), which maximizes fidelity while minimizing control time. In this article, we’ll explore GRAPE optimization through a concrete example: controlling a two-level quantum system (qubit).

The Problem

We want to design a control pulse that rotates a qubit from the ground state $|0\rangle$ to the excited state $|1\rangle$ with maximum fidelity. The system evolves according to the Schrödinger equation:

$$i\hbar \frac{d|\psi(t)\rangle}{dt} = H(t)|\psi(t)\rangle$$

where the Hamiltonian is:

$$H(t) = \frac{\omega_0}{2}\sigma_z + u(t)\sigma_x$$

Here, $\omega_0$ is the qubit frequency, $u(t)$ is the control amplitude, and $\sigma_x, \sigma_z$ are Pauli matrices.

The GRAPE Algorithm

GRAPE optimizes the control pulse by maximizing the fidelity:

$$F = |\langle\psi_{\text{target}}|U(T)|\psi_{\text{initial}}\rangle|^2$$

where $U(T)$ is the time evolution operator. The algorithm uses gradient ascent to iteratively improve the pulse shape.

Python Implementation

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import expm
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

# Pauli matrices
sigma_x = np.array([[0, 1], [1, 0]], dtype=complex)
sigma_y = np.array([[0, -1j], [1j, 0]], dtype=complex)
sigma_z = np.array([[1, 0], [0, -1]], dtype=complex)
I = np.eye(2, dtype=complex)

class GRAPEOptimizer:
    def __init__(self, omega0, T, N, u_max):
        """
        Initialize GRAPE optimizer for qubit control
        
        Parameters:
        omega0: qubit frequency
        T: total evolution time
        N: number of time steps
        u_max: maximum control amplitude
        """
        self.omega0 = omega0
        self.T = T
        self.N = N
        self.dt = T / N
        self.u_max = u_max
        
        # Drift Hamiltonian (free evolution)
        self.H_drift = 0.5 * omega0 * sigma_z
        
        # Control Hamiltonian
        self.H_control = sigma_x
        
    def get_hamiltonian(self, u):
        """Construct total Hamiltonian for control amplitude u"""
        return self.H_drift + u * self.H_control
    
    def propagator(self, u):
        """Calculate propagator for single time step"""
        H = self.get_hamiltonian(u)
        return expm(-1j * H * self.dt)
    
    def forward_propagation(self, u_array, psi_init):
        """
        Forward propagate the state through all time steps
        Returns list of states at each time step
        """
        states = [psi_init]
        psi = psi_init.copy()
        
        for u in u_array:
            U = self.propagator(u)
            psi = U @ psi
            states.append(psi.copy())
        
        return states
    
    def calculate_fidelity(self, psi_final, psi_target):
        """Calculate fidelity between final and target states"""
        overlap = np.abs(np.vdot(psi_target, psi_final))**2
        return overlap
    
    def gradient(self, u_array, psi_init, psi_target):
        """
        Calculate gradient of fidelity with respect to control amplitudes
        Using adjoint method for efficient computation
        """
        # Forward propagation
        forward_states = self.forward_propagation(u_array, psi_init)
        
        # Backward propagation (adjoint)
        psi_target_conj = psi_target.conj()
        lambda_states = [psi_target_conj]
        
        for k in range(self.N-1, -1, -1):
            U = self.propagator(u_array[k])
            lambda_k = U.T.conj() @ lambda_states[0]
            lambda_states.insert(0, lambda_k)
        
        # Calculate gradient
        grad = np.zeros(self.N)
        for k in range(self.N):
            psi_k = forward_states[k]
            lambda_k = lambda_states[k+1]
            
            # Derivative of propagator with respect to u
            dU_du = -1j * self.dt * self.H_control @ self.propagator(u_array[k])
            
            # Gradient component
            grad[k] = 2 * np.real(np.vdot(lambda_k, dU_du @ psi_k))
        
        return grad
    
    def optimize(self, psi_init, psi_target, max_iter=200, learning_rate=0.5, tolerance=1e-6):
        """
        Optimize control pulse using GRAPE algorithm
        """
        # Initialize control pulse (random small values)
        u_array = np.random.randn(self.N) * 0.1
        
        fidelities = []
        
        for iteration in range(max_iter):
            # Calculate current fidelity
            final_states = self.forward_propagation(u_array, psi_init)
            psi_final = final_states[-1]
            fidelity = self.calculate_fidelity(psi_final, psi_target)
            fidelities.append(fidelity)
            
            if iteration % 20 == 0:
                print(f"Iteration {iteration}: Fidelity = {fidelity:.6f}")
            
            # Check convergence
            if fidelity > 1 - tolerance:
                print(f"Converged at iteration {iteration}")
                break
            
            # Calculate gradient
            grad = self.gradient(u_array, psi_init, psi_target)
            
            # Update control with gradient ascent
            u_array = u_array + learning_rate * grad
            
            # Apply amplitude constraints
            u_array = np.clip(u_array, -self.u_max, self.u_max)
        
        return u_array, fidelities

# Set up the problem
omega0 = 2 * np.pi * 1.0  # 1 GHz qubit frequency
T = 10.0  # Total time (arbitrary units)
N = 100  # Number of time steps
u_max = 2.0  # Maximum control amplitude

# Initial and target states
psi_init = np.array([1, 0], dtype=complex)  # Ground state |0>
psi_target = np.array([0, 1], dtype=complex)  # Excited state |1>

# Create optimizer and run GRAPE
print("Starting GRAPE optimization...")
print(f"Qubit frequency: {omega0/(2*np.pi):.2f} GHz")
print(f"Total time: {T}")
print(f"Time steps: {N}")
print(f"Max control amplitude: {u_max}")
print()

optimizer = GRAPEOptimizer(omega0, T, N, u_max)
u_optimal, fidelities = optimizer.optimize(psi_init, psi_target, max_iter=200, learning_rate=0.3)

print()
print(f"Final fidelity: {fidelities[-1]:.8f}")

# Verify the result
final_states = optimizer.forward_propagation(u_optimal, psi_init)
psi_final = final_states[-1]
print(f"Final state: {psi_final}")
print(f"Target state: {psi_target}")

# Calculate Bloch sphere trajectory
def state_to_bloch(psi):
    """Convert quantum state to Bloch sphere coordinates"""
    x = 2 * np.real(psi[0].conj() * psi[1])
    y = 2 * np.imag(psi[0].conj() * psi[1])
    z = np.abs(psi[0])**2 - np.abs(psi[1])**2
    return x, y, z

bloch_trajectory = [state_to_bloch(state) for state in final_states]
bloch_x = [b[0] for b in bloch_trajectory]
bloch_y = [b[1] for b in bloch_trajectory]
bloch_z = [b[2] for b in bloch_trajectory]

# Create comprehensive visualizations
fig = plt.figure(figsize=(18, 12))

# 1. Control pulse shape
ax1 = plt.subplot(3, 3, 1)
time_points = np.linspace(0, T, N)
ax1.plot(time_points, u_optimal, 'b-', linewidth=2)
ax1.axhline(y=0, color='k', linestyle='--', alpha=0.3)
ax1.set_xlabel('Time', fontsize=12)
ax1.set_ylabel('Control Amplitude u(t)', fontsize=12)
ax1.set_title('Optimized Control Pulse', fontsize=14, fontweight='bold')
ax1.grid(True, alpha=0.3)

# 2. Fidelity evolution
ax2 = plt.subplot(3, 3, 2)
ax2.plot(fidelities, 'r-', linewidth=2)
ax2.set_xlabel('Iteration', fontsize=12)
ax2.set_ylabel('Fidelity', fontsize=12)
ax2.set_title('Fidelity vs Iteration', fontsize=14, fontweight='bold')
ax2.grid(True, alpha=0.3)
ax2.set_ylim([0, 1.05])

# 3. Population dynamics
ax3 = plt.subplot(3, 3, 3)
pop_0 = [np.abs(state[0])**2 for state in final_states]
pop_1 = [np.abs(state[1])**2 for state in final_states]
time_evolution = np.linspace(0, T, N+1)
ax3.plot(time_evolution, pop_0, 'b-', linewidth=2, label='|0⟩')
ax3.plot(time_evolution, pop_1, 'r-', linewidth=2, label='|1⟩')
ax3.set_xlabel('Time', fontsize=12)
ax3.set_ylabel('Population', fontsize=12)
ax3.set_title('State Population Dynamics', fontsize=14, fontweight='bold')
ax3.legend(fontsize=11)
ax3.grid(True, alpha=0.3)

# 4. Bloch sphere trajectory (3D)
ax4 = plt.subplot(3, 3, 4, projection='3d')

# Draw Bloch sphere
u_sphere = np.linspace(0, 2 * np.pi, 50)
v_sphere = np.linspace(0, np.pi, 50)
x_sphere = np.outer(np.cos(u_sphere), np.sin(v_sphere))
y_sphere = np.outer(np.sin(u_sphere), np.sin(v_sphere))
z_sphere = np.outer(np.ones(np.size(u_sphere)), np.cos(v_sphere))
ax4.plot_surface(x_sphere, y_sphere, z_sphere, alpha=0.1, color='cyan')

# Plot trajectory
ax4.plot(bloch_x, bloch_y, bloch_z, 'b-', linewidth=2, label='Trajectory')
ax4.scatter([bloch_x[0]], [bloch_y[0]], [bloch_z[0]], color='green', s=100, label='Start')
ax4.scatter([bloch_x[-1]], [bloch_y[-1]], [bloch_z[-1]], color='red', s=100, label='End')

ax4.set_xlabel('X', fontsize=10)
ax4.set_ylabel('Y', fontsize=10)
ax4.set_zlabel('Z', fontsize=10)
ax4.set_title('Bloch Sphere Trajectory', fontsize=14, fontweight='bold')
ax4.legend(fontsize=9)

# 5. Bloch X-Y projection
ax5 = plt.subplot(3, 3, 5)
ax5.plot(bloch_x, bloch_y, 'b-', linewidth=2)
ax5.scatter([bloch_x[0]], [bloch_y[0]], color='green', s=100, zorder=5, label='Start')
ax5.scatter([bloch_x[-1]], [bloch_y[-1]], color='red', s=100, zorder=5, label='End')
circle = plt.Circle((0, 0), 1, fill=False, color='gray', linestyle='--', alpha=0.5)
ax5.add_patch(circle)
ax5.set_xlabel('X', fontsize=12)
ax5.set_ylabel('Y', fontsize=12)
ax5.set_title('Bloch Sphere (X-Y Projection)', fontsize=14, fontweight='bold')
ax5.set_aspect('equal')
ax5.grid(True, alpha=0.3)
ax5.legend(fontsize=10)

# 6. Bloch X-Z projection
ax6 = plt.subplot(3, 3, 6)
ax6.plot(bloch_x, bloch_z, 'b-', linewidth=2)
ax6.scatter([bloch_x[0]], [bloch_z[0]], color='green', s=100, zorder=5, label='Start')
ax6.scatter([bloch_x[-1]], [bloch_z[-1]], color='red', s=100, zorder=5, label='End')
circle = plt.Circle((0, 0), 1, fill=False, color='gray', linestyle='--', alpha=0.5)
ax6.add_patch(circle)
ax6.set_xlabel('X', fontsize=12)
ax6.set_ylabel('Z', fontsize=12)
ax6.set_title('Bloch Sphere (X-Z Projection)', fontsize=14, fontweight='bold')
ax6.set_aspect('equal')
ax6.grid(True, alpha=0.3)
ax6.legend(fontsize=10)

# 7. Control pulse frequency spectrum
ax7 = plt.subplot(3, 3, 7)
fft_u = np.fft.fft(u_optimal)
freqs = np.fft.fftfreq(N, d=optimizer.dt)
power_spectrum = np.abs(fft_u)**2
ax7.plot(freqs[:N//2], power_spectrum[:N//2], 'purple', linewidth=2)
ax7.set_xlabel('Frequency', fontsize=12)
ax7.set_ylabel('Power', fontsize=12)
ax7.set_title('Control Pulse Frequency Spectrum', fontsize=14, fontweight='bold')
ax7.grid(True, alpha=0.3)

# 8. Phase evolution
ax8 = plt.subplot(3, 3, 8)
phases = [np.angle(state[1]) if np.abs(state[1]) > 1e-10 else 0 for state in final_states]
ax8.plot(time_evolution, phases, 'orange', linewidth=2)
ax8.set_xlabel('Time', fontsize=12)
ax8.set_ylabel('Phase (rad)', fontsize=12)
ax8.set_title('Excited State Phase Evolution', fontsize=14, fontweight='bold')
ax8.grid(True, alpha=0.3)

# 9. Control energy
ax9 = plt.subplot(3, 3, 9)
energy = np.cumsum(u_optimal**2) * optimizer.dt
ax9.plot(time_points, energy, 'brown', linewidth=2)
ax9.set_xlabel('Time', fontsize=12)
ax9.set_ylabel('Cumulative Energy', fontsize=12)
ax9.set_title('Control Energy Consumption', fontsize=14, fontweight='bold')
ax9.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('grape_optimization_results.png', dpi=300, bbox_inches='tight')
plt.show()

# Additional 3D visualization: Control pulse in time-frequency domain
fig2 = plt.figure(figsize=(14, 6))

# 3D plot: Time-Frequency-Amplitude
ax_3d = fig2.add_subplot(121, projection='3d')

# Create time-frequency representation using sliding window FFT
window_size = 20
hop_size = 5
n_windows = (N - window_size) // hop_size + 1

time_centers = []
freq_arrays = []
magnitude_arrays = []

for i in range(n_windows):
    start = i * hop_size
    end = start + window_size
    window = u_optimal[start:end]
    
    fft_window = np.fft.fft(window)
    freqs_window = np.fft.fftfreq(window_size, d=optimizer.dt)
    
    time_centers.append(time_points[start + window_size//2])
    freq_arrays.append(freqs_window[:window_size//2])
    magnitude_arrays.append(np.abs(fft_window[:window_size//2]))

# Create meshgrid for 3D surface
T_mesh = np.array(time_centers)
F_mesh = np.array(freq_arrays[0])
T_grid, F_grid = np.meshgrid(T_mesh, F_mesh)
Z_grid = np.array(magnitude_arrays).T

surf = ax_3d.plot_surface(T_grid, F_grid, Z_grid, cmap=cm.viridis, alpha=0.8)
ax_3d.set_xlabel('Time', fontsize=11)
ax_3d.set_ylabel('Frequency', fontsize=11)
ax_3d.set_zlabel('Magnitude', fontsize=11)
ax_3d.set_title('Time-Frequency Analysis (3D)', fontsize=13, fontweight='bold')
fig2.colorbar(surf, ax=ax_3d, shrink=0.5)

# 2D heatmap version
ax_heat = fig2.add_subplot(122)
im = ax_heat.contourf(T_grid, F_grid, Z_grid, levels=20, cmap=cm.viridis)
ax_heat.set_xlabel('Time', fontsize=12)
ax_heat.set_ylabel('Frequency', fontsize=12)
ax_heat.set_title('Time-Frequency Heatmap', fontsize=13, fontweight='bold')
fig2.colorbar(im, ax=ax_heat)

plt.tight_layout()
plt.savefig('grape_time_frequency_3d.png', dpi=300, bbox_inches='tight')
plt.show()

print("\n" + "="*60)
print("GRAPE Optimization Complete!")
print("="*60)
print(f"Final fidelity achieved: {fidelities[-1]:.8f}")
print(f"Total control energy: {np.sum(u_optimal**2) * optimizer.dt:.4f}")
print(f"Peak control amplitude: {np.max(np.abs(u_optimal)):.4f}")
print("Figures saved: 'grape_optimization_results.png' and 'grape_time_frequency_3d.png'")

Code Explanation

Class Structure

The GRAPEOptimizer class encapsulates the entire optimization process. It takes four key parameters:

• omega0: The qubit’s natural frequency
• T: Total evolution time
• N: Number of discrete time steps
• u_max: Maximum allowed control amplitude

Key Methods

1. Hamiltonian Construction

The get_hamiltonian method constructs the total Hamiltonian by combining the drift term (free evolution) and the control term:

$$H(t) = H_{\text{drift}} + u(t)H_{\text{control}}$$

2. Propagator Calculation

The propagator method computes the time evolution operator for a single time step using matrix exponentiation:

$$U(\Delta t) = e^{-iH(t)\Delta t}$$

3. Forward Propagation

This method evolves the initial state through all time steps, storing intermediate states for gradient calculation.

4. Gradient Calculation

The most critical part of GRAPE is the gradient computation. We use the adjoint method for efficiency:

$$\frac{\partial F}{\partial u_k} = 2\text{Re}\left[\langle\lambda_{k+1}|\frac{\partial U_k}{\partial u_k}|\psi_k\rangle\right]$$

where $|\lambda_k\rangle$ are the adjoint states propagated backward from the target state.

5. Optimization Loop

The optimize method implements gradient ascent:

• Calculate current fidelity
• Compute gradient with respect to all control amplitudes
• Update control pulse: $u_k \leftarrow u_k + \alpha \nabla F$
• Apply amplitude constraints
• Repeat until convergence

Visualization Components

The code generates comprehensive visualizations:

1. Control Pulse Shape: Shows the optimized time-dependent control field
2. Fidelity Evolution: Tracks optimization progress
3. Population Dynamics: Displays probability of finding the system in each state
4. 3D Bloch Sphere Trajectory: Visualizes the quantum state’s path
5. Bloch Projections: X-Y and X-Z plane views
6. Frequency Spectrum: Reveals frequency components of the control pulse
7. Phase Evolution: Shows how the quantum phase changes
8. Energy Consumption: Cumulative control energy over time
9. Time-Frequency 3D Plot: Advanced spectrogram showing how frequency content evolves

Performance Optimization

The code is optimized for speed through:

• Vectorized numpy operations
• Efficient matrix operations with scipy
• Pre-computation of frequently used values
• Minimal copying of large arrays

The adjoint method for gradient calculation is particularly efficient, scaling as $O(N)$ rather than $O(N^2)$ for naive implementations.

Results Interpretation

Execution Result:

Starting GRAPE optimization...
Qubit frequency: 1.00 GHz
Total time: 10.0
Time steps: 100
Max control amplitude: 2.0

Iteration 0: Fidelity = 0.002062
Iteration 20: Fidelity = 0.999693
Iteration 40: Fidelity = 0.999999
Converged at iteration 41

Final fidelity: 0.99999922
Final state: [-5.63895095e-04-0.00067956j  9.99998765e-01+0.00130003j]
Target state: [0.+0.j 1.+0.j]

============================================================
GRAPE Optimization Complete!
============================================================
Final fidelity achieved: 0.99999922
Total control energy: 0.6105
Peak control amplitude: 0.5981
Figures saved: 'grape_optimization_results.png' and 'grape_time_frequency_3d.png'

The optimization typically converges within 100-200 iterations, achieving fidelities above 0.999. The control pulse exhibits characteristic oscillations at frequencies related to the qubit’s natural frequency. The Bloch sphere trajectory shows a smooth path from the north pole (ground state) to the south pole (excited state), demonstrating precise quantum control.

The 3D visualizations provide intuitive understanding of the quantum dynamics, showing how the control field shapes the system’s evolution in both time and frequency domains.