Control hardware: pulse optimization under realistic experimental constraints

Highly flexible optimizer for hardware-limited signal generation and non-ideal control lines

The Q-CTRL optimization engine provides a large, modular collection of configuration options in order to support a wide variety of signal generators, control lines and quantum hardware. In this guide, we provide an overview of the Q-CTRL optimization engine, then provide a set of examples that demonstrate how this engine can be tailored to model realistic constraints on the signal generation, including transients and bandwidth limits. In the first section, we show how to deal with these issues by incorporating filters in the control signal. This can also be used to design a pulse where the natural filtering properties of the control line are known. In the second section, we present a simpler approach to deal with band-limited signals by constraining the rate of change in the control pulse. Finally, in the last section we show how to create optimizations using arbitrary basis functions.

This notebook contains preview features that are not currently available in BOULDER OPAL. Please contact us if you want a live demonstration for your hardware. Preview features are included through the package qctrlcore, which is not publicly available: cells with qctrlcore will not be executable.

Table of Contents

  1. Overview: the Q-CTRL optimization engine

  2. Smoothing control pulses using linear filters

    1. Low-pass sinc filter
    2. RC filter
  3. Optimization with band-limited pulses via bounded slew rates

  4. Optimization using arbitrary basis functions

Imports and initialization

# Essential imports
import numpy as np
import qctrlcore as qc
import tensorflow as tf

# Plotting imports and setup
from qctrlvisualizer import (
    get_qctrl_style,
    plot_controls,
    plot_smooth_controls,
)
import matplotlib.pyplot as plt
from matplotlib.ticker import EngFormatter
formatter = EngFormatter(unit='')
plt.style.use(get_qctrl_style())

# Define standard matrices
identity = np.array([[1, 0], [0, 1]])
sigma_x = np.array([[0, 1], [1, 0]])
sigma_y = np.array([[0, -1j], [1j, 0]])
sigma_z = np.array([[1, 0], [0, -1]])
sigma_m = np.array([[0, 1], [0, 0]])

# Define helper functions
def segments_to_dicts(durations, values):
    return [{'duration': d, 'value': v} for d, v in zip(durations, values)]
def samples_to_dicts(times, values):
    return [{'time': d, 'value': v} for d, v in zip(times, values)]

1. Overview of the Q-CTRL optimization engine

Standard optimal control techniques such as GRAPE have seen widespread success in recent years, but their flexibility is typically limited by enforcing very specific Hamiltonian structures. Specifically, a standard GRAPE implementation requires that the Hamiltonian be expressed as a sum of terms that are directly and independently controllable. This assumption, while valid in some cases, is not adequate for many systems. For example, if there is non-linear dependence between controllable parameters and Hamiltonian terms, or alternatively if multiple Hamiltonian terms depend on the same underlying controls, the system cannot be coerced into the form required by GRAPE. Similarly, with standard techniques, creating optimized controls that exhibit temporal structure such as time symmetry or bandwidth limits is possible only by implementing changes in the algorithm itself (as opposed to simply changing its inputs).

The Q-CTRL optimization engine resolves these issues by allowing the user to express their system Hamiltonians as almost-arbitrary functions of the controllable parameters. That is, rather than enforcing a specific Hamiltonian structure and allowing the user to define the values of specific operators in that structure, the Q-CTRL optimization engine allows the user to define the entire Hamiltonian. This flexibility is enabled by leveraging the TensorFlow automatic differentiation library. Using simple convenience functions provided by Q-CTRL, the user creates a representation of their system that defines a map from controllable parameters to Hamiltonian. The underlying structure of this map is a TensorFlow graph, which can be efficiently evaluated and differentiated. Once constructed, this mapping (or graph) is then inserted into the optimization cost function calculation. The resulting optimized controls thus achieve the desired objectives, and do so within the constraints imposed by the user-defined Hamiltonian structure.

Speaking concretely, in the code samples in this notebook, the create_optimization_specification functions create these user-defined graphs. The variable_factory object is used to create controllable, or optimizable, parameters (representing the available controls), while the returned OptimizationSpecification describes the cost function to minimize. All steps in between are entirely customizable by the user, although a typical flow is to:

  1. create "signals", or scalar-valued functions of time, which may be non-linear, have enforced temporal structure such as time symmetry, or more generally may depend arbitrarily on the controllable parameters;
  2. create "operators", or operator-valued functions of time, by modulating constant operators by signals;
  3. combine the operators into a single Hamiltonian operator;
  4. calculate the optimization cost function (typically an infidelity) from the Hamiltonian.

