注釈

当ページは tutorials/finance/09_credit_risk_analysis.ipynb. から生成されました。

IBM Quantum lab でインタラクティブに実行します。

信用リスク分析

はじめに

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

問題定義

このチュートリアルでは \(K\) 資産のポートフォリオの信用リスクを分析します。それぞれの資産 \(k\) のデフォルト確率は ガウス型条件付き独立 モデルに従います。すなわち、潜在乱数 \(Z\) からサンプルされたある値 \(z\) は標準正規分布に従い、資産:math:k のデフォルト確率は次で与えられます。

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

ここで \(F\) は \(Z\) の累積確率分布、\(p_k^0\) は資産 \(k\) におけるデフォルト確率、そして \(\rho_k\) は 資産 \(k\) の \(Z\) におけるデフォルト確率の感度です。こうして \(Z\) の具体形が与えられ、個別のデフォルトはそれぞれ独立とみなされます。

我々は総損失のリスク分析に興味があります。

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

ここで \(\lambda_k\) は資産 \(k\) の*デフォルトによる損失*を意味し、\(Z\), \(X_k(Z)\) は資産 \(k\) のデフォルト・イベントに対するベルヌーイ変数です。より正確には、我々は期待値 \(\mathbb{E}[L]\) 、\(L\) のバリュー・アット・リスク(VaR)と \(L\) の条件付きバリュー・アット・リスク(期待ショートフォールとも呼ばれる)に興味があります。ここで VaR と CVaR は次のように定義されます。

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

信頼レベル \(\alpha \in [0, 1]\)、かつ以下の通りです。

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

検討されたモデルの詳細については、例えば Regulatory Capital Modeling for Credit Risk. Marek Rutkowski, Silvio Tarca を参照してください。

この問題は次のパラメーターで定義されます。- \(n_z\) で表される \(Z\) を表現するための量子ビット数。 - \(z_{\text{max}}\) で表される \(Z\) の切り捨て値、すなわちZ は \(\{-z_{max}, ..., +z_{max}\}\) 内で、\(2^{n_z}\), \(k=1, ..., K\) で表されるそれぞれの資産のベース・デフォルト確率。 - \(\rho_k \in [0, 1)\) で表される、デフォルト確率の \(Z\) による感度。 - \(\lambda_k\) で表される、資産 \(k\) のデフォルト損失。 - \(\alpha \in [0, 1]\) VaR/CVaR の信頼レベル。

# 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

不確実性モデル

我々は不確実性モデルをロードする回路を構築します。これは標準正規分布に従う \(Z\) を表現する \(n_z\) 個の量子ビットの量子状態を作ることで達成されます。この状態は、\(K\) 量子ビットの二つ目の量子レジスター上で、単一量子ビットの Y 回転を制御するために使われます。ここで量子ビット \(k\) の \(|1\rangle\) 状態は資産 \(k\) のデフォルト・イベントを表します。その結果の量子状態は、次のように記述できます。

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

ここで \(z_i\) は、切り捨てられ丸められた \(Z\) の \(i\) 番目の値を表します。[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        ├
                         └─────────┘

シミュレーターを使って、\(|\Psi\rangle\) を構築する回路を検証し、次の厳密値を計算しましょう。- 期待損失 \(\mathbb{E}[L]\)。 - \(L\) のPDF と CDF。 - バリュー・アット・リスク \(VaR(L)\)。- 条件付きバリュー・アット・リスク \(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()

期待損失

期待損失を見積もるために、我々はまずはじめに加重合計演算子を適用し、個別損失を合計損失に足し上げます:

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

結果を表すために必要な量子ビット数は次の通りです。

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

量子レジスターの合計損失分布が分かると、 [Woerner2019] で述べられたテクニックを使って、合計損失 \(L \in \{0, ..., 2^{n_s}-1\}\) を演算子で目的量子ビットの増幅にマップすることができます。

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

これにより、期待損失を評価するために振幅推定を実行できます。

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

状態準備回路を作成します:

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

期待損失を見積もるために QAE を実行する前に、我々は目的関数の量子回路を直接シミュレートで検証し、目的量子ビットが \(|1\rangle\) 状態にいる確率を分析します。つまり、QAE の値は最終的に近似されます。

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

次に期待損失を推定するために、古典的モンテカルロ法よりも二次でスピードアップする QAE を実行します。

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

累積分布関数

(古典的な方法で効率的に推定できる) 期待損失の代わりに、損失の累積分布関数(CDF) を推定します。古典的には、資産デフォルトのすべての可能な組み合わせを評価するか、モンテカルロシミュレーションで多数の古典的サンプルを評価するか、どちらかを伴います。QAE に基づくアルゴリズムは将来大幅に解析をスピードアップする可能性があります。

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

結果の量子状態は次のように記述できます。

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

ここで、不確定性モデルの詳細を表す代わりに、合計された損失値とその確率を直接仮定します。

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)

再び、量子シミュレーションを使って量子回路を検証します。

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

次に、 QAE を実行して、与えられた \(x\) の CDF を推定します。

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

バリュー・アット・リスク

以下では、効果的に CDF を評価してバリュー・アット・リスクを推定するために、二分法検索と QAE を使います。

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

条件付きバリュー・アット・リスク

最後に、VaR 以上の値をとる条件付き損失の期待値CVaRを計算します。そうするために、次のように与えられる合計損失 \(L\) に依存する区分線形目的関数 \(f(L)\) を評価します。

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

再び、量子シミュレーションを使って量子回路を検証します。

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

次に、CVaRを推定するためにQAEを実行します。

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