Getting started with the Parameter class

The QLM comes with a class dubbed Parameter that allows to describe stochastic parameters (e.g of a Hamiltonian).

Let us consider a parameter $p(t)$. Generically, we decompose it as

$$p(t)= p_0 + p_d(t) + p_s(t)$$

with

• $p_0$ a constant
• $p_d(t)$ a deterministic function of time
• $p_s(t)$ a stochastic function of time

Let us say a few more words about the stochastic part. Since it is stochastic, a given realization at a given point in time will differ from another realization at another time. Therefore, stochastic noise is only characterized by the average over many noise realizations denoted as $\langle \dots \rangle$.
In the current implementation, $p_s(t)$ is a Gaussian stochastic process with an arbitrary user-specified power spectral density (PSD) $J(f)$:

$$ J(f) = \frac{1}{t_f} \left \langle |P_s(f)|^2 \right \rangle $$

where $P_s(f)$ is the Fourier transform of $p_s(t)$ and $J(f)$ is defined in a time window $[0, t_f]$.

Below we show how one can define a Parameter and its PSD on the QLM.

import matplotlib.pyplot as plt
import numpy as np
from qat.core.variables import Variable, heaviside, cos, sin, sqrt
from qat.core import Parameter

Power spectral density - pink noise

Let us take for example a pink noise (aka $1/f$) power spectral density:

f = Variable("f", float)
cutoff_frequency = 50
intensity = 1.0
ir_frequency = 1.
psd = heaviside(f, -cutoff_frequency, -ir_frequency) * intensity * 1/sqrt(f**2 + 0.001) +\
      heaviside(f, ir_frequency, cutoff_frequency) * intensity * 1/sqrt(f**2 + 0.001) +\
      heaviside(f, -ir_frequency, ir_frequency) * intensity * 1/sqrt(ir_frequency**2 + 0.001)

The psd object is an ArithExpression, i.e. a custom QLM object for arithmetic expressions. The functions that are used inside, heaviside, sqrt are simple overloads of the usual math functions.

Create a parameter

We can now define the parameter. Here, $p_0=0$, $p_d(t)=0$, so it is purely stochastic.

sampling_frequency = 200
param = Parameter(value=0,
                  psd=psd,
                  sampling_frequency=sampling_frequency)

Let us manipulate this object:

• param.generate_noise(tf) generates a new realization of the noise between time 0 and $t_f$.
• param.get_noise() returns a tuple of arrays - for the corresponding noise amplitudes and their respective times

# Generate the noise samples
tf = 1.0
nsamples = 10
param_vals = {}
for ind in range(nsamples):
    param.generate_noise(tf=tf)
    param_vals[ind] = param.get_noise()[0]
    
# Same time points for all noise realizations of the same Parameter
param_times = param.get_noise()[1]

Plot the different noise realizations

Let us plot the 10 realizations:

plt.figure(figsize=(5, 3.5), dpi=100)
for ind in range(nsamples):
    plt.plot(param_times, param_vals[ind])
plt.title("Different noise realizations")
plt.xlabel("time $t$")
plt.ylabel("Noise amplitude (arb. units)")
plt.show()

Investigate the noise

With the help of the Fast Fourtier Transform (FFT) we can average over the noise signals to check that they all stem from the required PSD. Of course, the more the noise samples, the better the match would be.

# Take the FFTs to compute the average PSD
avg_psd = {}
for ind in range(nsamples):
    avg_psd[ind] = 1/tf * np.abs(np.fft.fft(param_vals[ind]) / sampling_frequency)**2
avg = avg_psd[0]
for ind in range(1, nsamples):
    avg += avg_psd[ind]
avg /= nsamples

# Plot the PSDs
freqs = np.linspace(0, sampling_frequency, int(tf * sampling_frequency))
plt.figure(figsize=(6, 4), dpi=110)
plt.plot(freqs, avg, label = "averaged PSD")
plt.plot(freqs, intensity/freqs, label = "target PSD")
plt.xlim(0, 52)
plt.xlabel("frequency $f$")
plt.ylabel("Noise amplitude (arb. units)")
plt.title("Comparing averaged vs target PSD")
plt.legend(loc='best')
plt.show()

/tmp/ipykernel_8653/71965848.py:14: RuntimeWarning: divide by zero encountered in divide
  plt.plot(freqs, intensity/freqs, label = "target PSD")