Note

Cette page a été générée à partir de tutorials/finance/9_credit_risk_analysis.ipynb.

Exécuter en mode interactif dans le IBM Quantum lab.

Analyse du Risque de Crédit

Introduction

This tutorial shows how quantum algorithms can be used for credit risk analysis. More precisely, how Quantum Amplitude Estimation (QAE) can be used to estimate risk measures with a quadratic speed-up over classical Monte Carlo simulation. The tutorial is based on the following papers:

• Quantum Risk Analysis. Stefan Woerner, Daniel J. Egger. [Woerner2019]

• Credit Risk Analysis using Quantum Computers. Egger et al. (2019) [Egger2019]

A general introduction to QAE can be found in the following paper:

• Quantum Amplitude Amplification and Estimation. Gilles Brassard et al.

The structure of the tutorial is as follows:

1. Problem Definition

2. Uncertainty Model

3. Expected Loss

4. Cumulative Distribution Function

5. Value at Risk

6. Conditional Value at Risk

import numpy as np
import matplotlib.pyplot as plt

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

Définition du Problème

Dans ce tutoriel, nous voulons analyser le risque de crédit d’un portefeuille de \(K\) actifs. La probabilité par défaut de chaque actif \(k\) suit un modèle Gaussien d’indépendance conditionnelle, c’est-à-dire étant donné une valeur \(z\) échantillonnée à partir d’une variable aléatoire latente \(Z\) suivant une distribution normale standard, la probabilité par défaut de l’actif \(k\) est donnée par

\[p_k(z) = F\left( \frac{F^{-1}(p_k^0) - \sqrt{\rho_k}z}{\sqrt{1 - \rho_k}} \right)\]

où \(F\) indique la fonction de distribution cumulée de \(Z\), \(p_k^0\) est la probabilité par défaut de l’actif \(k\) pour \(z=0\) et \(\rho_k\) est la sensibilité de la probabilité par défaut de l’actif \(k\) par rapport à \(Z\). Ainsi, compte tenu d’une réalisation concrète de \(Z\) les événements par défaut individuels sont supposés être indépendants les uns des autres.

Nous nous intéressons à l’analyse des mesures de risque de perte totale

\[L = \sum_{k=1}^K \lambda_k X_k(Z)\]

où \(\lambda_k\) indique la perte par défaut de l’actif \(k\), et a donné \(Z\), \(X_k(Z)\) désigne une variable de Bernoulli représentant l’événement par défaut de l’actif \(k\). Plus précisément, nous nous interessons à la valeur espérée \(\mathbb{E}[L]\), la valeur à risque (VaR) de \(L\) et la valeur conditionnelle à risque de \(L\) (également appelée déficit attendu). Là où VaR et CVaR sont définis comme

\[\text{VaR}_{\alpha}(L) = \inf \{ x \mid \mathbb{P}[L <= x] \geq 1 - \alpha \}\]

avec le niveau de confiance \(\alpha \in [0, 1]\), et

\[\text{CVaR}_{\alpha}(L) = \mathbb{E}[ L \mid L \geq \text{VaR}_{\alpha}(L) ].\]

Pour plus de détails sur le modèle considéré, voir par exemple Regulatory Capital Modeling for Credit Risk. Marek Rutkowski, Silvio Tarca

Le problème est défini par les paramètres suivants : - nombre de qubits utilisés pour représenter \(Z\), notés \(n_z\) - valeur de troncature pour \(Z\), dénoté par \(z_{\text{max}}\), i. ., Z est supposé prendre \(2^{n_z}\) valeurs équidistantes dans \(\{-z_{max}, ... +z_{max}\}\) - les probabilités par défaut de base pour chaque ressource \(p_0^k \in (0, 1)\), \(k=1, ... K\) - Sensibilité des probabilités par défaut par rapport à \(Z\), notés \(\rho_k \in [0, 1)\) - perte par défaut pour la ressource \(k\), notée \(\lambda_k\) - niveau de confiance pour VaR / CVaR \(\alpha \in [0, 1]\).

# set problem parameters
n_z = 2
z_max = 2
z_values = np.linspace(-z_max, z_max, 2**n_z)
p_zeros = [0.15, 0.25]
rhos = [0.1, 0.05]
lgd = [1, 2]
K = len(p_zeros)
alpha = 0.05

