Nota

Esta página foi gerada a partir de `tutorials/finance/09_credit_risk_analysis.ipynb `__.

Execute interativamente no IBM Quantum lab.

Análise de Risco de Crédito

Introdução

Este tutorial mostra como algoritmos quânticos podem ser usados para a análise de risco de crédito. Mais precisamente, como a Estimativa da Amplitude Quântica (QAE) pode ser usada para estimar medidas de risco com uma aceleração quadrática sobre a simulação clássica de Monte Carlo. O tutorial é baseado nos seguintes artigos:

• Análise de Risco Quântico. Stefan Woerner, Daniel J. Ogger. [Woerner2019]

• Análise de Risco de Crédito utilizando Computadores Quânticos. Egger et al. (2019) <https://arxiv.org/abs/1907.03044> __ [Egger2019]

Uma introdução geral ao QAE pode ser encontrada no seguinte artigo:

• Amplificação de Amplitude Quântica e Estimativa. Gilles Brassard et al.

A estrutura do tutorial é a seguinte:

1. Definição do Problema

2. Modelo de Incerteza

3. Perda Esperada

4. Função de Distribuição Cumulativa

5. `Valor em Risco <#Value-at-Risk> ` __

6. Valor condicional em risco

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

Definição do Problema

Neste tutorial queremos analisar o risco de crédito de uma carteira de ativos de :math:` K . A probabilidade padrão de cada ativo :math: k ` segue um modelo *Gaussiana Condicional de Independência *, ou seja, dado um valor :math:` z ` exemplificado a partir de uma variável aleatória latente :math:` Z ` seguindo uma distribuição normal padrão, a probabilidade padrão do ativo :math:` k ` é dada por

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

onde \(F\) denota a função de distribuição cumulativa de \(Z\), \(p_k^0\) é a probabilidade padrão do ativo \(k\) para \(z=0\) e \(\rho_k\) é a sensibilidade da probabilidade padrão do ativo \(k\) com relação a \(Z\). Assim, dada uma realização concreta de \(Z\) os eventos individuais padrões assumem-se como independentes uns dos outros.

Estamos interessados em analisar medições de risco de perda total

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

onde \(\lambda_k\) denota a perda dada padrão do ativo \(k\), e dada \(Z\), \(X_k (Z)\) denota uma variável Bernoulli representando o evento padrão do ativo \(k\). Mais precisamente, estamos interessados no valor esperado \(\mathbb{E}[L]\), o Valor em Risco (VaR) de \(L\) e o Valor Condicional em Risco de \(L\) (também chamado de Shortfall Esperado). Onde VaR e CVaR são definidos como

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

com nível de confiança \(\alpha \in [0, 1]\), e

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

Para obter mais detalhes sobre o modelo considerado, veja, por exemplo, Regulatory Capital Modeling for Credit Risk. Marek Rutkowski, Silvio Tarca

O problema é definido pelos seguintes parâmetros: - número de qubits usados para representar \(Z\), denotado por \(n_z\) - valor de truncamento para \(Z\), denotado por \(z_{\text{max}}\), i. ., Z presume-se que tenha \(2^{n_z}\) valores equidistantes em \(\{-z_{max}, ... +z_{max}\}\) - as probabilidades padrões base para cada ativo \(p_0^k \in (0, 1)\), \(k=1, ... K\) - sensibilidades das probabilidades padrão em relação a \(Z\), denotado por \(\rho_k \in [0, 1)\) - perda padrão dada para o ativo \(k\), denotado por \(\lambda_k\) - nível de confiança do 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

Modelo de incerteza

Nós agora construímos um circuito que carrega o modelo de incerteza. Isto pode ser alcançado através da criação de um estado quântico num registo de \(n_z\) qubits que representa \(Z\) seguindo uma distribuição normal padrão. Este estado é então utilizado para controlar o único qubit Y-rotações em um segundo registo de qubit de \(K\) qubits, onde um estado :math:` | 1rangle` do qubit \(k\) representa o evento padrão do ativo \(k\). O estado quântico resultante pode ser escrito como

\[ |\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),\]

onde denotamos por \(z_i\) the \(i\)-th valor discreto e truncado \(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        ├
                         └─────────┘

Utilizamos agora o simulador para validar o circuito que constrói \(|\Psi\rangle\) e computar os valores exatos correspondentes para - perda esperada \(\mathbb{E}[L]\) - PDF e CDF de \(L\)- valor em risco \(VaR (L)\) e probabilidade correspondente - valor condicional em risco \(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()

Perda Esperada

Para estimar a perda esperada, primeiramente aplicamos um operador de soma ponderada para somar as perdas individuais à perda total:

\[\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}.\]

O número necessário de qubits para representar o resultado é dado por

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

Uma vez que temos a distribuição total de perda em um registro quântico, podemos utilizar as técnicas descritas em [Woerner2019] para mapear uma perda total \(L \in \ {0, ..., 2^{n_s}-1 \}\) à amplitude de um qubit objetivo por um operador

\[ | 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),\]

que permite executar uma amplitude estimada, do valor da perda esperada.

# 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
)

Crie o circuito de preparação de estado:

# 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          ├
                      └────────┘      └───────────┘

Antes de utilizarmos o QAE para estimar a perda esperada, validamos o circuito quântico representando a função objetiva apenas simulando-a diretamente e analisando a probabilidade de o qubit objetivo estar no estado \(|1\rangle\), ou seja, o valor QAE eventualmente se aproximará.

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

A seguir, executaremos QAE para estimar a perda esperada, com uma aceleração quadrática sobre a simulador de Monte Carlo clássico.

# 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]

Função De Distribuição Cumulativa

Em vez da perda esperada (que também poderia ser estimada eficientemente utilizando técnicas clássicas) estimamos agora a função de distribuição cumulativa (CDF) da perda. Classicamente, isso ou envolve a avaliação de todas as combinações possíveis de ativos padrões, ou muitas amostras clássicas em uma simulação de Monte Carlo. Os algoritmos baseados em QAE têm o potencial de acelerar significativamente essa análise no futuro.

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}\]

O estado quântico resultante pode ser escrito como

\[ \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,\]

onde assumimos diretamente os valores somados de perda e as probabilidades correspondentes em vez de apresentar os detalhes do modelo de incerteza.

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)

Novamente, primeiro usaremos o simulador quântico para validar o circuito quântico.

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

A seguir, executaremos QAE para estimar o CDF para um dado \(x\).

# 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]

Valor em Risco

A seguir usaremos um busca bipartida e QAE para avaliar eficientemente o CDF e estimar o valor de risco.

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

Valor condicional em Risco

Por último, computamos o CVaR, ou seja, o valor esperado da perda condicionado a ele ser maior ou igual ao VaR. Para isso, avaliamos uma função objetiva linear por partes \(f(L)\), dependente da perda total \(L\), que é dada por

\[\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>

Novamente, primeiro usaremos o simulador quântico para validar o circuito quântico.

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

A seguir executaremos QAE para estimar o 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.