# Quantum Fourier Iterative Amplitude Estimation#

Based on the paper arXiv:2305.01686

## Problem overview#

This tutorial implements an example of the application of the Quantum Fourier Iterative Amplitude Estimation (QFIAE) method, which allows to perform Quantum Monte Carlo integration of any one-dimensional function $$f(x)$$.

The QFIAE algorithm functioning is the following:

1. The function we want to integrate $$f(x)$$, the probability distribution $$p(x)$$ and the limits of the integral $$x_{min}$$, $$x_{max}$$ are defined

2. A Quantum Neural Network (QNN) is employed to fit the function $$f(x)$$ using a linear ansatz of the following form:

where $$\mathcal{L}^{(i)}$$ is the i-th layer, that for the one-dimensional case is composed of the trainable gates:

$\mathcal{A}^{(i)}=R_z(\theta_{i,2})R_z(\theta_{i,1})R_z(\theta_{i,0}),$

and the encoding gate:

$\mathcal{S}=R_z(x).$

Now, as demonstrated in arXiv:2008.08605 the expectation value of this quantum model corresponds to a universal 1-D Fourier series:

$\langle \textit{M} \rangle (x, {\vec{\theta}}) = \sum_{w \in \Omega} c_w e^{ixw}~.$
1. The Fourier coefficients of the quantum circuit that fits the function $$f(x)$$ are extracted using discrete Fourier transform

2. The Iterative Quantum Amplitude Estimation (IQAE) model of qibo, which corresponds with the following quantum circuit:

is used to estimate the integral of each trigonometric piece, $$\sin(\omega_i x)$$ or $$\cos(\omega_j x)$$, of the Fourier series.

1. The integrals obtained are summed with their corresponding coefficients to derive the final integral result.

A summary of the workflow of the algorithm is the following:

## QFIAE Demo#

Let’s define the problem:

We are interested in computing the following example integral:

$\int_0^1 dx(1+x^2),$

so the easiest way to define our problem variables would be:

• $$x_{min}$$: $$0$$

• $$x_{max}$$: $$1$$

• $$p(x)$$: $$dx$$

• $$f(x)$$: $$1+x^2$$

## Code#

First, we must decide the variables for our problem.

Regarding the definition of the problem:

• x_min_int: the lower limit of the problem integral.

• x_max_int: the upper limit of the problem integral.

Regarding the QNN:

• nlayers: number of layers of the QNN.

• nepochs: number of epochs, i.e. iterations of the AdamDescent Optimizer.

• loss_treshold: value of the desired loss function’s treshold.

• learning_rate: learning rate of the AdamDescent Optimizer.

• nbatches: number of batches in which divide the dataset.

• ndata_points: number of data points to define $$f(x)$$.

• x_max_fit: upper limit of the function we are going to fit.

• x_min_fit: lower limit of the function we are going to fit.

Regarding the IQAE:

• nqubits: number of qubits decided to build the good state in the IQAE algorithm.

• alpha: confidence level, the target probability is 1 - alpha, has values between 0 and 1.

• epsilon: target precision for estimation target $$a$$, has values between 0 and 0.5.

• nshots: number of shots used to perform IQAE.

• method_iqae: statistical method employed by IQAE. It will be chernoff or beta.

Regardind the results:

• qfiae_integral: estimated value of the integral using QFIAE.

• iae_error: statistical uncertainty in the integral estimation due to IQAE.

Include necessary packages:

[1]:

import numpy as np

from qibo import Circuit, gates
from qibo.models.iqae import IQAE


### 1. Defining $$f(x)$$, $$x_{min}$$, $$x_{max}$$#

[2]:

def f(x):
"""Target function f(x) = 1 + x ^ 2"""

y = 1 + x ** 2
return y

x_min_int = 0
x_max_int = 1


### 2. Defining and training the QNN#

We import the QNN: vqregressor_linear_ansatz, which is a modified version of the vqregressor found at https://github.com/qiboteam/qibo/tree/master/examples/vqregressor. In our modified version, we replace the ansatz function used in the original implementation with the one we have previously defined.

[3]:

from vqregressor_linear_ansatz import VQRegressor_linear_ansatz


We define the parameters of the QNN

[4]:

