# Noise characterization and reconstruction

Reconstructing noise spectra using shaped control pulses

The Q-CTRL Python package enables you to characterize and then reconstruct power spectral densities for noise processes affecting your system. These power spectral densities can provide useful information for designing robust controls that improve the performance of the system in the presence of noise. In this notebook we show how to characterize and reconstruct noise spectra using the Q-CTRL Python package.

## Theory of SVD-based noise reconstruction

Before presenting specific examples using the Q-CTRL Python package, we give an overview of the singular value decomposition (SVD)-based algorithm underlying the reconstruction.

The general idea is to use the system to probe the spectral properties of the noise process, using the following steps:

1. generate a set of control pulses that effect identity gates but are sensitive to noises at different frequencies,
2. run the control pulses on the system and record the resulting infidelities,
3. reconstruct the spectrum based on the control pulse sensitivities (specifically their filter functions) and measured infidelities.

More precisely, consider an arbitrary qubit system. We define the noise Hamiltonian, which adds perturbatively to the noise-free Hamiltonian, as a sum of individual noise channels:

$$H_\text{noise}(t) = \sum_{i=1}^n \beta_i(t)N_i(t),$$

where $N_i$ is a noise axis operator corresponding $i$-th noise channel.

We assume that we have some array of different control pulses (resulting from step 1 above). The filter function associated with the $i$-th noise channel and the $j$-th control is defined as:

$$F_{i,j}(\omega) = \frac{1}{D}\texttt{Tr}\left(\mathcal{F}[\tilde{N}_{i,j}]\mathcal{F}[\tilde{N}_{i,j}]^\dagger\right)$$

where $D$ is the Hilbert space dimension, $\mathcal{F}$ is the Fourier transform operator and $\tilde{N}_{i,j}:=U_j^\dagger N_iU_j$ is the $i$-th noise operator in the control frame (where $U_j$ is the time evolution operator due to the $j$-th control).

Using filter function theory, the infidelity measurement collected for the $j$-th control is:

$$I_j = \sum_{i=1}^n \frac{1}{2\pi} \int_{-\infty}^{\infty} S_i(\omega)F_{i,j}(\omega) d\omega,$$

where $S_i(\omega)$ is the power spectral density for the $i$-th noise channel. This integral can be approximated as a discrete summation:

$$\tilde{I}_j = \sum_{i=1}^n \frac{1}{2\pi} \sum_p S_i[\omega_p]F_{i,j}[\omega_p] \Delta\omega.$$

Accumulating the infidelities from all the controls into a vector, concatenating the noise spectral densities associated with all noises into a vector, and flattening the discretized filter functions into a matrix with rows corresponding to controls and columns corresponding to noises and frequencies, we have:

$$\tilde{\mathbf{I}} = \frac{\Delta\omega}{2\pi}\mathbf{FS} = \mathbf{F'S},$$

where $\mathbf{F'}:=\frac{\Delta\omega}{2\pi}\mathbf{F}$. The noise spectrum $\mathbf{S}$ is thus:

$$\mathbf{S} = \mathbf{F'}^{-1}\tilde{\mathbf{I}}.$$

The singular value decomposition (SVD) gives $\mathbf{F'} = \mathbf{UDV}^\dagger$ where $\mathbf{U}$, $\mathbf{V}^\dagger$ are unitary matrices and $\mathbf{D}$ is a diagonal matrix of singular values $s_i$. Therefore, the pseudoinverse is easily calculated as:

$$\mathbf{F'}^{-1} = \mathbf{VD}^{+}\mathbf{U}^\dagger,$$

where $\mathbf{D}^+$ is a diagonal matrix of reciprocal singular values $\frac{1}{s_i}$. Combining this pseudoinverse with the previous expression for $\mathbf S$, we have:

$$\mathbf{S} = \mathbf{VD}^{+}\mathbf{U}^\dagger\tilde{\mathbf{I}}.$$

There may be several infidelity values collected for each control pulse, which gives rise to a distribution of infidelities with standard deviation $\sigma_{j}$. This translates to uncertainty in the reconstructed noise spectrum estimation. The correlation matrix for infidelity is $\Sigma_I=\texttt{diag}(\sigma_0, \ldots, \sigma_C)$. Due to linear transformation, the corresponding uncertainty in the noise spectrum estimation is given by:

$$\Sigma_S = \mathbf{VD}^{+}\mathbf{U}^\dagger\Sigma_I\mathbf{UD}^{+}\mathbf{V}^\dagger.$$

The final issue we consider is numerical instability. Small singular values in the decomposition may give rise to instabilities during the inversion process, which may necessitate truncation of singular values. Singular values may be truncated manually, or using entropy-based selection. For the latter, the entropy of the normalized singular values is calculated as:

$$E = - \sum_i \frac{s_i}{\sum_i s_i} \texttt{log}_2\left(\frac{s_i}{\sum_i s_i}\right),$$

and then only the largest $\lfloor{2^E}\rfloor$ singular values are retained.

With this in theory in hand, we now show how to perform noise reconstruction using the Q-CTRL Python package.

## Imports and initialization

All usage of the Q-CTRL Python package begins by importing the qctrl package and starting a session.

# Essential imports
import numpy as np
from qctrl import Qctrl

# Plotting imports
import matplotlib.pyplot as plt
from qctrlvisualizer import (
plot_controls,
plot_filter_functions,
)

# Convenience imports
from tqdm import tqdm

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


## Worked example: amplitude noise

In this example we consider a one-qubit system driven by a controllable Rabi rate that is subject to amplitude noise. The Hamiltonian of the system is:

\begin{align*} H(t) = &\frac{1+\beta_\Omega(t)}{2} \Omega(t) \sigma_x, \end{align*}

where $\Omega(t)$ is the controllable Rabi rate, $\beta_\Omega(t)$ is a fractional time-dependent amplitude fluctuation process and $\sigma_x$ is the Pauli X matrix. We assume that the noise process $\beta_\Omega(t)$ consists of pink noise with a small Gaussian feature, as shown below.

def gaussian(frequencies, offset, width):
return np.exp(-0.5*(frequencies - offset)**2 / width**2) / (np.sqrt(2*np.pi)*width)

def pink(frequencies, frequency_cutoff, power):
return frequency_cutoff**(power-1)/(frequencies**power + frequency_cutoff**power)

frequencies = np.linspace(0, 0.5*1e6, 1000)

amplitude_noise = 0.5e-11*(
25*pink(frequencies=frequencies, frequency_cutoff=0.1*1e6, power=6)
+ gaussian(frequencies=frequencies, offset=0.2*1e6, width=10*1e3)
)

plt.plot(frequencies/1e6, amplitude_noise*1e6)
plt.fill_between(frequencies/1e6, 0, amplitude_noise*1e6, alpha=0.25)
plt.xlabel("Frequency (MHz)")
plt.ylabel("Power density (1/MHz)")
plt.title("Amplitude noise spectrum")
plt.show()


We will demonstrate how this spectrum, which would initially be unknown in practice, can be reconstructed using the Q-CTRL Python package.

As described in the introduction to this guide, the reconstruction proceeds in three stages:

1. generate a set of control pulses that are sensitive to noises at different frequencies
2. run the control pulses on the system and record the resulting infidelities
3. reconstruct the spectrum based on the control pulse sensitivities and measured infidelities.

The Q-CTRL Python package provides functionality to perform steps 1 and 3 of this procedure.

### Creating the noise characterization

The first step is to generate the set of control pulses that will be used to probe the system. This process is known as noise characterization. To perform a noise characterization, we start by creating a noise_characterization object, which requires:

• maximum pulse duration,
• minimum time step (which describes the resolution for pulse shaping),
• maximum Rabi rate,
• pulse count.

Care must be taken when choosing these values. Let the maximum pulse duration be $T$, the minimum time step $t$, maximum Rabi rate $\Omega_\text{max}$ and pulse count $m$.

Generally speaking, the maximum noise frequency to which a characterization pulse is sensitive is $0.5/t$ (the Nyquist frequency). Therefore, to reconstruct a spectrum defined up to some frequency $f_\text{max}$, the minimum time step should be chosen as $t\approx 0.5/f_\text{max}$.

The bandwidth of each filter function, and thus the resolution of the resulting reconstruction, is roughly $2/T$. Therefore, for a maximum pulse duration $T$, you should not expect to detect features of the noise spectrum any narrower than $2/T$.

It is important that the set of filter functions provides good coverage of the whole spectrum, with no gaps. Given that the bandwidth of each filter function is $2/T$ and the total frequency range is $0.5/t$, the number of pulses $m$ should thus be chosen on the order of $0.25T/t$. A significantly smaller value will lead to gaps in the coverage of the spectrum, while a significantly larger value does not yield any improvement in coverage and is thus unnecessary.

Finally, the maximum Rabi rate $\Omega_\text{max}$ generally controls the sensitivity of the pulses to noise at the appropriate frequency—a larger maximum Rabi rate leads to more prominent filter function peaks, and therefore a higher-quality reconstruction.

Note that if the noise spectrum is to be reconstructed in a logarithmic manner (where the desired precision decreases as the frequency increases), it may be beneficial to perform the reconstruction in separate frequency intervals. For example, if a reconstruction from 0Hz to 1MHz is desired, two individual reconstructions from 0Hz–1kHz (with 100 pulses and resolution 10Hz) and 0Hz–1MHz (with 100 pulses and resolution 10kHz) could be performed. In this way, different parameters can be used for the characterization pulses for each interval, in order to achieve higher precision at low frequencies without sacrificing computational efficiency. This approach is valid whenever the maximum Rabi rate can be chosen sufficiently high that each characterization pulse in each interval is effectively immune to noise from other intervals.

# Define system parameters
omega_max = 2*np.pi * 250 * 1e6 #Hz
maximum_duration = 100 * 1e-6 #s
minimum_timestep = 1 * 1e-6 #s
pulse_count = 30

# Define noise characterization object
noise_characterization = qctrl.factories.noise_characterizations.create(
maximum_duration=maximum_duration,
minimum_timestep=minimum_timestep,
upper_bound=omega_max,
pulse_count=pulse_count,
noise_type="amplitude")


### Calculating the noise characterization

The next step is to calculate the noise characterization pulses, which we achieve via the noise_characterizations service, as shown below.

noise_characterization_result = qctrl.services.noise_characterization.calculate(noise_characterization)

100%|██████████| 100/100 [00:13<00:00,  8.30it/s]


### Visualizing the noise characterization pulses

The noise_characterization_result object contains the pulses that must be run on the system. We first show how to extract and visualize (a subset of) the pulses.

# Print some segments from the first pulse.
print(f"First 5 segments from pulse #0: {noise_characterization_result.pulses[0].segments[:5]}")

# Visualize every 5th pulse.
plot_controls(plt.figure(), {f'Pulse #{n}': noise_characterization_result.pulses[n].segments
for n in range(0, len(noise_characterization_result.pulses), 5)})
plt.show()

First 5 segments from pulse #0: [{'value': 3417.077527377352, 'duration': 1e-06}, {'value': -4402.927409657532, 'duration': 1e-06}, {'value': 5138.897382456248, 'duration': 1e-06}, {'value': -4969.63580129028, 'duration': 1e-06}, {'value': 2633.7087506445405, 'duration': 1e-06}]


### Run the characterization pulses on your system to obtain infidelities

In a noise-free system, each of the generated control pulses would effect in identity gate. In the presence of noise, however, there will be some non-zero infidelity. The next step of the reconstruction process is to calculate these infidelities. This process requires an actual system, which we do not have in this sample case, so instead we perform simulations. This computation is quite slow, so we cache values and only re-calculate if necessary (for example if the characterization pulses produced above have changed).

use_cached_data = True

if not use_cached_data:
identity = np.array([[1., 0.],[0., 1.]], dtype=np.complex)
sigma_x = np.array([[0., 1.],[1., 0.]], dtype=np.complex)

infidelity_results = []
infidelity_results_map = {}
for pulse in noise_characterization_result.pulses[len(infidelity_results):]:
system = qctrl.factories.systems.create(
name='Simulation',
hilbert_space_dimension=2,
)

drive = qctrl.factories.shift_controls.create(
name='Rabi rate',
system=system,
operator=sigma_x/2,
)

drive_pulse = qctrl.factories.custom_pulses.create(
control=drive,
segments=pulse.segments,
)

noise = qctrl.factories.control_noises.create(
name='Noise',
control=drive,
)

noise_spectral_density = qctrl.factories.piecewise_linear_noise_spectral_densities.create(
noise=noise,
sampled_points=[{'frequency': f,
'power': p,
'power_uncertainty': 0.,
'weight': 0.} for f, p in zip(frequencies, amplitude_noise)],
)

target = qctrl.factories.targets.create(
system=system,
unitary_operator=identity,
projection_operator=identity,
)

simulation = qctrl.factories.colored_noise_simulations.create(
system=system,
point_times=np.array([drive_pulse.duration]),
trajectory_count=50,
initial_state_vector=np.array([1, 0]),
)

simulation_result = qctrl.services.colored_noise_simulations.run(system)

infidelity_results.append((
simulation_result.simulations[0].average_frames[0]['infidelity'],
simulation_result.simulations[0].average_frames[0]['infidelity_uncertainty'],
))
infidelity_results_map[int(pulse.segments[0]['value'])] = infidelity_results[-1]
print(f"Infidelity: {infidelity_results[-1][0]} +/- {infidelity_results[-1][1]}")

print(f"Please copy this into the 'else' block:\ninfidelity_results_map = {infidelity_results_map}")
else:
infidelity_results_map = {3417: (1.6340149368909707e-07, 2.4336828612080006e-08), -2374: (2.7707846994573515e-07, 5.936632974647917e-08), -4442: (1.5577812511047285e-07, 3.220018220634141e-08), -7516: (2.508948990120885e-05, 3.966711417340705e-06), 5502: (3.9630760586550904e-07, 6.838101791724435e-08), 7516: (5.931400968339773e-07, 1.3311310816516383e-07), -7199: (1.8868617641065731e-06, 4.015816150621983e-07), 1259: (7.459972922929659e-07, 1.4154162795920153e-07), -6601: (9.729681222059305e-07, 1.5020957324369353e-07), 0: (2.493767435933769e-06, 3.566894898101242e-07), 7144: (3.428805460476969e-06, 5.746546406213463e-07), 6266: (4.051377373348597e-06, 6.705933662645917e-07), -1254: (4.450774070854368e-06, 8.587943086167942e-07), -8306: (1.1663049093331158e-05, 1.8310993084846337e-06), 2699: (6.35314906142459e-05, 1.4490176111734967e-05), 4692: (0.0006652327745554066, 0.00011402243942975895), 7893: (0.0002647244443594521, 4.330043989434606e-05), -3489: (0.00023321050010069344, 4.785760200442452e-05), -7591: (0.00017584300291950238, 3.4407800725077583e-05), 12351: (0.007987752648739015, 0.001552254981375516), -1315: (0.0041368414707316, 0.0010405675631251351), 2393: (0.0016632460352267132, 0.00028770910107400293), -7632: (0.0032617428696319404, 0.0005265076931762977), 7515: (0.008875602386989473, 0.0012324078641273038), -5669: (0.002797561178130705, 0.00047046067399619866), 7919: (0.0009821285622321275, 0.0001787968221699354), 7669: (0.006235578931527894, 0.00104416836270733)}
keys = [int(pulse.segments[0]['value']) for pulse in noise_characterization_result.pulses]
assert set(keys) == set(infidelity_results_map.keys()), "Set use_cached_data=False and re-run this cell"
infidelity_results = [infidelity_results_map[key] for key in keys]


### Creating the noise reconstruction

With the infidelities in hand (which, we recall, in practice would be obtained from running the characterization pulses on an actual system rather than a simulation), we may proceed with the noise reconstruction. The next step is to create a noise_reconstruction object, which encapsulates the type of singular value truncation to use (in this case we use entropy-based truncation) and will contain the result of the reconstruction.

noise_reconstruction = qctrl.factories.entropy_noise_reconstructions.create()


### Creating a system for each characterization pulse

Next we need to populate a system object corresponding to each characterization experiment. This object is what combines the measured infidelities with the applied control pulses, in the context of a noise reconstruction. The general process for creating system, control and pulse objects is described in the Setup feature guide.

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

# Store information about each characterization pulse in a dictionary
pulses_info = [{} for _ in noise_characterization_result.pulses]

# Create the system, control and pulse objects corresponding to each characterization pulse
for index, (pulse_info, pulse, infidelity_result) in enumerate(zip(tqdm(pulses_info),
noise_characterization_result.pulses,
infidelity_results)):
# Create the system, providing the noise reconstruction and measured infidelity
system = qctrl.factories.systems.create(
name=f"System {index}",
hilbert_space_dimension=2,
noise_reconstruction=noise_reconstruction,
measured_infidelity=infidelity_result[0],
measured_infidelity_uncertainty=infidelity_result[1],
)

# Create the control, with custom pulse matching the characterization pulse
control = qctrl.factories.shift_controls.create(
system=system,
name="Rabi rate",
operator=sigma_x/2,
)
control_pulse = qctrl.factories.custom_pulses.create(
control=control,
segments=pulse.segments,
)

# Save objects for later use
pulse_info['system'] = system
pulse_info['control'] = control

100%|██████████| 30/30 [03:29<00:00,  6.77s/it]


### Creating the noises

For noise reconstruction, we use a special type of "reconstructed" noise object, which specifies (in addition to the standard fields like system and control) the details to use when reconstructing the noise. In particular, we must specify minimum frequency, maximum frequency and sample count for the reconstructed spectrum, and a sample count to use for the internal filter function calculation.

# Define system parameters
minimum_frequency = 1. #Hz (must be >0)
maximum_frequency = 0.5*1e6 #Hz
noise_spectral_density_point_count = 1000
filter_function_sample_count = 200

for index, pulse_info in enumerate(tqdm(pulses_info)):
pulse_info['noise'] = qctrl.factories.reconstructed_control_noises.create(
name="Amplitude",
control=pulse_info['control'],
system=pulse_info['system'],
minimum_frequency=minimum_frequency,
maximum_frequency=maximum_frequency,
filter_function_sample_count=filter_function_sample_count,
noise_spectral_density_point_count=noise_spectral_density_point_count,
)

100%|██████████| 30/30 [01:40<00:00,  3.27s/it]


### Creating the filter functions

Noise reconstruction requires computation of filter functions. We must therefore create filter_function objects, but need not actually compute them here (that happens internally within the reconstruction computation).

for index, pulse_info in enumerate(tqdm(pulses_info)):
qctrl.factories.filter_functions.create(
noise=pulse_info['noise'],
sample_count=filter_function_sample_count,
interpolated_frequencies=np.linspace(
minimum_frequency, maximum_frequency, noise_spectral_density_point_count),
)

100%|██████████| 30/30 [01:28<00:00,  3.20s/it]


### Calculating the noise reconstruction

We have now set up, for each characterization pulse, objects that describe the pulse used for the experiment, the resulting infidelity, and the requirements for the reconstructed noise spectral density. We may now proceed to perform the noise reconstruction, using the noise_reconstructions service.

noise_reconstruction_result = qctrl.services.noise_reconstructions.calculate(noise_reconstruction)

100%|██████████| 100/100 [01:58<00:00,  6.78s/it]


### Extracting and visualizing the filter functions

Before presenting the final reconstructed noise, we note that an intermediate step of the noise reconstruction computation is calculation of the filter functions associated with each characterization experiment. Extracting and visualizing these filter functions can provide valuable insights into the mechanism underlying the reconstruction. Here, we see that each characterization experiment is highly sensitive to noise only in the vicinity of a specific frequency. This property ensures that we can produce an accurate reconstruction even with only a small number of measurements.

filter_functions = [system.controls[0].noise.filter_functions[0]
for system in noise_reconstruction_result.systems]

# Plot all filter functions.
_, ax = plt.subplots(figsize=(10,5))
for filter_function in filter_functions:
ax.plot([point['frequency']*1e-6 for point in filter_function.interpolated_points],
[point['inverse_power'] for point in filter_function.interpolated_points])
plt.xlabel("Frequency (MHz)")
plt.ylabel("Inverse power")
plt.show()


### Extracting and visualizing the reconstructed noise spectral density

Finally, we may extract the reconstructed noise spectral density out of any of the "reconstructed" noise objects. Here we use the noise in the system corresponding to the first pulse, and plot the reconstructed spectrum against the original spectrum.

sampled_points = noise_reconstruction_result.systems[0].controls[0].noise.noise_spectral_density.sampled_points

# Print some samples from the reconstructed noise spectral density.
print(f"First 5 samples: {sampled_points[:5]}")

# Plot the entire noise spectral density, including uncertainties.
_, ax = plt.subplots(figsize=(10,5))
sampled_frequencies = np.array([point['frequency'] for point in sampled_points])
sampled_power_densities = np.array([point['power'] for point in sampled_points])
sampled_power_densities_uncertainties = np.array([point['power_uncertainty'] for point in sampled_points])

ax.plot(frequencies*1e-6, amplitude_noise*1e6, label="Original")
lines = ax.plot(sampled_frequencies*1e-6, sampled_power_densities*1e6, label="Reconstructed")
ax.fill_between(
frequencies*1e-6,
(sampled_power_densities-sampled_power_densities_uncertainties)*1e6,
(sampled_power_densities+sampled_power_densities_uncertainties)*1e6,
alpha=0.35,
hatch='||',
facecolor='none',
edgecolor=lines[0].get_color(),
linewidth=0)
ax.legend()
ax.set_xlabel("Frequency (MHz)")
ax.set_ylabel("Power density (1/MHz)")
plt.show()

First 5 samples: [{'power': 1.9996463662515084e-15, 'weight': 0.8756645538929789, 'frequency': 1.0, 'power_uncertainty': 5.745432299470451e-16}
{'power': 1.9977832339211176e-15, 'weight': 0.8759532161464888, 'frequency': 501.4994994995, 'power_uncertainty': 5.726640512070351e-16}
{'power': 1.9922237556768758e-15, 'weight': 0.8768144019150228, 'frequency': 1001.998998999, 'power_uncertainty': 5.670648877748995e-16}
{'power': 1.983012964462874e-15, 'weight': 0.878240633733722, 'frequency': 1502.4984984985, 'power_uncertainty': 5.578161524949501e-16}
{'power': 1.9702254098309196e-15, 'weight': 0.8802195127098325, 'frequency': 2002.997997998, 'power_uncertainty': 5.450362300238626e-16}]


### Summary

We see that the reconstructed noise spectral density matches the original noise spectral density reasonably closely, especially in terms of the large-scale trends. By taking more precise infidelity data (with lower uncertainties) we would see an even closer match between the two spectra. We have thus demonstrated how to use the Q-CTRL Python package to reconstruct the noise spectral density of an amplitude noise process affecting a single qubit system. Similar procedures may be employed to reconstruct noise spectral densities for other noise processes (for example stochastic dephasing noises) in other types of systems.

## Wiki

Comprehensive knowledge base of quantum control theory