Once this graph has been defined, it is passed to the optimize function in order to run the optimization and produce the optimized controls.

Return to Table of Contents

2. Smoothing control pulses using linear filters

In this section we present examples showing how linear time-invariant filters may be incorporated into an optimization. This is the relevant scenario when trying to model and design pulses that go through control lines in which electronic filters are naturally present. Failing to take the filter into account during the optimization can lead to poor results, since in that case the system used for the optimization does not accurately model reality. Adding filters is also an appropriate strategy to deliberately smooth out control pulses that would otherwise be badly implemented due to the characteristics of the signal generators themselves, such as bandwidth limits or specific resonances. To exemplify the use of filters, we consider a basic single-qubit system described by the Hamiltonian:

\begin{align*} H(t) &= \frac{1}{2} L(\alpha)(t)\sigma_{x} + \beta(t) \sigma_{z} \,, \end{align*}

where $\alpha(t)$ is a real time-dependent pulse, $\beta(t)$ is a small, slowly-varying stochastic dephasing noise process, and $ L$ is the filter applied to the $\alpha(t)$ pulse. The effect of the filter is to transform the control signal before it reaches the quantum system, via convolution with the filter impulse response.

The Q-CTRL optimization engine provides convenience functions for transforming a piecewise-constant signal by a linear time-invariant filter, and for producing impulse response functions representing such filters.

Low-pass sinc filter

In this example, we use an in-built Q-CTRL function to produce the sinc filter to smooth the control pulse. From the optimizations we output two sets of signals: the raw piecewise-constant signal and its smoothed version. The latter is the filtered signal that actually reaches the quantum system and performs the optimized gate.

# Define physical constraints
alpha_max = 2*np.pi * 8.5e6 #Hz
sinc_cutoff_frequency = 2*np.pi * 48e6 #Hz
segment_count = 60
duration = 250e-9 #s

def create_optimization_specification(variable_factory):
    # Create alpha(t) signal
    alpha = qc.create_pwc_signal(
        values=variable_factory.create_bounded(segment_count, -alpha_max, alpha_max),
        duration=duration)
    
    # Create filtered signal
    sinc_filter = qc.create_sinc_integral_function(sinc_cutoff_frequency)
    alpha_filtered = qc.convolve_pwc(alpha, sinc_filter)
    
    # Create control term
    control_sigma_x = qc.create_stf_operator(operator=qc.HermitianOperator(sigma_x),
                                             signal=alpha_filtered)
    
    # Create dephasing noise term
    dephasing = qc.create_constant_stf_operator(operator=qc.HermitianOperator(sigma_z/duration))

    # Create infidelity
    sample_times = np.linspace(0, duration, 150)
    infidelity = qc.create_infidelity_stf(sample_times=sample_times,
                                          hamiltonian=qc.get_stf_sum([control_sigma_x]),
                                          noise_operators=[dephasing],
                                          target=qc.Target(sigma_x))
    
    # For demonstration purposes, we also calculate the filter
    # impulse responses by convolving with a delta-like function.
    delta_function = qc.create_pwc_signal(
        values=tf.constant([1000/duration], dtype=tf.float64),
        duration=duration/1000)
    sinc_impulse_response = qc.convolve_pwc(delta_function, sinc_filter)
    filter_sample_times = np.linspace(-duration/2, duration/2, 1000)
    
    return qc.OptimizationSpecification(
        cost=infidelity,
        output_callable=lambda s: [
            {'unfiltered': {'$\\alpha$': segments_to_dicts(alpha.durations, s.run(alpha.values))},
             'filter': ('Sinc impulse response', (filter_sample_times,
                         s.run(sinc_impulse_response(filter_sample_times)))),
             'filtered': {'$L(\\alpha)$': 
                          samples_to_dicts(sample_times, s.run(alpha_filtered(sample_times)))}}])

# Run optimization
result = qc.optimize(create_optimization_specification)
print(f"Cost: {result['cost']}")
      
