{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# How to tune the learning rate of a stochastic optimization\n",
"**Configuring the stochastic optimizer by requesting the cost history from the optimization results**\n",
"\n",
"Boulder Opal provides a highly flexible optimization engine for general-purpose optimization.\n",
"It can be directly applied to model-based control optimization and for model-based system identification.\n",
"\n",
"The optimization engine from Boulder Opal allows you to tune the parameters of your optimization.\n",
"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.\n",
"\n",
"This user guide will show how to use the cost histories returned by the optimization engine to select the learning rate of the [Adam stochastic optimizer](https://en.wikipedia.org/wiki/Stochastic_gradient_descent#Adam).\n",
"For general instructions about how to use the stochastic optimizer, see the user guide [How to optimize controls robust to strong noise sources](https://docs.q-ctrl.com/boulder-opal/user-guides/how-to-optimize-controls-robust-to-strong-noise-sources).\n",
"For an example of how to use the cost history to tune parameters of a non-stochastic optimization, see the user guide [How to tune the parameters of an optimization](https://docs.q-ctrl.com/boulder-opal/user-guides/how-to-tune-the-parameters-of-an-optimization)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Summary workflow\n",
"### 1. Define the computational graph\n",
"\n",
"The Boulder Opal optimization engine expresses all optimization problems as [data flow graphs](https://docs.q-ctrl.com/boulder-opal/topics/understanding-graphs-in-boulder-opal), which you can create with the `qctrl.create_graph` function.\n",
"The methods of the graph object allow you to represent the mathematical structure of the problem that you want to solve.\n",
"\n",
"### 2. Execute optimization and retrieve cost history\n",
"\n",
"You can calculate a stochastic optimization from an input graph using the `qctrl.functions.calculate_stochastic_optimization` function.\n",
"This function allows 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.\n",
"You can request the history of costs by passing the [`cost_history_scope` argument](https://docs.q-ctrl.com/boulder-opal/references/qctrl/Types/HistoryScope.html) to the function.\n",
"\n",
"### 3. Compare the cost histories for different choices of the learning rate\n",
"\n",
"You can retrieve the list of costs by looking at the `cost_history` attribute of the optimization result.\n",
"By plotting these results, you can then verify which choices of parameters resulted in costs that stagnated, and which ones kept improving with each iteration.\n",
"You can then keep refining the choices, and requesting more iterations for the parameter values whose costs did not stagnate.\n",
"\n",
"### 4. Restart the optimization for the most successful choice of the learning rate\n",
"\n",
"Together with the output of the optimization, the stochastic optimizer returns the `state` of the optimizer at the last iteration of the optimization.\n",
"You can use this `state` as a replacement to the `optimizer` configuration argument when calling the stochastic optimization function, in order to restart the optimization from the same point where the previous optimization had stopped."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Worked example: Tuning the learning rate and iteration count for a stochastic optimization\n",
"\n",
"Consider the problem of finding the best drive pulse to create an $X$ gate for a single-qubit system whose Hamiltonian is the following:\n",
"$$ H = (1 + \\beta) \\frac{\\Omega_x(t) \\sigma_x + \\Omega_y(t) \\sigma_y}{2} + \\frac{\\zeta}{2} \\sigma_z , $$\n",
"where $\\sigma_x$, $\\sigma_y$, and $\\sigma_z$ are the Pauli matrices, $\\Omega_x(t)$ and $\\Omega_y(t)$ are two drive pulses that you want to optimize, $\\beta$ and $\\zeta$ are a quasi-static amplitude noise and dephasing noise, respectively, both sampled from a normal distribution.\n",
"\n",
"Suppose that you want to use the stochastic optimizer to find the best pulses to reduce the effect of the noise, but you don't know which value of the learning rate will yield the best result.\n",
"In this case, you can run the optimization for different values of the learning rate, and retrieve the history of the best cost for each of them.\n",
"By plotting these best costs as a function of the iteration, you can then see which choices of the parameters stagnated at a certain level and which ones kept improving.\n",
"For those that kept improving without stagnating, you can increase the iteration count to obtain better results.\n",
"\n",
"The following example shows how to make this comparison for different values of the `learning_rate`, plotting the results of the cost histories.\n",
"At the end, we show how to continue the optimization for the learning rate that resulted in the best cost."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from qctrlvisualizer import get_qctrl_style, plot_cost_history\n",
"from qctrl import Qctrl\n",
"\n",
"plt.style.use(get_qctrl_style())\n",
"\n",
"# Start a Boulder Opal session.\n",
"qctrl = Qctrl()"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# Define optimization parameters.\n",
"omega_max = 2 * np.pi * 25e3 # Hz\n",
"duration = 200e-6 # s\n",
"segment_count = 80\n",
"noise_count = 400"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# Create optimization graph.\n",
"graph = qctrl.create_graph()\n",
"\n",
"# Define optimizable pulses, start with square pulse.\n",
"i_signal = graph.utils.real_optimizable_pwc_signal(\n",
" segment_count=segment_count,\n",
" duration=duration,\n",
" maximum=omega_max,\n",
" minimum=-omega_max,\n",
" initial_values=np.ones(segment_count) * np.pi / duration,\n",
")\n",
"q_signal = graph.utils.real_optimizable_pwc_signal(\n",
" segment_count=segment_count,\n",
" duration=duration,\n",
" maximum=omega_max,\n",
" minimum=-omega_max,\n",
" initial_values=np.zeros(segment_count),\n",
")\n",
"# Create drive term of the Hamiltonian.\n",
"drive_term = (\n",
" i_signal * graph.pauli_matrix(\"X\") / 2.0 + q_signal * graph.pauli_matrix(\"Y\") / 2.0\n",
")\n",
"\n",
"# Create amplitude noise realizations.\n",
"amplitude_noise = graph.random_normal(\n",
" mean=0, standard_deviation=0.2, shape=(noise_count, 1)\n",
")\n",
"amplitude_noise_pwc = graph.constant_pwc(\n",
" constant=amplitude_noise, duration=duration, batch_dimension_count=1\n",
")\n",
"\n",
"# Create dephasing noise realizations.\n",
"dephasing_noise = graph.random_normal(\n",
" mean=0, standard_deviation=2e4, shape=(noise_count, 1)\n",
")\n",
"dephasing_noise_pwc = graph.constant_pwc(\n",
" constant=dephasing_noise, duration=duration, batch_dimension_count=1\n",
")\n",
"dephasing_term = dephasing_noise_pwc * graph.pauli_matrix(\"Z\")\n",
"\n",
"# Define batch of infidelities.\n",
"hamiltonian = (1 + amplitude_noise_pwc) * drive_term + dephasing_term\n",
"infidelities = graph.infidelity_pwc(\n",
" hamiltonian=hamiltonian, target=graph.target(graph.pauli_matrix(\"X\"))\n",
")\n",
"\n",
"# Calculate cost as the average of the infidelities.\n",
"cost = graph.sum(infidelities) / noise_count\n",
"cost.name = \"cost\""
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
" 0%| | 0/100 [00:00"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Plot best cost history for each value of the learning rate.\n",
"plot_cost_history(\n",
" [result.cost_history.historical_best for result in results.values()],\n",
" labels=[str(learning_rate) for learning_rate in results.keys()],\n",
" y_axis_log=True,\n",
")\n",
"plt.suptitle(\"History of best average infidelities vs learning rates\")\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
" 0%| | 0/100 [00:00"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Get the complete history of best cost values.\n",
"complete_cost_history = np.minimum.accumulate(\n",
" best_result.cost_history.historical_best\n",
" + continued_results.cost_history.historical_best\n",
")\n",
"\n",
"# Plot complete graph of cost history.\n",
"plot_cost_history(list(complete_cost_history), y_axis_log=True)\n",
"plt.suptitle(\"Complete cost history for best learning rate\")\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"| Package | version |\n",
"| ----------------- | ----------- |\n",
"| Python | 3.9.12 |\n",
"| matplotlib | 3.5.1 |\n",
"| numpy | 1.21.5 |\n",
"| scipy | 1.7.3 |\n",
"| qctrl | 19.1.0 |\n",
"| qctrlcommons | 17.1.1 |\n",
"| qctrltoolkit | 1.5.0 |\n",
"| qctrlvisualizer | 3.2.1 |\n"
]
}
],
"source": [
"from qctrl.utils import print_environment_related_packages\n",
"\n",
"print_environment_related_packages()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.12"
}
},
"nbformat": 4,
"nbformat_minor": 4
}