Modèle d’incertitude

Nous construisons maintenant un circuit qui permet de charger le modèle d’incertitude. Ceci peut être réalisé en créant un état quantique dans un registre de \(n_z\) qubits qui représente \(Z\) suivant une distribution normale standard. Cet état est ensuite utilisé pour contrôler des rotations en Y sur un second registre de qubits de \(K\) de qubit \(k\) représente l’événement par défaut de l’asset \(k\). L’état quantique résultant peut être écrit comme

\[ |\Psi\rangle = \sum_{i=0}^{2^{n_z}-1} \sqrt{p_z^i} |z_i \rangle \bigotimes_{k=1}^K \left( \sqrt{1 - p_k(z_i)}|0\rangle + \sqrt{p_k(z_i)}|1\rangle\right),\]

où nous notons \(z_i\) la \(i\)-ème valeur de la valeur discrétisée et tronquée \(Z\) [Egger2019].

from qiskit.finance.applications import GaussianConditionalIndependenceModel as GCI
u = GCI(n_z, z_max, p_zeros, rhos)

u.draw()

     ┌───────┐┌─────────┐┌─────────┐
q_0: ┤0      ├┤0        ├┤0        ├
     │  P(X) ││         ││         │
q_1: ┤1      ├┤1 LinRot ├┤1        ├
     └───────┘│         ││  LinRot │
q_2: ─────────┤2        ├┤         ├
              └─────────┘│         │
q_3: ────────────────────┤2        ├
                         └─────────┘

Nous utilisons maintenant le simulateur pour valider le circuit qui construit l’état \(|\Psi\rangle\) et calculer les valeurs exactes correspondantes pour - la perte attendue \(\mathbb{E}[L]\) - PDF et CDF de \(L\) - la valeur à risque \(VaR(L)\) et la probabilité correspondante - la valeur conditionnelle à risque \(CVaR(L)\)

# run the circuit and analyze the results
job = execute(u, backend=Aer.get_backend('statevector_simulator'))

# analyze uncertainty circuit and determine exact solutions
p_z = np.zeros(2**n_z)
p_default = np.zeros(K)
values = []
probabilities = []
num_qubits = u.num_qubits
for i, a in enumerate(job.result().get_statevector()):

    # get binary representation
    b = ('{0:0%sb}' % num_qubits).format(i)
    prob = np.abs(a)**2

    # extract value of Z and corresponding probability
    i_normal = int(b[-n_z:], 2)
    p_z[i_normal] += prob

    # determine overall default probability for k
    loss = 0
    for k in range(K):
        if b[K - k - 1] == '1':
            p_default[k] += prob
            loss += lgd[k]
    values += [loss]
    probabilities += [prob]

values = np.array(values)
probabilities = np.array(probabilities)

expected_loss = np.dot(values, probabilities)

losses = np.sort(np.unique(values))
pdf = np.zeros(len(losses))
for i, v in enumerate(losses):
    pdf[i] += sum(probabilities[values == v])
cdf = np.cumsum(pdf)

i_var = np.argmax(cdf >= 1-alpha)
exact_var = losses[i_var]
exact_cvar = np.dot(pdf[(i_var+1):], losses[(i_var+1):])/sum(pdf[(i_var+1):])

print('Expected Loss E[L]:                %.4f' % expected_loss)
print('Value at Risk VaR[L]:              %.4f' % exact_var)
print('P[L <= VaR[L]]:                    %.4f' % cdf[exact_var])
print('Conditional Value at Risk CVaR[L]: %.4f' % exact_cvar)

Expected Loss E[L]:                0.6409
Value at Risk VaR[L]:              2.0000
P[L <= VaR[L]]:                    0.9591
Conditional Value at Risk CVaR[L]: 3.0000

# plot loss PDF, expected loss, var, and cvar
plt.bar(losses, pdf)
plt.axvline(expected_loss, color='green', linestyle='--', label='E[L]')
plt.axvline(exact_var, color='orange', linestyle='--', label='VaR(L)')
plt.axvline(exact_cvar, color='red', linestyle='--', label='CVaR(L)')
plt.legend(fontsize=15)
plt.xlabel('Loss L ($)', size=15)
plt.ylabel('probability (%)', size=15)
plt.title('Loss Distribution', size=20)
plt.xticks(size=15)
plt.yticks(size=15)
plt.show()

