How to tune the parameters of an optimization

Defining parameters of the optimization using the cost history and early halt conditions

Boulder Opal provides a highly flexible optimization engine for general-purpose optimization. It can be directly applied to model-based control optimization and for model-based system identification, where the model of the system is represented by an acyclic graph. You can adjust the parameters in the graph representation to customize many aspects of the optimization, such as the maximum value of the optimizable pulses, and the number of segments used in the optimization.

When tuning the parameters of an optimization, you can use the additional information provided by the history of costs obtained during the course of the optimization to see which choices of these parameters are the best for your purposes. By halting the optimization before it runs for too long and extracting the histories, you can determine which choices of parameters have values of costs that are still improving in time, and which ones stagnate at a high value.

This user guide will show how to use the cost histories returned by the optimization engine to select the number of segments for a non-stochastic optimization of a large system. For general instructions about how to perform optimization of a large system, see the user guide How to optimize controls on large sparse Hamiltonians. For an example of how to use the cost history to tune parameters of a stochastic optimization, see the user guide How to tune the learning rate of a stochastic optimization.

Summary workflow

1. Define a function to create the computational graphs

The Boulder Opal optimization engine expresses all optimization problems as data flow graphs, which you can create with the qctrl.create_graph function. The methods of the graph object allow you to represent the mathematical structure of the problem that you want to solve. You can define a function that creates a graph that represents your system, and which takes the parameter that you want to vary as an argument.

2. Execute optimization and retrieve cost history

You can calculate an optimization from an input graph using the qctrl.functions.calculate_optimization or the qctrl.functions.calculate_stochastic_optimization functions. Both functions allow you to retrieve the history of all the costs for each step of the optimization, or the history of the best value of the costs up to each iteration. You can request the history of costs by passing the argument cost_history_scope to these functions.

If you are interested in just testing an optimization that might take a while to finish, you can limit the number of iterations with the parameter max_iteration_count (for qctrl.functions.calculate_optimization) and iteration_count (for qctrl.function.calculate_stochastic_optimization). This will allow you to retrieve the beginning of the cost history, without spending too much time in the optimization.

3. Compare the cost histories for different choices of parameters

You can retrieve the list of costs by extracting the values from the cost_history attribute of the optimization result. By plotting these results, you can then verify which choices of parameters resulted in costs that stagnated, and which ones kept improving. You can then keep refining the choices, and requesting more iterations for the parameter values that did not stagnate.

Worked example: Tuning the number of segments in an optimization with six Rydberg atoms

Consider a chain of 6 atoms, each one modeled as a qubit, with $|0\rangle$ representing the ground state and $|1\rangle$ representing a Rydberg state. The total Hamiltonian of the system, is given by

$$ H = \frac{\Omega(t)}{2} \sum_{i=1}^6 \sigma_x^{(i)} - \sum_{i=1}^6 (\Delta(t) +\delta_i) n_i + \sum_{i<j}^6 \frac{V}{\vert i-j \vert^6}n_i n_j , $$

with $\sigma_x^{(i)} = \vert 0 \rangle \langle 1 \vert_i + \vert 1 \rangle \langle 0 \vert_i$, and $n_i = \vert 1\rangle\langle 1 \vert_i$. The fixed system parameters are the interaction strength between excited atoms, $V$, and the local energy shifts, $\{ \delta_i \}$. The control parameters are the effective coupling strength to the Rydberg state and the detuning, $\Omega(t)$ and $\Delta(t)$ respectively.

Suppose that you want to use the optimizer to find the best pulses to obtain a GHZ state from an initial ground state, but you don't know which choice for the number of segments will yield the best result. In this case, you can run the optimization for different values of the number of segments, and retrieve the history of the best cost for each of them. By plotting these best costs as a function of the iteration, you can then see which choices of the parameters stagnated at a certain value of the cost, and which ones kept improving. For those that kept improving without stagnating, you can increase the iteration count to obtain better results.

import matplotlib.pyplot as plt
import numpy as np
import scipy
from time import time
from qctrlvisualizer import get_qctrl_style