# Note that the chosen number of layers is intentionally larger than necessary to fit the function.
# This configuration ensures high accuracy while keeping a reasonable number of epochs.
nlayers = 10
nepochs = 100
loss_treshold = 10^-4
learning_rate = 0.05
nbatches = 1
ndata_points = 200
# Note they can be different than x_min,x_max, but [x_max,x_min] ⊆ [x_max_fit,x_min_fit]
x_max_fit = 1
x_min_fit = -1

[5]:

np.random.seed(27)
vqr = VQRegressor_linear_ansatz(layers=nlayers,
ndata=ndata_points,
function=f, xmax=x_max_fit,
xmin=x_min_fit)
# With the selected configuration, the training usually takes around 30'
vqr.train_with_psr(
epochs=nepochs,
learning_rate=learning_rate,
batches=nbatches,
J_treshold=loss_treshold,
)
vqr.show_predictions("Predictions of the QNN after training", False)

[Qibo 0.1.16|INFO|2023-07-31 08:50:42]: Using qibojit (numba) backend on /CPU:0

Iteration  1  epoch  1  | loss:  0.7513285431118756
Iteration  2  epoch  2  | loss:  0.49078195900193
Iteration  3  epoch  3  | loss:  0.25449817735514246
Iteration  4  epoch  4  | loss:  0.1246609331106105
Iteration  5  epoch  5  | loss:  0.0656780490580047
Iteration  6  epoch  6  | loss:  0.033157170041590786
Iteration  7  epoch  7  | loss:  0.021080634318318345
Iteration  8  epoch  8  | loss:  0.02302933284076279
Iteration  9  epoch  9  | loss:  0.02804566383739096
Iteration  10  epoch  10  | loss:  0.02993298809989462
Iteration  11  epoch  11  | loss:  0.02809272996633005
Iteration  12  epoch  12  | loss:  0.024260018563301875
Iteration  13  epoch  13  | loss:  0.020466737677239978
Iteration  14  epoch  14  | loss:  0.017961154451515637
Iteration  15  epoch  15  | loss:  0.01635525503721851
Iteration  16  epoch  16  | loss:  0.013824870329752614
Iteration  17  epoch  17  | loss:  0.009480363571353416
Iteration  18  epoch  18  | loss:  0.005428317914456559
Iteration  19  epoch  19  | loss:  0.004747379983842302
Iteration  20  epoch  20  | loss:  0.007261736229780217
Iteration  21  epoch  21  | loss:  0.008972805823541971
Iteration  22  epoch  22  | loss:  0.006995705828705413
Iteration  23  epoch  23  | loss:  0.0031234835775331904
Iteration  24  epoch  24  | loss:  0.0010645149525778657
Iteration  25  epoch  25  | loss:  0.001660608449799947
Iteration  26  epoch  26  | loss:  0.0024468184819136764
Iteration  27  epoch  27  | loss:  0.0019337945857293712
Iteration  28  epoch  28  | loss:  0.0012357754610404163
Iteration  29  epoch  29  | loss:  0.0013015455646036862
Iteration  30  epoch  30  | loss:  0.0016682446926897167
Iteration  31  epoch  31  | loss:  0.0016495496176469587
Iteration  32  epoch  32  | loss:  0.0012520977035885572
Iteration  33  epoch  33  | loss:  0.0009722650541584987
Iteration  34  epoch  34  | loss:  0.0010738639185073821
Iteration  35  epoch  35  | loss:  0.001237623835436679
Iteration  36  epoch  36  | loss:  0.001079914927835038
Iteration  37  epoch  37  | loss:  0.0008087925855009583
Iteration  38  epoch  38  | loss:  0.0008412140573143836
Iteration  39  epoch  39  | loss:  0.0010243360494146807
Iteration  40  epoch  40  | loss:  0.0009260517614442422
Iteration  41  epoch  41  | loss:  0.0006083674444202678
Iteration  42  epoch  42  | loss:  0.00048552087540479593
Iteration  43  epoch  43  | loss:  0.0006110987021632932
Iteration  44  epoch  44  | loss:  0.0006559683517174265
Iteration  45  epoch  45  | loss:  0.0005150053878726019
Iteration  46  epoch  46  | loss:  0.0004114721142009153
Iteration  47  epoch  47  | loss:  0.0004471177670137205
Iteration  48  epoch  48  | loss:  0.0004770075457458069
Iteration  49  epoch  49  | loss:  0.00040857499834427296
Iteration  50  epoch  50  | loss:  0.00033020904524591827
Iteration  51  epoch  51  | loss:  0.0003179087176282798
Iteration  52  epoch  52  | loss:  0.00032425348191094495
Iteration  53  epoch  53  | loss:  0.0003081269055732997
Iteration  54  epoch  54  | loss:  0.00029739364175727137
Iteration  55  epoch  55  | loss:  0.0002972067254557964
Iteration  56  epoch  56  | loss:  0.0002790839752707141
Iteration  57  epoch  57  | loss:  0.000250994009314296
Iteration  58  epoch  58  | loss:  0.00024068099758678817
Iteration  59  epoch  59  | loss:  0.00024012049242896222
Iteration  60  epoch  60  | loss:  0.00023072309010057932
Iteration  61  epoch  61  | loss:  0.0002204255806400206
Iteration  62  epoch  62  | loss:  0.00021693211012467254
Iteration  63  epoch  63  | loss:  0.00021032053408251615
Iteration  64  epoch  64  | loss:  0.00019814302410227347
Iteration  65  epoch  65  | loss:  0.00018970104609909337
Iteration  66  epoch  66  | loss:  0.0001849409887409433
Iteration  67  epoch  67  | loss:  0.000178175085771177
Iteration  68  epoch  68  | loss:  0.00017235266528996997
Iteration  69  epoch  69  | loss:  0.00016920670127216307
Iteration  70  epoch  70  | loss:  0.00016352665275946328
Iteration  71  epoch  71  | loss:  0.0001558330982688153
Iteration  72  epoch  72  | loss:  0.00015122871287318634
Iteration  73  epoch  73  | loss:  0.0001481650541488665
Iteration  74  epoch  74  | loss:  0.0001435783238101642
Iteration  75  epoch  75  | loss:  0.000139808928386417
Iteration  76  epoch  76  | loss:  0.00013763074842478883
Iteration  77  epoch  77  | loss:  0.00013405673992693997
Iteration  78  epoch  78  | loss:  0.00012925759824435786
Iteration  79  epoch  79  | loss:  0.00012568131156775818
Iteration  80  epoch  80  | loss:  0.000122960138626043
Iteration  81  epoch  81  | loss:  0.00012002635696736961
Iteration  82  epoch  82  | loss:  0.00011738809629738198
Iteration  83  epoch  83  | loss:  0.00011504380983272437
Iteration  84  epoch  84  | loss:  0.00011252787350936655
Iteration  85  epoch  85  | loss:  0.00011014441030141236
Iteration  86  epoch  86  | loss:  0.00010794102346639878
Iteration  87  epoch  87  | loss:  0.00010557841843900766
Iteration  88  epoch  88  | loss:  0.00010325318195926471
Iteration  89  epoch  89  | loss:  0.00010112572290540777
Iteration  90  epoch  90  | loss:  9.901680049760298e-05
Iteration  91  epoch  91  | loss:  9.702508371139033e-05
Iteration  92  epoch  92  | loss:  9.525568074239291e-05
Iteration  93  epoch  93  | loss:  9.354329673964097e-05
Iteration  94  epoch  94  | loss:  9.188064428852769e-05
Iteration  95  epoch  95  | loss:  9.032782168927708e-05
Iteration  96  epoch  96  | loss:  8.881823164035069e-05
Iteration  97  epoch  97  | loss:  8.739132343365943e-05
Iteration  98  epoch  98  | loss:  8.607781395924396e-05
Iteration  99  epoch  99  | loss:  8.477887801842459e-05
Iteration  100  epoch  100  | loss:  8.350302085692635e-05