# Visualize controls, filters and filtered controls
for d in result['output']:
    print("\n\n")
    plot_controls(plt.figure(), d['unfiltered'])
    plt.title("Unfiltered")
    plt.show()
    
    plt.figure(figsize=(7, 2))
    plt.plot(*d['filter'][1])
    plt.xlabel("Time (s)")
    plt.ylabel(f"{d['filter'][0]}\n(1/s)")
    plt.title("Filter")
    ax = plt.gca()
    ax.xaxis.set_major_formatter(formatter)
    ax.yaxis.set_major_formatter(formatter)
    plt.show()
               
    plot_smooth_controls(plt.figure(), d['filtered'])
    plt.title("Filtered")
    plt.show()
Cost: 1.924286990360414e-07



Unfiltered (top) and filtered (bottom) control amplitudes as a function of time. The filter impulse response is shown in the middle panel.

Return to Table of Contents

RC filter

Here we demonstrate how to implement a custom filter using the RC filter as an example. This is done by defining a function that calculates the integral of the filter impulse response. In the RC case, the impulse response is given by $\frac{e^{-t/T}}{T}{\cal {H}}(t)$, with ${\mathcal H}$(t) being the Heaviside step function and $T$ the time constant. The integral from $0$ to $t$ is therefore $(1-e^{-t/T}){\cal H}(t)$ and is directly implemented in the function rc_filter in the cell below. As in the sinc filter case, the filter is incorporated in the optimizer and we output the raw and the filtered versions of the control.

# Define physical constraints
alpha_max = 2*np.pi * 8.5e6 #Hz
nu = 2*np.pi * 6e6 #Hz
rc_time_constant = 20e-9 #s
segment_count = 50
duration = 250e-9 #s

def rc_filter(start, end):
    T = rc_time_constant
    integral = lambda t: tf.where_v2(t >= 0, 1-tf.exp(-t/T), 0.)
    return integral(end) - integral(start)

def create_optimization_specification(variable_factory):
    # Create alpha_1(t) signal
    alpha = qc.create_pwc_signal(
        values=variable_factory.create_bounded(segment_count, -alpha_max, alpha_max),
        duration=duration)
    # Create filtered signals
    alpha_filtered = qc.convolve_pwc(alpha, rc_filter)
    # Create control term
    control_sigma_x = qc.create_stf_operator(operator=qc.HermitianOperator(sigma_x),
                                             signal=alpha_filtered)
    # Create dephasing noise term
    dephasing = qc.create_constant_stf_operator(operator=qc.HermitianOperator(sigma_z/duration))
    # Create infidelity
    sample_times = np.linspace(0, duration, 150)
    infidelity = qc.create_infidelity_stf(sample_times=sample_times,
                                          hamiltonian=qc.get_stf_sum([control_sigma_x]),
                                          noise_operators=[dephasing],
                                          target=qc.Target(sigma_x))
    
    # For demonstration purposes, we also calculate the filter
    # impulse responses by convolving with a delta-like function.
    delta_function = qc.create_pwc_signal(values=tf.constant([1000/duration], dtype=tf.float64),
                                          duration=duration/1000)
    rc_impulse_response = qc.convolve_pwc(delta_function, rc_filter)
    filter_sample_times = np.linspace(-duration/2, duration/2, 1000)
    
    return qc.OptimizationSpecification(
        cost=infidelity,
        output_callable=lambda s: [
            {'unfiltered': {'$\\alpha$': segments_to_dicts(alpha.durations, s.run(alpha.values))},
             'filter': ('RC impulse response', (filter_sample_times,
                         s.run(rc_impulse_response(filter_sample_times)))),
             'filtered': {'$R(\\alpha)$': 
                          samples_to_dicts(sample_times, s.run(alpha_filtered(sample_times)))}}])

# Run optimization
result = qc.optimize(create_optimization_specification)
print(f"Cost: {result['cost']}")

# Visualize controls, filters and filtered controls
for d in result['output']:
    print("\n\n")
    plot_controls(plt.figure(), d['unfiltered'])
    plt.title("Unfiltered")
    plt.show()
    
    plt.figure(figsize=(7, 2))
    plt.plot(*d['filter'][1])
    plt.xlabel("Time (s)")
    plt.ylabel(f"{d['filter'][0]}\n(1/s)")
    plt.title("Filter")
    ax = plt.gca()
    ax.xaxis.set_major_formatter(formatter)
    ax.yaxis.set_major_formatter(formatter)
    plt.show()
               
    plot_smooth_controls(plt.figure(), d['filtered'])
    plt.title("Filtered")
    plt.show()