plt.style.use(get_qctrl_style())

from qctrl import Qctrl

# Starting a session with the API.
qctrl = Qctrl()
# Define optimization parameters.
duration = 1.1e-6  # s
omega_max = 5.0e6 * (2 * np.pi)  # Hz
delta_range = 15.0e6 * (2 * np.pi)  # Hz

# Define physical parameters of the quantum system.
interaction_strength = 24.0e6 * (2 * np.pi)  # Hz
local_shifts = np.zeros((6,))
local_shifts[[0, -1]] = -4.5e6 * (2 * np.pi)  # Hz
# Build the terms of the Hamiltonian.
h_omega_values = [0.5] * 2**6 * 6
h_omega_indices = []
h_delta_values = []
h_delta_indices = []
h_fixed_values = []
h_fixed_indices = []

for j in range(2**6):
    state = np.array([int(m) for m in np.binary_repr(j, width=6)])

    # Global shift term (big delta).
    n_exc = np.sum(state)  # Number of excited atoms.
    h_delta_values.append(-n_exc)
    h_delta_indices.append([j, j])

    # Omega term.
    for i in range(6):
        # Couple "state" and "state with qubit i flipped".
        h_omega_indices.append([j, j ^ 2 ** (5 - i)])

    # Local shift term (small deltas).
    h_fixed = -np.sum(local_shifts * state)

    # Interaction term.
    for d in range(1, 6):
        # Number of excited atoms at distance d.
        n_d = np.dot(((state[:-d] - state[d:]) == 0), state[:-d])
        h_fixed += interaction_strength * n_d / d**6

    if h_fixed != 0.0:
        h_fixed_values.append(h_fixed)
        h_fixed_indices.append([j, j])

H_omega = scipy.sparse.coo_matrix(
    (h_omega_values, np.array(h_omega_indices).T), shape=[2**6, 2**6]
)
H_delta = scipy.sparse.coo_matrix(
    (h_delta_values, np.array(h_delta_indices).T), shape=[2**6, 2**6]
)
H_fixed = scipy.sparse.coo_matrix(
    (h_fixed_values, np.array(h_fixed_indices).T), shape=[2**6, 2**6]
)
def get_recommended_k(segment_count, H_omega, H_delta, H_fixed):
    """
    Return the appropriate Krylov subspace dimension (k) for a given set of parameters, as
    required for the integration. This is obtained for the Hamiltonian with the largest spectral range.
    """
    graph = qctrl.create_graph()

    maximal_range_hamiltonian = (
        omega_max * H_omega - delta_range * H_delta + H_fixed
    ).tocoo()

    # Calculate spectral range.
    spectral_range = graph.spectral_range(operator=maximal_range_hamiltonian)

    # Calculate the appropriate Krylov subspace dimension k.
    graph.estimated_krylov_subspace_dimension_lanczos(
        spectral_range=spectral_range,
        error_tolerance=1e-6,
        duration=duration,
        maximum_segment_duration=duration / segment_count,
        name="k",
    )

    result = qctrl.functions.calculate_graph(graph=graph, output_node_names=["k"])

    return result.output["k"]["value"]


