# How to perform Hamiltonian parameter estimation using a large amount of measured data

Estimate Hamiltonian model parameters using measured data and the graph-based stochastic optimization engine

The Q-CTRL Python package stochastic optimization engine 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.

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. You can read our 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.

## Summary workflow for Hamiltonian parameter estimation

### 1. Perform measurements to probe the system

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.

### 2. Build a graph-based stochastic optimization encoding the problem

Represent the cost function to be minimized using a graph. Define optimization parameters in the graph using qctrl.operations.optimization_variable and sample your data using qctrl.operations.random_choices. Execute the stochastic optimization using qctrl.functions.calculate_stochastic_optimization by assigning the optimization variables to output nodes of the graph.

## Worked example: Identifying parameters in a multi-qubit Hamiltonian

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.

We consider a $N=5$ qubit system whose Hamiltonian contains coupling, detuning, and interaction terms: $$H = \alpha \, \Omega(t) \sum_{k=1}^N \sigma_x^{(k)} + \delta \sum_{k=1}^N \sigma_z^{(k)} + \gamma \bigotimes_{k=1}^N \sigma_z^{(k)} \, ,$$ where $\sigma_i^{(k)}$ is the $i$ Pauli matrix acting on the $k$-th qubit, $\alpha \, \Omega(t)$ is a Rabi coupling due to external pulses, $\delta$ is the detuning in each qubit, and $\gamma$ parametrizes the interaction strength between qubits. 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$).

The task of the optimizer in this case will be to determine the parameters $\alpha$, $\delta$, and $\gamma$. 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. 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.

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. 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. 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.

import matplotlib.pyplot as plt
import numpy as np

from qctrl import Qctrl

# Define standard matrices
sigma_x = np.array([[0, 1], [1, 0]], dtype=complex)
sigma_z = np.array([[1, 0], [0, -1]], dtype=complex)

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

# Parameters to estimate
actual_alpha = 2 * np.pi * 0.83e6  # Hz
actual_delta = 2 * np.pi * 0.31e6  # Hz
actual_gamma = 2 * np.pi * 0.11e6  # Hz

# System parameters
qubit_count = 5

# Create Hamiltonian operators
rabi_coupling_operator = 0
detuning_operator = 0
interaction_operator = 1
for k in range(qubit_count):
# Kronecker products of the identity matrices before/after the n-th operator
pre_identity = np.eye(2 ** k)
post_identity = np.eye(2 ** (qubit_count - k - 1))

# Add (Id₂ ⊗ … ⊗ Id₂) ⊗ σx ⊗ (Id₂ ⊗ … ⊗ Id₂) term
rabi_coupling_operator += np.kron(np.kron(pre_identity, sigma_x), post_identity)

# Add (Id₂ ⊗ … ⊗ Id₂) ⊗ σz ⊗ (Id₂ ⊗ … ⊗ Id₂) term
detuning_operator += np.kron(np.kron(pre_identity, sigma_z), post_identity)

# Product of all σz operators
interaction_operator = np.kron(interaction_operator, sigma_z)

# Initial state |00000>, shape=[N=2^qubit_count]
initial_state = np.zeros([2 ** qubit_count])
initial_state[0] = 1.0

# Measurement parameters
shot_count = 100  # number of projective measurements shots to take
measurement_error = 1.0 / shot_count

# Dataset parameters
dataset_size = 5000  # Size of the entire dataset
batch_size = 100  # Size for batch to use at each the stochastic optimizer iteration

# Pulses parameters
duration = 10.0e-6  # s
segment_count = 10  # number of piecewise-constant segments

# Generate values for all Rabi couplings pulses in the dataset, shape [D=5000,T=10]
omega_dataset = np.random.uniform(low=0.0, high=1.0, size=[dataset_size, segment_count])

# Define graph to simulate ground-state populations
with qctrl.create_graph() as graph:

# Create signal, batch of D pwc signals
omega_signal = qctrl.operations.pwc_signal(values=omega_dataset, duration=duration)

# Create Hamiltonian, batch of D [N,N] pwc operators
hamiltonian = (
omega_signal * actual_alpha * rabi_coupling_operator
+ actual_delta * detuning_operator
+ actual_gamma * interaction_operator
)

# Calculate final unitary evolution operator, shape=[D,N,N]
unitary = qctrl.operations.time_evolution_operators_pwc(
hamiltonian=hamiltonian, sample_times=np.array([duration])
)[:, -1]

# Evolve intial state, shape=[D,N,1]
final_state = unitary @ initial_state[:, None]

# Calculate final ground-state populations, shape=[D,N]
populations = qctrl.operations.abs(final_state[:, :, 0]) ** 2
populations.name = "populations"

# Execute the graph and extract the populations
result = qctrl.functions.calculate_graph(graph=graph, output_node_names=["populations"])
calculated_populations = result.output["populations"]["value"]

# Take shot_count projective measurements of the system.
measured_populations = []
for pops in calculated_populations:

measurements = list(np.random.choice(2 ** qubit_count, size=shot_count, p=pops))
measured_populations.append(measurements.count(0) / shot_count)

measurement_dataset = np.array(measured_populations)

Your task calculate_graph (action_id="526248") has started. You can use the qctrl.get_result method to retrieve previous results.



### Estimating the parameters of the Hamiltonian

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$:

$$C = \sum_i \frac{[P_i-p_i(\alpha, \delta, \gamma)]^2}{2(\Delta P_i)^2} \, .$$

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.

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:

with qctrl.create_graph() as graph:

# Sample a batch of the omega/measurement dataset
(omega_values_batch, measurement_batch) = qctrl.operations.random_choices(
data=[omega_dataset, measurement_dataset], sample_count=batch_size
)

# Parameters to be estimated
alpha = qctrl.operations.optimization_variable(
count=1, lower_bound=0, upper_bound=10.0e6
)[0]
alpha.name = "alpha"
delta = qctrl.operations.optimization_variable(
count=1, lower_bound=0, upper_bound=5.0e6
)[0]
delta.name = "delta"
gamma = qctrl.operations.optimization_variable(
count=1, lower_bound=0, upper_bound=2.0e6
)[0]
gamma.name = "gamma"

# Create signal, batch of D pwc signals
omega_signal = qctrl.operations.pwc_signal(
values=omega_values_batch, duration=duration
)

# Create Hamiltonian, batch of D N×N pwc operators
hamiltonian = (
omega_signal * alpha * rabi_coupling_operator
+ delta * detuning_operator
+ gamma * interaction_operator
)

# Calculate final unitary evolution operator, shape=[D,N,N]
unitary = qctrl.operations.time_evolution_operators_pwc(
hamiltonian=hamiltonian, sample_times=np.array([duration])
)[:, -1]

# Evolve intial state, shape=[D,N,1]
state = unitary @ initial_state[:, None]

# Calculate final ground-state populations, shape=[D,N]
populations = qctrl.operations.abs(state[:, :, 0]) ** 2
calculated_measurements = populations[:, 0]

# Create the negative log-likelihood cost node
cost = 0.5 * qctrl.operations.sum(
((calculated_measurements - measurement_batch) / measurement_error) ** 2.0
)
cost.name = "cost"

# Calculate Hessian
hessian = qctrl.operations.hessian_matrix(
cost, [alpha, delta, gamma], name="hessian"
)


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 10 optimizations in parallel with the qctrl.parallel() context manager, and then take the one yielding the best result.

with qctrl.parallel():
parallel_results = [
qctrl.functions.calculate_stochastic_optimization(
graph=graph,
cost_node_name="cost",
output_node_names=["alpha", "delta", "gamma", "hessian"],
)
for _ in range(10)
]

result = min(parallel_results, key=lambda r: r.best_cost)

Your task calculate_stochastic_optimization (action_id="526249") has started. You can use the qctrl.get_result method to retrieve previous results.
Your task calculate_stochastic_optimization (action_id="526250") has started. You can use the qctrl.get_result method to retrieve previous results.
Your task calculate_stochastic_optimization (action_id="526251") has started. You can use the qctrl.get_result method to retrieve previous results.
Your task calculate_stochastic_optimization (action_id="526252") has started. You can use the qctrl.get_result method to retrieve previous results.
Your task calculate_stochastic_optimization (action_id="526253") has started. You can use the qctrl.get_result method to retrieve previous results.
Your task calculate_stochastic_optimization (action_id="526254") has started. You can use the qctrl.get_result method to retrieve previous results.
Your task calculate_stochastic_optimization (action_id="526255") has started. You can use the qctrl.get_result method to retrieve previous results.
Your task calculate_stochastic_optimization (action_id="526256") has started. You can use the qctrl.get_result method to retrieve previous results.
Your task calculate_stochastic_optimization (action_id="526257") has started. You can use the qctrl.get_result method to retrieve previous results.
Your task calculate_stochastic_optimization (action_id="526258") has started. You can use the qctrl.get_result method to retrieve previous results.



We can now retrieve the estimated values of the Hamiltonian parameters from the optimization results, along with the Hessian of the cost function, which allows us to calculate the uncertainties associated to the estimates.

# Calculate 2-sigma uncertainties (with 95% precision)
hessian = result.best_output["hessian"]["value"]
uncertainties = 2.0 * np.sqrt(np.diag(np.linalg.inv(hessian)))
alpha_uncertainty, delta_uncertainty, gamma_uncertainty = uncertainties

# Print parameter estimates
print(f"alpha:\t   actual =  {actual_alpha / 1e6:.4f} MHz")
print(
f"\testimated = ({result.best_output['alpha']['value'] / 1e6:.4f} "
f"± {alpha_uncertainty / 1e6:.4f}) MHz"
)
print()
print(f"delta:\t   actual =  {actual_delta / 1e6:.4f} MHz")
print(
f"\testimated = ({result.best_output['delta']['value'] / 1e6:.4f} "
f"± {delta_uncertainty / 1e6:.4f}) MHz"
)
print()
print(f"gamma:\t   actual =  {actual_gamma / 1e6:.4f} MHz")
print(
f"\testimated = ({result.best_output['gamma']['value'] / 1e6:.4f} "
f"± {gamma_uncertainty / 1e6:.4f}) MHz"
)

alpha:	   actual =  5.2150 MHz
estimated = (5.2167 ± 0.0031) MHz

delta:	   actual =  1.9478 MHz
estimated = (1.9461 ± 0.0022) MHz

gamma:	   actual =  0.6912 MHz
estimated = (0.6871 ± 0.0047) MHz


The three values we have estimated for $\alpha$, $\delta$, and $\gamma$ correspond to the actual parameters present in the Hamiltonian defined above. The errors in the parameter estimates correspond to 2 times the standard deviation and thus have 95% reliability.