### 3. Obtaining the Fourier series from the QNN#

[6]:

# We import the qibo function that obtains the Fourier series out of a
# quantum circuit that returns an expectation value
from qibo.models.utils import fourier_coefficients
# We get the optimized quantum circuit and the normalized Fourier coefficients
quantum_circuit_optimized = vqr.one_prediction
normalized_coeffs = fourier_coefficients(quantum_circuit_optimized,
n_inputs=1, degree=nlayers, lowpass_filter=True)
# We obtain the non-normalized coeffs
coeffs = normalized_coeffs * vqr.norm


Let’s see how they look like

[7]:

print(coeffs)

[ 9.91228857e-01+0.j         -9.14152432e-04-0.04336848j
3.40828521e-01-0.02305663j -1.77033953e-01+0.00110917j
-3.63049787e-01+0.12160816j  7.33018369e-02-0.05831041j
2.28576955e-01-0.1387953j  -2.13031210e-03+0.16034618j
-1.57394048e-01-0.03133756j  5.23539080e-02-0.02270211j
5.30676569e-03-0.00145169j  5.30676569e-03+0.00145169j
5.23539080e-02+0.02270211j -1.57394048e-01+0.03133756j
-2.13031210e-03-0.16034618j  2.28576955e-01+0.1387953j
7.33018369e-02+0.05831041j -3.63049787e-01-0.12160816j
-1.77033953e-01-0.00110917j  3.40828521e-01+0.02305663j
-9.14152432e-04+0.04336848j]