def get_optimization_graph(segment_count, krylov_subspace_dimension):
    """
    Create an optimization graph for the requested number of segments.
    """
    graph = qctrl.create_graph()

    # Define pulses.
    omega_segments = graph.pwc_signal(
        values=graph.optimization_variable(
            count=segment_count, lower_bound=0.0, upper_bound=omega_max
        ),
        duration=duration,
        name="omega",
    )
    delta_segments = graph.pwc_signal(
        values=graph.optimization_variable(
            count=segment_count, lower_bound=-delta_range, upper_bound=delta_range
        ),
        duration=duration,
        name="delta",
    )
    fixed_segment = graph.pwc_signal(
        values=np.array([1]), duration=duration, name="fixed"
    )

    shift_delta = graph.sparse_pwc_operator(signal=delta_segments, operator=H_delta)
    shift_omega = graph.sparse_pwc_operator(signal=omega_segments, operator=H_omega)
    shift_fixed = graph.sparse_pwc_operator(signal=fixed_segment, operator=H_fixed)

    # Define initial state
    initial_state = graph.fock_state(2**6, 0)

    # Create target state.
    idx_even = sum([2**n for n in range(0, 6, 2)])
    idx_odd = sum([2**n for n in range(1, 6, 2)])
    target_state = np.zeros([2**6])
    target_state[[idx_even, idx_odd]] = 1.0 / np.sqrt(2.0)

    # Evolve the initial state.
    evolved_state = graph.state_evolution_pwc(
        initial_state=initial_state,
        hamiltonian=graph.sparse_pwc_sum([shift_fixed, shift_omega, shift_delta]),
        krylov_subspace_dimension=krylov_subspace_dimension,
    )

    # Calculate state infidelity.
    infidelity = graph.state_infidelity(target_state, evolved_state, name="infidelity")

    return graph
results = {}
segment_counts = [2**n for n in range(4, 10)]
calculation_times = []

# Run optimizations.
for segment_count in segment_counts:
    krylov_subspace_dimension = get_recommended_k(
        segment_count, H_omega, H_delta, H_fixed
    )
    graph = get_optimization_graph(segment_count, krylov_subspace_dimension)
    time_reference = time()
    results[segment_count] = qctrl.functions.calculate_optimization(
        graph=graph,
        cost_node_name="infidelity",
        output_node_names=["infidelity"],
        optimization_count=1,
        max_iteration_count=300,
        cost_history_scope="HISTORICAL_BEST",
    )
    calculation_times.append(time() - time_reference)
Your task calculate_graph (action_id="936186") has completed.
Your task calculate_optimization (action_id="936188") has started. You can use the `qctrl.get_result` method to retrieve previous results.
Your task calculate_optimization (action_id="936188") has completed.
Your task calculate_graph (action_id="936225") has completed.
Your task calculate_optimization (action_id="936226") has started. You can use the `qctrl.get_result` method to retrieve previous results.
Your task calculate_optimization (action_id="936226") has completed.
Your task calculate_graph (action_id="936279") has completed.
Your task calculate_optimization (action_id="936281") has started. You can use the `qctrl.get_result` method to retrieve previous results.
Your task calculate_optimization (action_id="936281") has completed.
Your task calculate_graph (action_id="936344") has completed.
Your task calculate_optimization (action_id="936346") has started. You can use the `qctrl.get_result` method to retrieve previous results.
Your task calculate_optimization (action_id="936346") has completed.
Your task calculate_graph (action_id="936435") has completed.
Your task calculate_optimization (action_id="936436") has started. You can use the `qctrl.get_result` method to retrieve previous results.
Your task calculate_optimization (action_id="936436") has completed.
Your task calculate_graph (action_id="936440") has completed.
Your task calculate_optimization (action_id="936442") has started. You can use the `qctrl.get_result` method to retrieve previous results.
Your task calculate_optimization (action_id="936442") has completed.
# Plot history of best costs.

fig = plt.figure(figsize=(10, 5))
fig.suptitle("History of best cost values")

for segment_count, result in results.items():
    plt.semilogy(result.cost_history.historical_best[0], label=f"{segment_count}")

plt.xlabel("Iteration count")
plt.ylabel("Best cost")

plt.legend(title="Segment count", loc="upper left", bbox_to_anchor=(1.0, 1.0))
plt.show()
# Plot final costs and calculation times.

costs = [result.cost for result in results.values()]

fig, axes = plt.subplots(1, 2, figsize=(15, 5))
labels = [str(s) for s in segment_counts]
axes[0].bar(labels, costs, log=True)
axes[1].bar(labels, calculation_times)

axes[0].set_title("Best cost achieved")
axes[0].set_xlabel("Segment count")
axes[0].set_ylabel("Final cost")

axes[1].set_title("Calculation time")
axes[1].set_xlabel("Segment count")
axes[1].set_ylabel("Runtime (s)")

plt.show()