Computing the CDF with QAE in Qiskit

Question

I want to load a lognormal distribution and then use an IntegerComparator to flip a qubit ($|0\rangle$ to $|1\rangle$) if its value is less than a threshold. Then I want to use an Quantum Amplitude Estimation algorithm to calculate the probability of measuring $|1\rangle$.

My code so far is:

import matplotlib.pyplot as plt
import numpy as np

from qiskit import Aer, QuantumCircuit, QuantumRegister, execute
from qiskit.aqua.algorithms import IterativeAmplitudeEstimation
from qiskit.circuit.library import LogNormalDistribution, IntegerComparator



num_uncertainty_qubits = 3

S = 100
vol = 0.4
r = 0.04
T = 3*(30/365)


mu = np.log(S) + (r-0.5*vol**2)*T
sigma = vol*np.sqrt(T)

mean = np.exp(mu - 0.5*sigma**2)
variance = (np.exp(sigma**2)-1)*np.exp(2*mu + sigma**2)
stddev = np.sqrt(variance)

low = np.maximum(0, mean-3*stddev)
high = mean + 3*stddev

uncertainty_model = LogNormalDistribution(num_uncertainty_qubits, mu=mu, sigma=sigma**2, bounds=(low, high))

# 3 qubit LogNormalDistribution model
uncertainty_model = LogNormalDistribution(3, mu=mu, sigma=sigma, bounds=(low, high))

# function to create the quantum circuit of the IntegerComparator and the uncertainty model
# x_eval is the threshold below which the qubit should be flipped
def get_cdf_circuit(x_eval):
    qr_state = QuantumRegister(uncertainty_model.num_qubits, 'state')
    qr_obj = QuantumRegister(1, 'obj')
    qr_comp = QuantumRegister(2, 'compare')
    state_preparation = QuantumCircuit(qr_state, qr_obj, qr_comp)
    
    state_preparation.append(uncertainty_model, qr_state)
    comparator = IntegerComparator(uncertainty_model.num_qubits, x_eval, geq=False)
    
    state_preparation.append(comparator, qr_state[:]+qr_obj[:]+qr_comp[:])
    
    return state_preparation

# Function to implement the Amplitude Estimation algorithm

def run_ae_for_cdf(x_eval, epsilon=0.01, alpha=0.05, simulator='qasm_simulator'):

    state_preparation = get_cdf_circuit(x_eval)

    ae_var = IterativeAmplitudeEstimation(state_preparation=state_preparation,
                                          epsilon=epsilon, alpha=alpha,
                                          objective_qubits=[len(qr_state)]) 
    result_var = ae_var.run(quantum_instance=Aer.get_backend(simulator), shots=100)
    
    return result_var['estimation']

Broadly speaking, I want to flip the objective qubit to $|1\rangle$ if the state is less than or equal to x_eval.

On running this, I'm getting the same amplitude estimation every time, even when I use different threshold values.

I'm missing something. Please help me out.

PS: https://imgur.com/a/qAjzzEz This is the link to the QuantumCircuit. P(X) is the LogNormalDistribution and Cmp is the comparator.

PPS: Using LinearAmplitudeFunction to compare floating numbers:

def get_comparator(threshold, num_qubits, low, high):
    breakpoints = [low, threshold]
    offsets = [0,0]
    slopes = [0,0]
    f_min = 1
    f_max = 0
    objective = LinearAmplitudeFunction(
        num_qubits,
        slopes,
        offsets,
        domain=(low, high),
        image=(f_min, f_max),
        breakpoints=breakpoints
    )
    return objective

Answer 1

I think it would be also useful to see your Quantum Circuit using %matplotlib inline
and qc.draw('mpl') to see whether all gates are correctly connected. I had a similar problem with QAE and figured out by using this method that the order of qubits I tried to append the IntegerComparator onto was wrong. Cheers

Answer 2

I am also working with IterativeAmplitudeEstimation in qiskit. I was able to run QAE for 'Credit Analysis with Quantum Computing' as described in

https://qiskit.org/documentation/tutorials/finance/09_credit_risk_analysis.html

by preparing the quantum circuits from scratch. But when using the libraries as described on the page I got many errors (see below). I guess some of the libraries have simply just moved?

From which library are you loading LogNormalDistribution? I use

from qiskit.circuit.library import IntegerComparator
from qiskit.aqua.algorithms.amplitude_estimators.iqae import IterativeAmplitudeEstimation
from qiskit.aqua.components.uncertainty_models import NormalDistribution
from qiskit.aqua.components.uncertainty_models import LogNormalDistribution

I for example get the following error for LogNormalDistribution when I load the model from qiskit.aqua.components.uncertainty_models

__init__() got an unexpected keyword argument 'bounds'

or 'LogNormalDistribution' object has no attribute 'num_qubits'

I hope to get the libraries running to help you dig into your question.

Answer 3

I had to exchange len(qr_state) with num_uncertainty_qubits in the code you posted

    # in run_ae_for_cdf
    ae_var = IterativeAmplitudeEstimation(state_preparation=state_preparation,
                                          epsilon=epsilon, alpha=alpha,
                                          objective_qubits=[num_uncertainty_qubits])

but then for me it seems to work just fine.

The output is

# executing your script with python -i and then:
>>> run_ae_for_cdf(1)
0.03904972139433151
>>> run_ae_for_cdf(3)
0.35743729093195115
>>> run_ae_for_cdf(6)
0.8437041835749457
>>> run_ae_for_cdf(70)
0.9999743710958622

You have 3 qubits in your system so the largest integer you can represent is 7 ($ = 2^3 - 1$) so for values $\geq 8$ the CDF should be 1. That's because the CDF $F(x)$ is the probability of measuring a state smaller than $x$: $$ F(x) = \mathrm{Pr}(\text{measure any state } |i\rangle \text{ with } i \leq x) $$ Since $i$ can reach 7 at most, we should see $$ F(x) = \begin{cases} \text{some increasing function, if } x \leq 7 \\ 1 \text{, if } x > 7 \end{cases} $$

With your code you can generate the following plot

>>> from matplotlib import pyplot as plt
>>> i = list(range(20))
>>> f = [run_ae_for_cdf(val) for val in i]
>>> plt.axhline(1, color="k", linestyle="-")  # the max we reach
>>> plt.axvline(8, color="k", linestyle=":")  # from here on it should be 1
>>> plt.xlabel("x")
>>> plt.ylabel("CDF, $F(x)$")
>>> plt.show()

which seems to meet what we expected.