# plot results for Z
plt.plot(z_values, p_z, 'o-', linewidth=3, markersize=8)
plt.grid()
plt.xlabel('Z value', size=15)
plt.ylabel('probability (%)', size=15)
plt.title('Z Distribution', size=20)
plt.xticks(size=15)
plt.yticks(size=15)
plt.show()

# plot results for default probabilities
plt.bar(range(K), p_default)
plt.xlabel('Asset', size=15)
plt.ylabel('probability (%)', size=15)
plt.title('Individual Default Probabilities', size=20)
plt.xticks(range(K), size=15)
plt.yticks(size=15)
plt.grid()
plt.show()

Perte attendue

Pour estimer la perte attendue, nous appliquons d’abord un opérateur de somme pondérée pour obtenir la perte totale à partir des pertes individuelles :

\[\mathcal{S}: |x_1, ..., x_K \rangle_K |0\rangle_{n_S} \mapsto |x_1, ..., x_K \rangle_K |\lambda_1x_1 + ... + \lambda_K x_K\rangle_{n_S}.\]

Le nombre qubits requis pour représenter le résultat est donné par

\[n_s = \lfloor \log_2( \lambda_1 + ... + \lambda_K ) \rfloor + 1.\]

Une fois que nous avons la distribution des pertes totales dans un registre quantique, nous pouvons utiliser les techniques décrites dans [Woerner2019] pour faire correspondre une perte totale :math:`L in { 0, …, 2 ^{n_s}-1 } ` à l’amplitude d’un qubit cible par un opérateur

\[ | L \rangle_{n_s}|0\rangle \mapsto | L \rangle_{n_s} \left( \sqrt{1 - L/(2^{n_s}-1)}|0\rangle + \sqrt{L/(2^{n_s}-1)}|1\rangle \right),\]

ce qui permet d’exécuter l’estimation de l’amplitude pour évaluer l’espérance de la perte.

# add Z qubits with weight/loss 0
from qiskit.circuit.library import WeightedAdder
agg = WeightedAdder(n_z + K, [0]*n_z + lgd)

from qiskit.circuit.library import LinearAmplitudeFunction

# define linear objective function
breakpoints = [0]
slopes = [1]
offsets = [0]
f_min = 0
f_max = sum(lgd)
c_approx = 0.25

objective = LinearAmplitudeFunction(
    agg.num_sum_qubits,
    slope=slopes,
    offset=offsets,
    # max value that can be reached by the qubit register (will not always be reached)
    domain=(0, 2**agg.num_sum_qubits-1),
    image=(f_min, f_max),
    rescaling_factor=c_approx,
    breakpoints=breakpoints
)

Création du circuit de préparation d’état :

# define the registers for convenience and readability
qr_state = QuantumRegister(u.num_qubits, 'state')
qr_sum = QuantumRegister(agg.num_sum_qubits, 'sum')
qr_carry = QuantumRegister(agg.num_carry_qubits, 'carry')
qr_obj = QuantumRegister(1, 'objective')

# define the circuit
state_preparation = QuantumCircuit(qr_state, qr_obj, qr_sum, qr_carry, name='A')

# load the random variable
state_preparation.append(u.to_gate(), qr_state)

# aggregate
state_preparation.append(agg.to_gate(), qr_state[:] + qr_sum[:] + qr_carry[:])

# linear objective function
state_preparation.append(objective.to_gate(), qr_sum[:] + qr_obj[:])

# uncompute aggregation
state_preparation.append(agg.to_gate().inverse(), qr_state[:] + qr_sum[:] + qr_carry[:])

# draw the circuit
state_preparation.draw()

             ┌───────┐┌────────┐      ┌───────────┐
    state_0: ┤0      ├┤0       ├──────┤0          ├
             │       ││        │      │           │
    state_1: ┤1      ├┤1       ├──────┤1          ├
             │  P(X) ││        │      │           │
    state_2: ┤2      ├┤2       ├──────┤2          ├
             │       ││        │      │           │
    state_3: ┤3      ├┤3       ├──────┤3          ├
             └───────┘│  adder │┌────┐│  adder_dg │