Now we check that Fourier series obtained it’s correct and also fits the function

[8]:

# We import some functions to obtain the Fourier series from the coefficients and to plot these results
from utils_qfiae import fourier_series, plot_data_fourier
# We plot the original function and the Fourier approximation
y_fourier = fourier_series(coeffs, vqr.features)
plot_fourier, accuracy_qnn = plot_data_fourier(vqr.labels * vqr.norm, vqr.features, y_fourier)


### 4. Integrating every trigonometric piece#

First we should rewrite the Fourier coefficients in a way that will be easy to integrate the trigonometric functions Sine and Cosine

[9]:

# We import some functions to convert exponential-like Fourier coefficients into trigonometric-like and to
# convert this trigonometric coeffs into an array that provides insight about if we get a Sine or a Cosine
from utils_qfiae import exp_fourier_to_trig, trig_to_final_array
trig_coeffs = exp_fourier_to_trig(coeffs)
fourier_array = trig_to_final_array(trig_coeffs)


Let’s visualize what we have and note that in the first columnn 0 stands for Cosine and 1 stands for Sine, in the second column we have the frequencies and in the third one the coefficients

[10]:

print(fourier_array)

[[ 0.00000000e+00  0.00000000e+00  9.91228857e-01]
[ 0.00000000e+00  1.00000000e+00 -1.82830486e-03]
[ 0.00000000e+00  2.00000000e+00  6.81657041e-01]
[ 0.00000000e+00  3.00000000e+00 -3.54067907e-01]
[ 0.00000000e+00  4.00000000e+00 -7.26099574e-01]
[ 0.00000000e+00  5.00000000e+00  1.46603674e-01]
[ 0.00000000e+00  6.00000000e+00  4.57153911e-01]
[ 0.00000000e+00  7.00000000e+00 -4.26062419e-03]
[ 0.00000000e+00  8.00000000e+00 -3.14788097e-01]
[ 0.00000000e+00  9.00000000e+00  1.04707816e-01]
[ 0.00000000e+00  1.00000000e+01  1.06135314e-02]
[ 1.00000000e+00  1.00000000e+00  8.67369616e-02]
[ 1.00000000e+00  2.00000000e+00  4.61132575e-02]
[ 1.00000000e+00  3.00000000e+00 -2.21834947e-03]
[ 1.00000000e+00  4.00000000e+00 -2.43216317e-01]
[ 1.00000000e+00  5.00000000e+00  1.16620823e-01]
[ 1.00000000e+00  6.00000000e+00  2.77590597e-01]
[ 1.00000000e+00  7.00000000e+00 -3.20692363e-01]
[ 1.00000000e+00  8.00000000e+00  6.26751169e-02]
[ 1.00000000e+00  9.00000000e+00  4.54042108e-02]
[ 1.00000000e+00  1.00000000e+01  2.90337531e-03]]


Before implementeing the IQAE model, we should define the quantum circuits $$\mathcal{A}$$ and $$\mathcal{Q}$$.

The circuit $$\mathcal{A}$$ is defined as:

