How to perform model-based optimization using a Fourier basis

Create optimized pulses using CRAB techniques

Boulder Opal exposes a highly-flexible optimization engine for general-purpose gradient-based optimization. In chopped random basis (CRAB) optimization, pulses are defined via optimizable linear combinations from a set of basis functions, which can greatly reduce the dimensionality of the optimization search space. Traditionally, a randomized Fourier basis is used, although the same technique has also seen success with other bases, for example Slepian functions.

Summary workflow

1. Define basis function for signal composition in the graph

The Q-CTRL optimization engine provides a convenience function, qctrl.operations.real_fourier_pwc_signal, 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, as shown in the How to perform model-based optimization with user-defined basis functions guide.

2. Execute graph-based optimization

With the graph object created, an optimization can be run using the qctrl.functions.calculate_optimization function. The cost, the outputs, and the graph must be provided. The function returns the results of the optimization.

Worked example: CRAB optimization on a qutrit

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 + (1+\beta(t)) \left(\gamma(t) a + \gamma^*(t) a^\dagger \right) + \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, $\beta$ is a small, slowly-varying stochastic amplitude noise process, and $a = |0 \rangle \langle 1 | + \sqrt{2} |1 \rangle \langle 2 |$.

import matplotlib.pyplot as plt
import numpy as np
from qctrlvisualizer import plot_controls

from qctrl import Qctrl

# Starting a session with the API
qctrl = Qctrl()

First, define the operators and parameters for the optimization:

# Define standard matrices
transmon_levels = 3
a = np.diag(np.sqrt(np.arange(1, transmon_levels)), k=1)
ad = a.T.conj()
ada = ad @ a
ad2a2 = ad @ ada @ a
hadamard = np.array(
    [[1.0, 1.0, 0], [1.0, -1.0, 0], [0, 0, np.sqrt(2)]], dtype=complex
) / np.sqrt(2)
qubit_projector = np.pad(np.eye(2), ((0, 1), (0, 1)), mode="constant")

# 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

Create and execute the optimization graph:

# Create graph object
with qctrl.create_graph() as graph:
    # 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 = qctrl.operations.real_fourier_pwc_signal(
        duration=duration, segment_count=segment_count, randomized_frequency_count=10
    gamma_q = qctrl.operations.real_fourier_pwc_signal(
        duration=duration, segment_count=segment_count, optimizable_frequency_count=10
    gamma = qctrl.operations.complex_value(
        gamma_i * gamma_max, gamma_q * gamma_max, name="gamma"

    # Create alpha(t) signal
    alpha = qctrl.operations.real_fourier_pwc_signal(

    # Create anharmonicity term
    anharmonicity = ad2a2 * chi / 2

    # Create drive term
    drive = qctrl.operations.pwc_operator_hermitian_part(2 * gamma * a)

    # Create clock shift term
    shift = alpha * ada / 2

    # Create target operator in the qubit subspace
    target_operator =, filter_function_projector=qubit_projector

    # Create infidelity
    infidelity = qctrl.operations.infidelity_pwc(
        hamiltonian=anharmonicity + drive + shift,

# Run the optimization
optimization_result = qctrl.functions.calculate_optimization(
    cost_node_name="infidelity", output_node_names=["alpha", "gamma"], graph=graph

print(f"\nOptimized cost:\t{optimization_result.cost:.3e}")

# Plot the optimized controls
        "$\\alpha$": optimization_result.output["alpha"],
        "$\\gamma$": optimization_result.output["gamma"],
Your task calculate_optimization has started.
Your task calculate_optimization has completed.

Optimized cost:	9.952e-08