objective_0: ─────────┤        ├┤2   ├┤           ├
                      │        ││    ││           │
      sum_0: ─────────┤4       ├┤0 F ├┤4          ├
                      │        ││    ││           │
      sum_1: ─────────┤5       ├┤1   ├┤5          ├
                      │        │└────┘│           │
    carry_0: ─────────┤6       ├──────┤6          ├
                      └────────┘      └───────────┘

Avant d’utiliser QAE pour estimer la perte attendue, nous validons le circuit quantique représentant la fonction objectif en le simulant directement et en analysant la probabilité que le qubit objectif soit dans l’état \(|1\rangle\) c’est-à-dire que la valeur QAE sera approchée.

job = execute(state_preparation, backend=Aer.get_backend('statevector_simulator'))

# evaluate resulting statevector
value = 0
for i, a in enumerate(job.result().get_statevector()):
    b = ('{0:0%sb}' % (len(qr_state) + 1)).format(i)[-(len(qr_state) + 1):]
    am = np.round(np.real(a), decimals=4)
    if np.abs(am) > 1e-6 and b[0] == '1':
        value += am**2

print('Exact Expected Loss:   %.4f' % expected_loss)
print('Exact Operator Value:  %.4f' % value)
print('Mapped Operator value: %.4f' % objective.post_processing(value))

Exact Expected Loss:   0.6409
Exact Operator Value:  0.3906
Mapped Operator value: 0.6640

Ensuite, nous exécutons QAE pour estimer la perte attendue avec une accélération quadratique par rapport à la simulation classique de Monte Carlo.

# set target precision and confidence level
epsilon = 0.01
alpha = 0.05

# construct amplitude estimation
ae = IterativeAmplitudeEstimation(state_preparation=state_preparation,
                                  epsilon=epsilon, alpha=alpha,
                                  objective_qubits=[len(qr_state)],
                                  post_processing=objective.post_processing)
result = ae.run(quantum_instance=Aer.get_backend('qasm_simulator'), shots=100)

# print results
conf_int = np.array(result['confidence_interval'])
print('Exact value:    \t%.4f' % expected_loss)
print('Estimated value:\t%.4f' % result['estimation'])
print('Confidence interval: \t[%.4f, %.4f]' % tuple(conf_int))

Exact value:            0.6409
Estimated value:        0.6913
Confidence interval:    [0.6224, 0.7601]

Fonction de distribution cumulée

Au lieu de l’espérance de la perte (qui pourrait aussi être estimée efficacement à l’aide de techniques classiques), nous estimons à présent la fonction de distribution cumulative (CDF) de la perte. En termes classiques, il s’agit soit d’évaluer toutes les combinaisons possibles d’actifs défaillants, soit de nombreux échantillons classiques dans une simulation de Monte Carlo. Les algorithmes basés sur la QAE ont le potentiel d’accélérer considérablement cette analyse à l’avenir.

To estimate the CDF, i.e., the probability \(\mathbb{P}[L \leq x]\), we again apply \(\mathcal{S}\) to compute the total loss, and then apply a comparator that for a given value \(x\) acts as

\[\begin{split} \mathcal{C}: |L\rangle_n|0> \mapsto \begin{cases} |L\rangle_n|1> & \text{if}\quad L \leq x \\ |L\rangle_n|0> & \text{if}\quad L > x. \end{cases}\end{split}\]

L’état quantique obtenu peut être écrit comme ceci

\[ \sum_{L = 0}^{x} \sqrt{p_{L}}|L\rangle_{n_s}|1\rangle + \sum_{L = x+1}^{2^{n_s}-1} \sqrt{p_{L}}|L\rangle_{n_s}|1\rangle,\]

où nous adoptons directement la somme des valeurs de perte et des probabilités correspondantes au lieu de présenter les détails du modèle d’incertitude.

The CDF(\(x\)) equals the probability of measuring \(|1\rangle\) in the objective qubit and QAE can be directly used to estimate it.

# set x value to estimate the CDF
x_eval = 2

comparator = IntegerComparator(agg.num_sum_qubits, x_eval + 1, geq=False)
comparator.draw()

  state_0: ──■──────────────■──
             │              │
  state_1: ──┼────■─────────┼──
             │  ┌─┴─┐┌───┐  │
compare_0: ──┼──┤ X ├┤ X ├──┼──
           ┌─┴─┐└─┬─┘└───┘┌─┴─┐
     a0_0: ┤ X ├──■───────┤ X ├
           └───┘          └───┘