$|\psi\rangle=\mathcal{A}|0\rangle_{n+1}=\sqrt{a}|\tilde{\psi}_1\rangle|1\rangle +\sqrt{1-a}|\tilde{\psi}_0\rangle|0\rangle,$

where $$a\in[0,1]$$ is the parameter to be estimated, while $$|\tilde{\psi}_1\rangle$$ and $$|\tilde{\psi}_0\rangle$$ are the $$n$$-qubit good and bad states, respectively.

From a more practical point of view, this circuit can be decomposed in two pieces:

$\mathcal{A}=\left(\mathcal{P}\otimes \mathcal{I}_1\right)\mathcal{R},$

where $$\mathcal{I}_1$$ is the identity matrix acting over the last qubit.

So, to get $$\mathcal{A}$$ we have to define the operator $$\mathcal{P}$$, which encodes the probability distribution $$p(x)$$, and the operator $$\mathcal{R}$$, which encodes the function we are going to integrate, $$\sin(x)^2$$, into the interval $$[b_{min},b_{max}]$$

[11]:

def operator_p(circ, list_of_qubits):
"""Generating the uniform probability distribution p(x)=1/2^n"""

for qubit in list_of_qubits:

def operator_r(circ, list_of_qubits, qubit_to_measure, b_max, b_min):
"""Encoding the function Sin(x)^2 to be integrated in the interval [b_min,b_max]"""

number_of_qubits = qubit_to_measure
gates.RY(q=qubit_to_measure, theta=(b_max - b_min) / 2 ** number_of_qubits + 2 * b_min)
)
for qubit in list_of_qubits:
gates.CRY(qubit, qubit_to_measure,
(b_max - b_min) / 2 ** (number_of_qubits-qubit-1))
)


Then, the circuit $$\mathcal{A}$$ will be:

[12]:

def create_circuit_a(n_qubits, b_max, b_min):
"""Creation of circuit A"""

circuit_a = Circuit(n_qubits + 1)
list_of_qubits = list(range(n_qubits))
operator_p(circuit_a, list_of_qubits)
operator_r(circuit_a, list_of_qubits, n_qubits, b_max, b_min)
return circuit_a


Let’s see how the circuit A looks like for $$n_{qubits}$$=4

[13]:

circuit_a_prob = create_circuit_a(n_qubits = 4,b_max = x_max_int, b_min = x_min_int)
print(circuit_a_prob.draw())

q0: ─H──o───────────
q1: ─H──|──o────────
q2: ─H──|──|──o─────
q3: ─H──|──|──|──o──
q4: ─RY─RY─RY─RY─RY─


On the other hand, the circuit $$\mathcal{Q}$$ is defined as:

$\mathcal{Q}=-\mathcal{A} S_0 \mathcal{A}^{-1} S_\chi,$

where the operator $$S_0$$ tags with a negative sign the $$|0\rangle_{n+1}$$ state and does nothing to the other states, and the operator $$S_\chi$$ changes the sign to the good state, i.e. $$S_\chi |\tilde{\psi}_1\rangle|1\rangle=-|\tilde{\psi}_1\rangle|1\rangle$$.

So, let’s define $$S_0$$ and $$S_\chi$$ to get $$\mathcal{Q}$$:

[14]:

def operator_s_0(circ, qubit_to_measure):
"""Tagging with a negative sign any state different to |0>_{n+1}. Note that this is
the opposite of what S_0 is supposed to do. The explanation is that the operator

def operator_s_chi(circ, qubit_to_measure):
"""Computing reflection operator (I - 2|0><0|)"""

number_of_qubits = qubit_to_measure
list_of_qubits = list(range(number_of_qubits))
for qubit in list_of_qubits:
for qubit in list_of_qubits:


Then the circuit $$\mathcal{Q}$$ will be:

[15]:

def create_circuit_q(n_qubits, circuit_a):
"""Creation of circuit Q"""

circuit_q = Circuit(n_qubits + 1)
operator_s_0(circuit_q, n_qubits)
*range(0, n_qubits + 1, 1)))
operator_s_chi(circuit_q, n_qubits)
*range(0, n_qubits + 1, 1)))
return circuit_q


Let’s see how the circuit Q looks like for $$n_{qubits}$$=4

[16]:

print(create_circuit_q(n_qubits=4,circuit_a=circuit_a_prob).draw())

q0: ────────────o──H──X─o─X─H──o───────────
q1: ─────────o──|──H──X─o─X─H──|──o────────
q2: ──────o──|──|──H──X─o─X─H──|──|──o─────
q3: ───o──|──|──|──H──X─o─X─H──|──|──|──o──
q4: ─Z─RY─RY─RY─RY─RY─X─Z─X─RY─RY─RY─RY─RY─


Since we already have the circuits $$\mathcal{A}$$ and $$\mathcal{Q}$$, we can proceed to execute the IQAE model.

But first let’s do the last transformation of the Fourier coefficients

[17]:

from utils_qfiae import coeffs_array_to_class
coeffs_class = coeffs_array_to_class(fourier_array)


Then, we will create the following function that will integrate each sine and cosine and then sum every term to obtain the final Monte Carlo integral.

To create this integrator we will use the following mathematical properties:

$$\frac{1}{x_{min}-x_{min}}\int_{x_{min}}^{x_{max}}dx\sin(x \omega/2)^2=\frac{1}{b_{max}-b_{min}} \int_{b_{min}}^{b_{max}}dx\sin(x)^2, \tag{1}$$

if $$b_{max}=x_{max}\omega / 2$$ and $$b_{min}=x_{min} \omega / 2$$.

$$\sin(x) = \cos (x - \pi/2) \tag{2}$$

$$\cos(x) = 1 - 2 \sin(x/2) ^ 2 \tag{3}$$

[18]:

def qfiae_integrator(coeffs_f, x_max, x_min, nbit, shots, alpha, epsilon, method):
"""Takes the Fourier coefficients of the desired function and calculates the integral of the function
in the specified interval."""

mc_integral = 0
error_iqae = 0
if coeffs_f[0].angle == 0:
mc_integral += coeffs_f[0].coeff * (x_max - x_min)
for n in range(1,len(coeffs_f)):
angle = coeffs_f[n].angle
if coeffs_f[n].function == 0:  # cosine
# Here we use Eq.(1)
a_circuit = create_circuit_a(n_qubits=nbit,
b_max=x_max * angle / 2,
b_min=x_min * angle / 2)

if coeffs_f[n].function == 1:  # sine
# Here we use Eq.(2)
a_circuit = create_circuit_a(n_qubits=nbit,
b_max=x_max * angle / 2 - np.pi / 4,
b_min=x_min * angle / 2 - np.pi / 4)

q_circuit = create_circuit_q(n_qubits=nbit, circuit_a=a_circuit)
result = IQAE(a_circuit, q_circuit, alpha, epsilon, n_shots=shots, method=method).execute()
a_estimated = result.estimation
# This a_estimated corresponds to the integral of the right part of Eq.(1)
error = result.epsilon_estimated
# Here we use Eq.(3) and multiply the result by (x_max - x_min) to get rid of the normalization
mc_integral += (1 - 2 * a_estimated) * coeffs_f[n].coeff * (x_max - x_min)
error_iqae += (2 * abs(coeffs_f[n].coeff) * error * (x_max - x_min))**2

return mc_integral , np.sqrt(error_iqae)


Finally, we define the parameters for the IQAE model and execute the qfiae_integrator function

[19]:

nqubits = 4
nshots = 1024
alpha = 0.05
epsilon = 0.005
method_iqae = "chernoff"
qfiae_integral, iae_error = qfiae_integrator(coeffs_class, x_max=x_max_int, x_min=x_min_int,
nbit=nqubits, shots=nshots, alpha=alpha,
epsilon=epsilon, method=method_iqae)


Once the QFIAE algorithm has been performed we can visualize and compare the results with the ones obtained by classical calculation

[20]:

from scipy.integrate import quad
classical_integral, classical_error = quad(f, x_min_int, x_max_int)
print("The result obtained with QFIAE is", qfiae_integral,
" ± ", iae_error)
print("The result obtained with classical numerical integration is",
classical_integral, " ± ", classical_error)

The result obtained with QFIAE is 1.3368772283309165  ±  0.00781780414187758
The result obtained with classical numerical integration is 1.3333333333333333  ±  1.4802973661668752e-14