{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# How to perform parameter estimation with a large amount of data\n",
"**Estimate Hamiltonian model parameters using measured data and the graph-based stochastic optimization engine**\n",
"\n",
"The [stochastic optimization engine](https://docs.q-ctrl.com/references/qctrl/Functions/calculate_stochastic_optimization.html) from Boulder Opal provides a large, modular collection of configuration options that allows it to be used for a range of tasks in quantum control, including estimating Hamiltonian parameters in a model using measured data. \n",
"\n",
"In this notebook, we show how this engine can be used to determine the parameters that characterize aspects of a hardware system using a large amount of measured data.\n",
"You can read our [Characterizing your hardware using system identification in Boulder Opal](https://docs.q-ctrl.com/boulder-opal/topics/characterizing-your-hardware-using-system-identification-in-boulder-opal) topic to learn more about parameter estimation, as well as which optimization engine to use depending on your system size or the amount of measured data you have.\n",
"\n",
"\n",
"## Summary workflow\n",
"\n",
"### 1. Perform measurements to probe the system\n",
"In general the output of probe measurements will be estimates of various parameters, such as entries in a tomographic reconstruction of a state. Suitable probe pulses must be crafted which give access to terms relevant to a model being captured in order to constrain the optimization search space.\n",
"\n",
"### 2. Build a graph-based stochastic optimization encoding the problem\n",
"Represent the cost function to be minimized using a graph. Define optimization parameters in the graph using the `graph.optimization_variable` operation and sample your data using `graph.random_choices`. Execute the stochastic optimization using `qctrl.functions.calculate_stochastic_optimization` by assigning the optimization variables to output nodes of the graph.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example: Identifying parameters in a multi-qubit Hamiltonian\n",
"\n",
"In this example, we will consider how to estimate some constant parameters in a multi-qubit Hamiltonian. By preparing the system in its ground state, letting it evolve under different time-dependent pulses, and measuring, we can obtain information about the system parameters.\n",
"\n",
"We consider a $N=5$ qubit system whose Hamiltonian contains coupling, detuning, and interaction terms:\n",
"\\begin{equation}\n",
"H =\n",
"\\alpha \\, \\Omega(t) \\sum_{k=1}^N \\sigma_x^{(k)}\n",
"+\n",
"\\delta \\sum_{k=1}^N \\sigma_z^{(k)}\n",
"+\n",
"\\gamma \\bigotimes_{k=1}^N \\sigma_z^{(k)} ,\n",
"\\end{equation}\n",
"where $\\sigma_i^{(k)}$ is the $i$ Pauli matrix acting on the $k$-th qubit,\n",
"$\\alpha \\, \\Omega(t)$ is a Rabi coupling due to external pulses,\n",
"$\\delta$ is the detuning in each qubit,\n",
"and $\\gamma$ parametrizes the interaction strength between qubits.\n",
"The Rabi coupling, $\\alpha \\, \\Omega(t)$ has two components, namely unitless input values that the experimentalist is able to change ($\\Omega(t)$) and a scaling factor relating these to the physical coupling occurring in the system ($\\alpha$).\n",
"\n",
"The task of the optimizer in this case will be to determine the parameters $\\alpha$, $\\delta$, and $\\gamma$.\n",
"It can achieve this by relying on experimental data obtained by applying $\\Omega(t)$ pulses and observing the change in ground state population, which will depend on the Hamiltonian parameters to be estimated.\n",
"The sets of measured points are then provided to the optimization engine (together with the pulses used to generate them) which finds the parameters most likely to have generated that series of points.\n",
"\n",
"In this example we will rely on a large set of simulated measurement data ($\\Omega(t)$ signals and corresponding final ground state populations), which we generate first.\n",
"With the system initially in the ground state, $|0\\ldots0\\rangle$, we will send piecewise-constant controls for $\\Omega(t)$ with a fixed duration of 10 µs.\n",
"After the evolution, we measure an approximate ground state population by averaging over a large set of `shot_count` projective measurements on the final state."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"\n",
"from qctrl import Qctrl\n",
"\n",
"# Start a Boulder Opal session\n",
"qctrl = Qctrl()"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# Parameters to estimate\n",
"actual_alpha = 2 * np.pi * 0.83e6 # Hz\n",
"actual_delta = 2 * np.pi * 0.31e6 # Hz\n",
"actual_gamma = 2 * np.pi * 0.11e6 # Hz\n",
"\n",
"# System parameters\n",
"qubit_count = 5\n",
"\n",
"# Measurement parameters\n",
"shot_count = 100 # number of projective measurements shots to take\n",
"measurement_error = 1.0 / shot_count\n",
"\n",
"# Dataset parameters\n",
"dataset_size = 5000 # Size of the entire dataset\n",
"batch_size = 100 # Size for batch to use at each the stochastic optimizer iteration\n",
"\n",
"# Pulses parameters\n",
"duration = 10.0e-6 # s\n",
"segment_count = 10 # number of piecewise-constant segments"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
" 0%| | 0/100 [00:00\n",
"initial_state = graph.fock_state(2**qubit_count, 0)\n",
"\n",
"# Create signal, batch of D pwc signals\n",
"omega_signal = graph.pwc_signal(values=omega_dataset, duration=duration)\n",
"\n",
"# Create Hamiltonian operators\n",
"\n",
"# Add (Id₂ ⊗ … ⊗ Id₂) ⊗ σx ⊗ (Id₂ ⊗ … ⊗ Id₂) terms\n",
"rabi_coupling_operator = sum(\n",
" graph.pauli_kronecker_product([(\"X\", k)], qubit_count) for k in range(qubit_count)\n",
")\n",
"\n",
"# Add (Id₂ ⊗ … ⊗ Id₂) ⊗ σz ⊗ (Id₂ ⊗ … ⊗ Id₂) terms\n",
"detuning_operator = sum(\n",
" graph.pauli_kronecker_product([(\"Z\", k)], qubit_count) for k in range(qubit_count)\n",
")\n",
"\n",
"# Product of all σz operators\n",
"interaction_operator = graph.pauli_kronecker_product(\n",
" [(\"Z\", i) for i in range(qubit_count)], qubit_count\n",
")\n",
"\n",
"# Create Hamiltonian, batch of D [N,N] pwc operators\n",
"hamiltonian = (\n",
" omega_signal * actual_alpha * rabi_coupling_operator\n",
" + actual_delta * detuning_operator\n",
" + actual_gamma * interaction_operator\n",
")\n",
"\n",
"# Calculate final unitary evolution operator, shape=[D,N,N]\n",
"unitary = graph.time_evolution_operators_pwc(\n",
" hamiltonian=hamiltonian, sample_times=np.array([duration])\n",
")[:, -1]\n",
"\n",
"# Evolve intial state, shape=[D,N,1]\n",
"final_state = unitary @ initial_state[:, None]\n",
"\n",
"# Calculate final ground-state populations, shape=[D,N]\n",
"populations = graph.abs(final_state[:, :, 0]) ** 2\n",
"populations.name = \"populations\"\n",
"\n",
"# Execute the graph and extract the populations\n",
"result = qctrl.functions.calculate_graph(graph=graph, output_node_names=[\"populations\"])\n",
"calculated_populations = result.output[\"populations\"][\"value\"]\n",
"\n",
"# Take shot_count projective measurements of the system.\n",
"measured_populations = []\n",
"for pops in calculated_populations:\n",
" measurements = list(np.random.choice(2**qubit_count, size=shot_count, p=pops))\n",
" measured_populations.append(measurements.count(0) / shot_count)\n",
"\n",
"measurement_dataset = np.array(measured_populations)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Estimating the parameters of the Hamiltonian\n",
"\n",
"The parameter estimation for the Hamiltonian terms is performed by finding the form of the evolution that most closely matches the measured data points. The optimization will seek to minimize a cost function parameterized by $\\alpha$, $\\delta$ and $\\gamma$:\n",
"\n",
"$$ C = \\sum_i \\frac{[P_i-p_i(\\alpha, \\delta, \\gamma)]^2}{2(\\Delta P_i)^2} . $$\n",
"\n",
"Doing so returns the best choice of parameters that generates the original dynamics of the system, and also allows us to calculate the precision of the estimated parameters. This is done by using the Cramér–Rao bound to identify the Hessian of the cost function with the inverse of the covariance matrix for the variables estimated.\n",
"\n",
"The graph used for the optimization is very similar to the one we created to simulate the evolution and generate the measurement data, but with some key differences:\n",
"- we define optimization variables for the parameters to be estimated with `graph.optimization_variable`.\n",
"- we select a different subset of the whole data set at each optimization iteration by using `graph.random_choices`.\n",
"- we define a cost node with the negative log-likelihood (and its Hessian)."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"graph = qctrl.create_graph()\n",
"\n",
"# Define initial state |00000>\n",
"initial_state = graph.fock_state(2**qubit_count, 0)\n",
"\n",
"# Sample a batch of the omega/measurement dataset\n",
"(omega_values_batch, measurement_batch) = graph.random_choices(\n",
" data=[omega_dataset, measurement_dataset], sample_count=batch_size\n",
")\n",
"\n",
"# Parameters to be estimated\n",
"alpha = graph.optimization_variable(count=1, lower_bound=0, upper_bound=10.0e6)[0]\n",
"alpha.name = \"alpha\"\n",
"delta = graph.optimization_variable(count=1, lower_bound=0, upper_bound=5.0e6)[0]\n",
"delta.name = \"delta\"\n",
"gamma = graph.optimization_variable(count=1, lower_bound=0, upper_bound=2.0e6)[0]\n",
"gamma.name = \"gamma\"\n",
"\n",
"# Create signal, batch of D pwc signals\n",
"omega_signal = graph.pwc_signal(values=omega_values_batch, duration=duration)\n",
"\n",
"# Create Hamiltonian operators\n",
"\n",
"# Add (Id₂ ⊗ … ⊗ Id₂) ⊗ σx ⊗ (Id₂ ⊗ … ⊗ Id₂) term\n",
"rabi_coupling_operator = sum(\n",
" graph.pauli_kronecker_product([(\"X\", k)], qubit_count) for k in range(qubit_count)\n",
")\n",
"\n",
"# Add (Id₂ ⊗ … ⊗ Id₂) ⊗ σz ⊗ (Id₂ ⊗ … ⊗ Id₂) term\n",
"detuning_operator = sum(\n",
" graph.pauli_kronecker_product([(\"Z\", k)], qubit_count) for k in range(qubit_count)\n",
")\n",
"\n",
"# Product of all σz operators\n",
"interaction_operator = graph.pauli_kronecker_product(\n",
" [(\"Z\", i) for i in range(qubit_count)], qubit_count\n",
")\n",
"\n",
"# Create Hamiltonian, batch of D N×N pwc operators\n",
"hamiltonian = (\n",
" omega_signal * alpha * rabi_coupling_operator\n",
" + delta * detuning_operator\n",
" + gamma * interaction_operator\n",
")\n",
"\n",
"# Calculate final unitary evolution operator, shape=[D,N,N]\n",
"unitary = graph.time_evolution_operators_pwc(\n",
" hamiltonian=hamiltonian, sample_times=np.array([duration])\n",
")[:, -1]\n",
"\n",
"# Evolve intial state, shape=[D,N,1]\n",
"state = unitary @ initial_state[:, None]\n",
"\n",
"# Calculate final ground-state populations, shape=[D,N]\n",
"populations = graph.abs(state[:, :, 0]) ** 2\n",
"calculated_measurements = populations[:, 0]\n",
"\n",
"# Create the negative log-likelihood cost node\n",
"cost = 0.5 * graph.sum(\n",
" ((calculated_measurements - measurement_batch) / measurement_error) ** 2.0\n",
")\n",
"cost.name = \"cost\"\n",
"\n",
"# Calculate Hessian\n",
"hessian = graph.hessian(cost, [alpha, delta, gamma], name=\"hessian\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can finally execute the optimization by passing the graph to `qctrl.functions.calculate_stochastic_optimization`, along the names of the cost node and those nodes whose values we want to retrieve. In order to avoid retrieving a local minimum, we run five optimizations in parallel with the `qctrl.parallel` context manager, and then take the one yielding the best result."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
" 0%| | 0/100 [00:00