Cost: 2.185015267315081e-09



Unfiltered (top) and filtered (bottom) control amplitudes as a function of time. The filter impulse response is shown in the middle panel.

Return to Table of Contents

3. Optimization with band-limited pulses via bounded slew rates

In this section, we show how to optimize a system in which the rates of change of the controls are limited. Using this constraint can help to ensure that optimized controls can be reliably implemented on physical hardware (which may be subject to bandwidth limits, for example). We consider a standard single-qubit system subject to dephasing noise:

\begin{align*} H(t) &= \frac{1}{2} \alpha_1(t)\sigma_{x} + \frac{1}{2} \alpha_2(t)\sigma_{z} + \beta(t) \sigma_{z} \,, \end{align*}

where $\alpha_1(t)$ and $\alpha_2(t)$ are real time-dependent pulses and $\beta(t)$ is a small, slowly-varying stochastic dephasing noise process. In this case, we enforce a maximum slew rate constraint on $\alpha_k(t)$ to cap the variation between adjacent segment values.

To implement a bounded slew rate, we create signal values using a Q-CTRL helper function, create_anchored_difference_bounded_variables. This function produces values that are constrained to satisfy the slew rate requirement, and in addition are anchored to zero at the start and end of the gate.

# Define physical constraints
alpha_max = 2*np.pi * 8.5e6 #Hz
max_slew_rate = alpha_max / 10
segment_count = 250
duration = 400e-9 #s

def create_optimization_specification(variable_factory):
    # Create alpha_1(t) signal
    alpha_1_values = qc.create_anchored_difference_bounded_variables(
        variable_factory, segment_count, -alpha_max, alpha_max, max_slew_rate)
    alpha_1 = qc.create_pwc_signal(values=alpha_1_values, duration=duration)
    # Create alpha_2(t) signal
    alpha_2_values = qc.create_anchored_difference_bounded_variables(
        variable_factory, segment_count, -alpha_max, alpha_max, max_slew_rate)
    alpha_2 = qc.create_pwc_signal(values=alpha_2_values, duration=duration)
    
    # Create control terms
    control_x = qc.create_pwc_operator(
        operator=qc.HermitianOperator(sigma_x/2),
        signal=alpha_1)
    control_z = qc.create_pwc_operator(
        operator=qc.HermitianOperator(sigma_z/2),
        signal=alpha_2)
    
    # Create dephasing noise term
    dephasing = qc.create_constant_pwc_operator(operator=qc.HermitianOperator(sigma_z/duration),
                                                duration=duration)
       
    # Create infidelity
    infidelity = qc.create_infidelity(hamiltonian=qc.get_pwc_operator_sum([control_x, control_z]),
                                      noise_operators=[dephasing],
                                      target=qc.Target(sigma_y))
    
    return qc.OptimizationSpecification(cost=infidelity,
        output_callable=lambda s: {
            '$\\alpha_1$': segments_to_dicts(alpha_1.durations, s.run(alpha_1.values)),
            '$\\alpha_2$': segments_to_dicts(alpha_2.durations, s.run(alpha_2.values))})

# Run optimization
result = qc.optimize(create_optimization_specification)
print(f"Cost: {result['cost']}")

# Visualize optimized pulses
plot_controls(plt.figure(), result['output'])
plt.show()
Cost: 2.678665249092355e-10

Pulse amplitudes obtained from a bound slew optimization for $\alpha_1(t)$ (top) and $\alpha_2(t)$ (bottom).

Return to Table of Contents

4. Optimization using arbitrary basis functions

In many situations, using an appropriate set of basis functions can greatly reduce the dimensionality of the search space in an optimization problem. In Chopped RAndom Basis (CRAB) optimization, for example, pulses are defined via optimizable linear combinations from the basis functions set. Traditionally, a randomized Fourier basis is used, although the same technique has also seen success with other bases, for example Slepian functions. In this example, we perform a CRAB optimization (in the Fourier basis) of a qutrit system in which we effect a single-qubit gate while minimizing leakage out of the computational subspace. The system is described by the following Hamiltonian:

\begin{align*} H(t) = & \frac{\chi}{2} (a^\dagger)^2 a^2 + \gamma(t) a + \gamma^*(t) a^\dagger + \frac{\alpha(t)}{2} a^\dagger a \,, \end{align*}

where $\chi$ is the anharmonicity, $\gamma(t)$ and $\alpha(t)$ are, respectively, complex and real time-dependent pulses, and $a = |0 \rangle \langle 1 | + \sqrt{2} |1 \rangle \langle 2 |$.

The Q-CTRL optimization engine provides a convenience function for creating optimizable signals in a Fourier basis, suitable for use in a CRAB optimization. Other bases are supported by the framework, but require the user to manually provide operations that compute the appropriate linear combinations.

# Define standard matrices
a = np.array([[0., 1., 0.],
              [0., 0., np.sqrt(2)],
              [0., 0., 0.]], dtype=np.complex)
ada = np.matmul(a.T, a)
ad2a2 = np.matmul(np.matmul(a.T,a.T), np.matmul(a,a))
hadamard = np.array([[1., 1.,0],
                     [1., -1.,0],
                     [0, 0, np.sqrt(2)]], dtype=np.complex)/np.sqrt(2)
qubit_projector = np.array([[1, 0, 0],
                            [0, 1, 0],
                            [0, 0, 0]], dtype=np.complex)

# Define physical constraints
chi = 2 * np.pi * -300.0 * 1e6 #Hz
gamma_max = 2 * np.pi * 30e6 #Hz
alpha_max = 2 * np.pi * 30e6 #Hz
segment_count = 200
duration = 100e-9 #s

def create_optimization_specification(variable_factory):
    # Create gamma(t) signal in Fourier bases. To demonstrate the full
    # flexibility, we show how to use both randomized and optimizable
    # basis elements. Elements with fixed frequencies may be chosen too.
    gamma_i = qc.create_real_fourier_pwc_signal(variable_factory=variable_factory,
                                                duration=duration,
                                                segments_count=segment_count,
                                                randomized_frequencies_count=10)
    gamma_q = qc.create_real_fourier_pwc_signal(variable_factory=variable_factory,
                                                duration=duration,
                                                segments_count=segment_count,
                                                optimizable_frequencies_count=10)
    gamma = qc.TensorPwc(durations=gamma_i.durations,
                         values=tf.complex(gamma_i.values, gamma_q.values)*gamma_max)
    
    # Create alpha(t) signal.
    alpha = qc.create_real_fourier_pwc_signal(variable_factory=variable_factory,
                                              duration=duration,
                                              segments_count=segment_count,
                                              initial_coefficient_lower_bound=-alpha_max,
                                              initial_coefficient_upper_bound=alpha_max,
                                              optimizable_frequencies_count=10)
    
    # Create anharmonicity term
    anharmonicity = qc.create_constant_pwc_operator(operator=qc.HermitianOperator(ad2a2 * chi/2),
                                                    duration=duration)
    # Create drive term
    drive = qc.get_pwc_operator_hermitian_part(qc.create_pwc_operator(
        operator=qc.NonHermitianOperator(2*a), signal=gamma))
    # Create clock shift term
    shift = qc.create_pwc_operator(operator=qc.HermitianOperator(ada/2), signal=alpha)
       
    # Create infidelity
    infidelity = qc.create_infidelity(
        hamiltonian=qc.get_pwc_operator_sum([anharmonicity, drive, shift]),
        target=qc.Target(hadamard, projector=qubit_projector))
    
    return qc.OptimizationSpecification(
        cost=infidelity, output_callable=lambda s: {
            '$\gamma$': segments_to_dicts(gamma.durations, s.run(gamma.values)),
            '$\\alpha$': segments_to_dicts(alpha.durations, s.run(alpha.values))})

# Run optimization
result = qc.optimize(create_optimization_specification)
print(f"Cost: {result['cost']}")

# Visualize optimized pulses
plot_controls(plt.figure(), result['output'])
plt.show()
Cost: 6.094458271377334e-12

Modulus (top) and angle (middle) of the optimized complex pulse $\gamma(t)$. The bottom pannel shows the amplitude of the control term $\alpha(t)$.

Wiki

Comprehensive knowledge base of quantum control theory

Explore