def get_cdf_circuit(x_eval):
    # define the registers for convenience and readability
    qr_state = QuantumRegister(u.num_qubits, 'state')
    qr_sum = QuantumRegister(agg.num_sum_qubits, 'sum')
    qr_carry = QuantumRegister(agg.num_carry_qubits, 'carry')
    qr_obj = QuantumRegister(1, 'objective')
    qr_compare = QuantumRegister(1, 'compare')

    # define the circuit
    state_preparation = QuantumCircuit(qr_state, qr_obj, qr_sum, qr_carry, name='A')

    # load the random variable
    state_preparation.append(u, qr_state)

    # aggregate
    state_preparation.append(agg, qr_state[:] + qr_sum[:] + qr_carry[:])

    # comparator objective function
    comparator = IntegerComparator(agg.num_sum_qubits, x_eval + 1, geq=False)
    state_preparation.append(comparator, qr_sum[:] + qr_obj[:] + qr_carry[:])

    # uncompute aggregation
    state_preparation.append(agg.inverse(), qr_state[:] + qr_sum[:] + qr_carry[:])

    return state_preparation

state_preparation = get_cdf_circuit(x_eval)

Encore une fois, nous utilisons d’abord la simulation quantique pour valider le circuit quantique.

job = execute(state_preparation, backend=Aer.get_backend('statevector_simulator'))

state_preparation.draw()

             ┌───────┐┌────────┐        ┌───────────┐
    state_0: ┤0      ├┤0       ├────────┤0          ├
             │       ││        │        │           │
    state_1: ┤1      ├┤1       ├────────┤1          ├
             │  P(X) ││        │        │           │
    state_2: ┤2      ├┤2       ├────────┤2          ├
             │       ││        │        │           │
    state_3: ┤3      ├┤3       ├────────┤3          ├
             └───────┘│  adder │┌──────┐│  adder_dg │
objective_0: ─────────┤        ├┤2     ├┤           ├
                      │        ││      ││           │
      sum_0: ─────────┤4       ├┤0     ├┤4          ├
                      │        ││  cmp ││           │
      sum_1: ─────────┤5       ├┤1     ├┤5          ├
                      │        ││      ││           │
    carry_0: ─────────┤6       ├┤3     ├┤6          ├
                      └────────┘└──────┘└───────────┘

# evaluate resulting statevector
var_prob = 0
for i, a in enumerate(job.result().get_statevector()):
    b = ('{0:0%sb}' % (len(qr_state) + 1)).format(i)[-(len(qr_state) + 1):]
    prob = np.abs(a)**2
    if prob > 1e-6 and b[0] == '1':
        var_prob += prob
print('Operator CDF(%s)' % x_eval + ' = %.4f' % var_prob)
print('Exact    CDF(%s)' % x_eval + ' = %.4f' % cdf[x_eval])

Operator CDF(2) = 0.9591
Exact    CDF(2) = 0.9591

Ensuite, nous exécutons QAE pour estimer le CDF pour un \(x\) donné.

# set target precision and confidence level
epsilon = 0.01
alpha = 0.05

# construct amplitude estimation
ae_cdf = IterativeAmplitudeEstimation(state_preparation=state_preparation,
                                      epsilon=epsilon, alpha=alpha,
                                      objective_qubits=[len(qr_state)])
result_cdf = ae_cdf.run(quantum_instance=Aer.get_backend('qasm_simulator'), shots=100)

# print results
conf_int = np.array(result_cdf['confidence_interval'])
print('Exact value:    \t%.4f' % cdf[x_eval])
print('Estimated value:\t%.4f' % result_cdf['estimation'])
print('Confidence interval: \t[%.4f, %.4f]' % tuple(conf_int))

Exact value:            0.9591
Estimated value:        0.9595
Confidence interval:    [0.9584, 0.9605]

Valeur à Risque

Dans ce qui suit, nous utilisons une recherche de bisection et QAE pour évaluer efficacement le CDF pour estimer la valeur à risque.

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

    # construct amplitude estimation
    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']

def bisection_search(objective, target_value, low_level, high_level, low_value=None, high_value=None):
    """
    Determines the smallest level such that the objective value is still larger than the target
    :param objective: objective function
    :param target: target value
    :param low_level: lowest level to be considered
    :param high_level: highest level to be considered
    :param low_value: value of lowest level (will be evaluated if set to None)
    :param high_value: value of highest level (will be evaluated if set to None)
    :return: dictionary with level, value, num_eval
    """

    # check whether low and high values are given and evaluated them otherwise
    print('--------------------------------------------------------------------')
    print('start bisection search for target value %.3f' % target_value)
    print('--------------------------------------------------------------------')
    num_eval = 0
    if low_value is None:
        low_value = objective(low_level)
        num_eval += 1
    if high_value is None:
        high_value = objective(high_level)
        num_eval += 1

    # check if low_value already satisfies the condition
    if low_value > target_value:
        return {'level': low_level, 'value': low_value, 'num_eval': num_eval, 'comment': 'returned low value'}
    elif low_value == target_value:
        return {'level': low_level, 'value': low_value, 'num_eval': num_eval, 'comment': 'success'}

    # check if high_value is above target
    if high_value < target_value:
        return {'level': high_level, 'value': high_value, 'num_eval': num_eval, 'comment': 'returned low value'}
    elif high_value == target_value:
        return {'level': high_level, 'value': high_value, 'num_eval': num_eval, 'comment': 'success'}

    # perform bisection search until
    print('low_level    low_value    level    value    high_level    high_value')
    print('--------------------------------------------------------------------')
    while high_level - low_level > 1:

        level = int(np.round((high_level + low_level) / 2.0))
        num_eval += 1
        value = objective(level)

        print('%2d           %.3f        %2d       %.3f    %2d            %.3f' \
              % (low_level, low_value, level, value, high_level, high_value))

        if value >= target_value:
            high_level = level
            high_value = value
        else:
            low_level = level
            low_value = value

    # return high value after bisection search
    print('--------------------------------------------------------------------')
    print('finished bisection search')
    print('--------------------------------------------------------------------')
    return {'level': high_level, 'value': high_value, 'num_eval': num_eval, 'comment': 'success'}

# run bisection search to determine VaR
objective = lambda x: run_ae_for_cdf(x)
bisection_result = bisection_search(objective, 1-alpha, min(losses)-1, max(losses), low_value=0, high_value=1)
var = bisection_result['level']

--------------------------------------------------------------------
start bisection search for target value 0.950
--------------------------------------------------------------------
low_level    low_value    level    value    high_level    high_value
--------------------------------------------------------------------
-1           0.000         1       0.751     3            1.000
 1           0.751         2       0.960     3            1.000
--------------------------------------------------------------------
finished bisection search
--------------------------------------------------------------------

print('Estimated Value at Risk: %2d' % var)
print('Exact Value at Risk:     %2d' % exact_var)
print('Estimated Probability:    %.3f' % bisection_result['value'])
print('Exact Probability:        %.3f' % cdf[exact_var])

Estimated Value at Risk:  2
Exact Value at Risk:      2
Estimated Probability:    0.960
Exact Probability:        0.959

Valeur à Risque Conditionnelle

Enfin, nous calculons le CVaR, c’est-à-dire que la valeur attendue de la perte est conditionnelle à ce qu’elle soit supérieure ou égale à la VaR. Pour ce faire, nous évaluons une fonction objectif linéaire par morceaux \(f(L)\), dépendant de la perte totale \(L\), qui est donnée par

\[\begin{split} f(L) = \begin{cases} 0 & \text{if}\quad L \leq VaR \\ L & \text{if}\quad L > VaR. \end{cases}\end{split}\]

To normalize, we have to divide the resulting expected value by the VaR-probability, i.e. \(\mathbb{P}[L \leq VaR]\).

# define linear objective
breakpoints = [0, var]
slopes = [0, 1]
offsets = [0, 0]  # subtract VaR and add it later to the estimate
f_min = 0
f_max = 3 - var
c_approx = 0.25

cvar_objective = LinearAmplitudeFunction(
    agg.num_sum_qubits,
    slopes,
    offsets,
    domain=(0, 2**agg.num_sum_qubits - 1),
    image=(f_min, f_max),
    rescaling_factor=c_approx,
    breakpoints=breakpoints
)

cvar_objective.draw()

        ┌─────────┐┌──────┐┌─────────┐┌─────────┐
q138_0: ┤0        ├┤0     ├┤0        ├┤0        ├
        │         ││      ││         ││         │
q138_1: ┤1 LinRot ├┤1     ├┤1 LinRot ├┤1        ├
        │         ││      ││         ││         │
q139_0: ┤2        ├┤  cmp ├┤2        ├┤  cmp_dg ├
        └─────────┘│      │└────┬────┘│         │
  a4_0: ───────────┤2     ├─────■─────┤2        ├
                   │      │           │         │
  a4_1: ───────────┤3     ├───────────┤3        ├
                   └──────┘           └─────────┘

# define the registers for convenience and readability
qr_state = QuantumRegister(u.num_qubits, 'state')
qr_sum = QuantumRegister(agg.num_sum_qubits, 'sum')
qr_carry = QuantumRegister(agg.num_carry_qubits, 'carry')
qr_obj = QuantumRegister(1, 'objective')
qr_work = QuantumRegister(cvar_objective.num_ancillas - len(qr_carry), 'work')

# define the circuit
state_preparation = QuantumCircuit(qr_state, qr_obj, qr_sum, qr_carry, qr_work, name='A')

# load the random variable
state_preparation.append(u, qr_state)

# aggregate
state_preparation.append(agg, qr_state[:] + qr_sum[:] + qr_carry[:])

# linear objective function
state_preparation.append(cvar_objective, qr_sum[:] + qr_obj[:] + qr_carry[:] + qr_work[:])

# uncompute aggregation
state_preparation.append(agg.inverse(), qr_state[:] + qr_sum[:] + qr_carry[:])

<qiskit.circuit.instructionset.InstructionSet at 0x7fdc6a68f110>

Encore une fois, nous utilisons d’abord la simulation quantique pour valider le circuit quantique.

job = execute(state_preparation, backend=Aer.get_backend('statevector_simulator'))

# evaluate resulting statevector
value = 0
for i, a in enumerate(job.result().get_statevector()):
    b = ('{0:0%sb}' % (len(qr_state) + 1)).format(i)[-(len(qr_state) + 1):]
    am = np.round(np.real(a), decimals=4)
    if np.abs(am) > 1e-6 and b[0] == '1':
        value += am**2

# normalize and add VaR to estimate
value = cvar_objective.post_processing(value)
d = (1.0 - bisection_result['value'])
v = value / d if d != 0 else 0
normalized_value = v + var
print('Estimated CVaR: %.4f' % normalized_value)
print('Exact CVaR:     %.4f' % exact_cvar)

Estimated CVaR: 3.3121
Exact CVaR:     3.0000

Ensuite, nous exécutons QAE pour estimer la CVaR.

# set target precision and confidence level
epsilon = 0.01
alpha = 0.05

# construct amplitude estimation
ae_cvar = IterativeAmplitudeEstimation(state_preparation=state_preparation,
                                       epsilon=epsilon, alpha=alpha,
                                       objective_qubits=[len(qr_state)],
                                       post_processing=cvar_objective.post_processing)
result_cvar = ae_cvar.run(quantum_instance=Aer.get_backend('qasm_simulator'), shots=100)

# print results
d = (1.0 - bisection_result['value'])
v = result_cvar['estimation'] / d if d != 0 else 0
print('Exact CVaR:    \t%.4f' % exact_cvar)
print('Estimated CVaR:\t%.4f' % (v + var))

Exact CVaR:     3.0000
Estimated CVaR: 3.2474

import qiskit.tools.jupyter
%qiskit_version_table
%qiskit_copyright

Version Information

Qiskit Software	Version
Qiskit	None
Terra	0.17.0.dev0+4ada179
Aer	0.6.1
Ignis	0.5.0.dev0+470d8cc
Aqua	0.9.0.dev0+5a88b59
IBM Q Provider	0.8.0
System information
Python	3.7.7 (default, May 6 2020, 04:59:01) [Clang 4.0.1 (tags/RELEASE_401/final)]
OS	Darwin
CPUs	2
Memory (Gb)	16.0
Tue Oct 20 11:03:45 2020 CEST

This code is a part of Qiskit

© Copyright IBM 2017, 2020.

This code is licensed under the Apache License, Version 2.0. You may
obtain a copy of this license in the LICENSE.txt file in the root directory
of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.

Any modifications or derivative works of this code must retain this
copyright notice, and modified files need to carry a notice indicating
that they have been altered from the originals.