From 2d03960fbf032fc3be9157de10133ee7609f77b6 Mon Sep 17 00:00:00 2001 From: CREVIOS <48938983+CREVIOS@users.noreply.github.com> Date: Mon, 3 Aug 2026 16:00:05 +0600 Subject: [PATCH] contest2: quantum risk analysis application (OriginQCup_DU_Fanta) Adds a complete Value at Risk and Expected Shortfall workflow built from two components that already exist in the library: QCmp.int_comparator supplies the reversible predicate L >= t, and QAE.IQAE estimates the probability that it holds. The tail probability P(L >= t) therefore costs O(1/eps) oracle queries instead of the O(1/eps^2) samples a classical Monte Carlo simulation needs. The design point is that both risk measures reduce to tail probabilities, so no second estimation primitive is required. For an integer-valued loss E[(L - t)^+] = sum_{k > t} P(L >= k) ES_alpha = VaR_alpha + E[(L - VaR_alpha)^+] / P(L >= VaR_alpha) which makes Value at Risk a bisection over tail probabilities, costing O(log N) amplitude estimations rather than the O(N) of a linear scan, and Expected Shortfall a short sum of them. Both identities are verified against exact values in the tests. LossDistribution.from_credit_portfolio builds the exact loss distribution of a portfolio of independent obligors by convolution; the test suite checks it against brute-force enumeration of all default combinations. On the demo portfolio (five obligors, 8 qubits) the estimated Value at Risk matches the exact value at confidence levels 0.90, 0.95 and 0.99, and Expected Shortfall is within 0.07 loss levels, using 16 to 18 amplitude estimations. The VaR bisection compares tail probabilities against 1 - alpha, so the estimation accuracy has to be small relative to that gap. At alpha=0.99 with eps=0.01 the true P(L >= 5) = 0.01056 lies within eps of the 0.01 decision threshold and the level is under-estimated. recommended_epsilon returns 0.2 * (1 - alpha) and value_at_risk warns when eps is too coarse. Contents: - quantum_risk.py core module - example_credit_risk.py credit portfolio demonstration - Test_quantum_risk.py 26 unit tests, all against exact classical values - README.md bilingual documentation No end-to-end quantum speedup is claimed: the quadratic advantage applies to the query complexity of the tail-probability subroutine only. QAE.IQAE currently supports the CPU simulator backend, so no hardware run is included. --- contest2/OriginQCup_DU_Fanta/README.md | 164 ++++++ .../OriginQCup_DU_Fanta/Test_quantum_risk.py | 205 ++++++++ .../example_credit_risk.py | 161 ++++++ contest2/OriginQCup_DU_Fanta/quantum_risk.py | 479 ++++++++++++++++++ 4 files changed, 1009 insertions(+) create mode 100644 contest2/OriginQCup_DU_Fanta/README.md create mode 100644 contest2/OriginQCup_DU_Fanta/Test_quantum_risk.py create mode 100644 contest2/OriginQCup_DU_Fanta/example_credit_risk.py create mode 100644 contest2/OriginQCup_DU_Fanta/quantum_risk.py diff --git a/contest2/OriginQCup_DU_Fanta/README.md b/contest2/OriginQCup_DU_Fanta/README.md new file mode 100644 index 00000000..0929362c --- /dev/null +++ b/contest2/OriginQCup_DU_Fanta/README.md @@ -0,0 +1,164 @@ +# 量子风险分析：基于振幅估计的 VaR 与预期损失 + +**参赛队伍**：DU_Fanta ｜ **任务类型**：创新应用 ｜ **关联 Issue**：#13 + +--- + +## 一、这个应用做什么 + +在险价值（Value at Risk, VaR）和预期损失（Expected Shortfall, ES）是巴塞尔协议下最核心的两个风险指标。银行通常用蒙特卡洛模拟计算它们：要把误差降到 `eps`，需要 `O(1/eps^2)` 次采样。 + +本应用把 `pyqpanda_alg` 中已有的两个组件组合成一条完整的风险计量流水线： + +| 组件 | 作用 | +| --- | --- | +| `QCmp.int_comparator` | 构造可逆判定电路 `L >= t` | +| `QAE.IQAE` | 估计该判定成立的概率，即尾部概率 `P(L >= t)` | + +迭代振幅估计只需 `O(1/eps)` 次 oracle 查询，相对经典蒙特卡洛在**该子过程**上具有平方级优势。 + +### 关键设计：整条流水线只依赖尾部概率 + +本应用的核心观察是：**VaR 与 ES 都可以只用尾部概率表示**，因此一个量子原语就足够，不需要为 ES 单独设计带线性幅度旋转的电路。 + +对整数取值的损失 `L`： + +``` +E[(L - t)^+] = sum_{k > t} P(L >= k) + +ES_alpha = VaR_alpha + E[(L - VaR_alpha)^+] / P(L >= VaR_alpha) +``` + +于是： + +- **VaR** ：对 `t` 做二分搜索，只需 `O(log N)` 次振幅估计（而非线性扫描的 `O(N)` 次）； +- **ES** ：VaR 之上若干个尾部概率求和即可。 + +两个恒等式均在 `Test_quantum_risk.py` 中对照精确值验证。 + +--- + +## 二、文件说明 + +| 文件 | 内容 | +| --- | --- | +| `quantum_risk.py` | 核心模块：`LossDistribution`、`QuantumRiskAnalyzer` | +| `example_credit_risk.py` | 信贷组合的完整演示 | +| `Test_quantum_risk.py` | 26 个单元测试，全部对照经典精确值 | + +--- + +## 三、快速开始 + +```bash +pip install pyqpanda3 pyqpanda_alg numpy +cd contest2/OriginQCup_DU_Fanta +python example_credit_risk.py # 运行演示 +python -m pytest Test_quantum_risk.py # 运行测试 +``` + +最小示例： + +```python +from quantum_risk import LossDistribution, QuantumRiskAnalyzer, recommended_epsilon + +# 5 个债务人，各自的违约概率与整数敞口 +dist = LossDistribution.from_credit_portfolio( + default_probabilities=[0.08, 0.15, 0.04, 0.22, 0.10], + exposures=[2, 1, 3, 1, 2], + unit=1e5, +) + +confidence = 0.99 +analyzer = QuantumRiskAnalyzer(dist, epsilon=recommended_epsilon(confidence)) +summary = analyzer.report(confidence) + +print(summary['var_monetary']) # 500000.0 +print(summary['expected_shortfall_monetary']) # 约 530000 +``` + +--- + +## 四、实测结果 + +演示所用组合：5 个债务人，总敞口 90 万元，期望损失 8.5 万元；损失分布占 4 个量子比特（16 个取值），加上比较器共 **8 个量子比特**。 + +### 尾部概率（量子估计 vs 精确值，eps = 0.01） + +| 阈值 | 精确值 | 量子估计 | 绝对误差 | +| ---: | ---: | ---: | ---: | +| 1 | 0.47299 | 0.48212 | 0.00913 | +| 2 | 0.23135 | 0.23830 | 0.00695 | +| 3 | 0.10074 | 0.10271 | 0.00198 | +| 4 | 0.03092 | 0.02816 | 0.00276 | +| 5 | 0.01056 | 0.01015 | 0.00041 | +| 6 | 0.00278 | 0.00247 | 0.00032 | + +### 风险指标 + +| 置信度 | eps | VaR（量子） | VaR（精确） | ES（量子） | ES（精确） | ES 误差 | 振幅估计次数 | +| ---: | ---: | ---: | ---: | ---: | ---: | ---: | ---: | +| 0.90 | 0.0100 | 3 | 3 | 3.3785 | 3.4459 | 0.0673 | 18 | +| 0.95 | 0.0100 | 3 | 3 | 3.4660 | 3.4459 | 0.0201 | 18 | +| 0.99 | 0.0020 | 5 | 5 | 5.2986 | 5.3256 | 0.0271 | 16 | + +三个置信度下 VaR 均与精确值完全一致。 + +### 查询次数对比 + +| 精度 eps | 经典蒙特卡洛采样数 | +| ---: | ---: | +| 0.05 | 738 | +| 0.01 | 18,445 | +| 0.005 | 73,778 | +| 0.001 | 1,844,440 | + +经典方法为 `O(1/eps^2)`，迭代振幅估计为 `O(1/eps)`。 + +--- + +## 五、一个必须说明的精度条件 + +VaR 的二分搜索需要判断 `P(L >= t+1) <= 1 - alpha`。**如果估计精度 `eps` 与 `1 - alpha` 相当，这个比较就会被估计噪声主导**，返回的 VaR 可能相差一档。 + +开发过程中在 `alpha = 0.99`、`eps = 0.01` 下确实观察到了该现象：真实 `P(L >= 5) = 0.01056`，与判定阈值 `0.01` 的差距小于 `eps`，导致 VaR 被低估为 4。 + +因此模块提供 `recommended_epsilon(alpha) = 0.2 * (1 - alpha)`，并在 `eps` 过粗时发出 `RuntimeWarning`。上表中 `alpha = 0.99` 使用 `eps = 0.002` 后结果正确。 + +--- + +## 六、本应用不主张的内容 + +- **不主张端到端的量子加速。** 平方级优势存在于尾部概率这一子过程的查询复杂度上。整条流水线还包含经典的分布构造与二分搜索；在模拟器上，经典计算本身更快。 +- **不主张已可用于生产。** 演示规模为 4 个损失比特（16 个取值）。真实组合需要更多比特、相关性违约模型，以及在含噪声硬件上的误差缓解。 +- **未在真实量子硬件上运行。** `QAE.IQAE` 目前仅支持 CPU 模拟器后端。 + +--- + +## English summary + +**Quantum risk analysis: Value at Risk and Expected Shortfall via amplitude estimation.** + +This application composes two existing `pyqpanda_alg` components — `QCmp.int_comparator` (the reversible predicate `L >= t`) and `QAE.IQAE` (amplitude estimation) — into a complete risk-measurement workflow. + +The key design point is that **both VaR and Expected Shortfall reduce to tail probabilities**, so one quantum primitive suffices. For integer-valued losses: + +``` +E[(L - t)^+] = sum_{k > t} P(L >= k) +ES_alpha = VaR_alpha + E[(L - VaR_alpha)^+] / P(L >= VaR_alpha) +``` + +VaR therefore comes from a bisection over tail probabilities in `O(log N)` amplitude estimations instead of an `O(N)` scan, and ES from a short sum of them. No separate estimation primitive with linear-amplitude rotations is needed. + +**Results** on a 5-obligor credit portfolio (8 qubits: 4 loss + 4 comparator): VaR matches the exact value at all three confidence levels, Expected Shortfall to within 0.07 loss levels, using 16–18 amplitude estimations. + +**Accuracy condition.** The VaR bisection compares tail probabilities against `1 - alpha`, so `epsilon` must be small relative to that gap. This was observed concretely at `alpha = 0.99, epsilon = 0.01`, where the true `P(L >= 5) = 0.01056` sits within `epsilon` of the `0.01` decision threshold and VaR was under-estimated. The module exposes `recommended_epsilon(alpha) = 0.2 * (1 - alpha)` and warns when `epsilon` is too coarse. + +**Not claimed:** end-to-end quantum speedup (the quadratic advantage is in the query complexity of the tail-probability subroutine only), production readiness, or any hardware run — `QAE.IQAE` currently supports the CPU simulator backend only. + +**Tests:** 26 unit tests, every quantum estimate checked against the exact classical value. + +## References + +1. S. Woerner, D. J. Egger. *Quantum risk analysis.* npj Quantum Information **5**, 15 (2019). +2. D. Grinko, J. Gacon, C. Zoufal, S. Woerner. *Iterative quantum amplitude estimation.* npj Quantum Information **7**, 52 (2021). diff --git a/contest2/OriginQCup_DU_Fanta/Test_quantum_risk.py b/contest2/OriginQCup_DU_Fanta/Test_quantum_risk.py new file mode 100644 index 00000000..ec6a4980 --- /dev/null +++ b/contest2/OriginQCup_DU_Fanta/Test_quantum_risk.py @@ -0,0 +1,205 @@ +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Tests for the quantum risk analysis application. + +Every quantum estimate is checked against the exact value computed classically +from the same distribution, so the tests verify the workflow rather than merely +exercising it. +""" + +import itertools + +import numpy as np +import pytest + +from quantum_risk import ( + LossDistribution, + QuantumRiskAnalyzer, + classical_monte_carlo_samples, + recommended_epsilon, +) + +SIMPLE_PROBABILITIES = [0.30, 0.25, 0.18, 0.12, 0.08, 0.04, 0.02, 0.01] +PORTFOLIO_DEFAULTS = [0.10, 0.20, 0.05, 0.30] +PORTFOLIO_EXPOSURES = [1, 2, 1, 3] + + +class Test_LossDistribution: + + def test_normalises_and_pads_to_power_of_two(self): + dist = LossDistribution([1.0, 1.0, 1.0]) + assert dist.probabilities.size == 4, "概率向量应补齐到 2 的幂" + assert dist.num_qubits == 2, "3 个取值需要 2 个量子比特" + assert abs(dist.probabilities.sum() - 1.0) < 1e-12, "概率应归一化" + assert dist.probabilities[3] == 0.0, "补齐位置的概率应为 0" + + def test_tail_probability_matches_direct_sum(self): + dist = LossDistribution(SIMPLE_PROBABILITIES) + probabilities = np.asarray(SIMPLE_PROBABILITIES) / sum(SIMPLE_PROBABILITIES) + for threshold in range(len(SIMPLE_PROBABILITIES) + 1): + expected = probabilities[threshold:].sum() + assert abs(dist.tail_probability(threshold) - expected) < 1e-12, \ + f"阈值 {threshold} 的尾部概率不正确" + + def test_expected_shortfall_matches_conditional_mean(self): + dist = LossDistribution(SIMPLE_PROBABILITIES) + for confidence in (0.5, 0.8, 0.9, 0.95): + var = dist.value_at_risk(confidence) + mask = dist.levels >= var + expected = (dist.levels[mask] * dist.probabilities[mask]).sum() / \ + dist.tail_probability(var) + assert abs(dist.expected_shortfall(confidence) - expected) < 1e-12, \ + f"置信度 {confidence} 的期望损失不等于条件均值" + + def test_credit_portfolio_matches_brute_force(self): + """卷积得到的组合损失分布应与穷举所有违约组合一致。""" + dist = LossDistribution.from_credit_portfolio(PORTFOLIO_DEFAULTS, PORTFOLIO_EXPOSURES) + brute = np.zeros(sum(PORTFOLIO_EXPOSURES) + 1) + for indicators in itertools.product([0, 1], repeat=len(PORTFOLIO_DEFAULTS)): + probability, loss = 1.0, 0 + for indicator, default, exposure in zip(indicators, PORTFOLIO_DEFAULTS, + PORTFOLIO_EXPOSURES): + probability *= default if indicator else (1.0 - default) + loss += exposure * indicator + brute[loss] += probability + padded = np.zeros(dist.probabilities.size) + padded[:brute.size] = brute + assert np.max(np.abs(dist.probabilities - padded)) < 1e-12, "组合损失分布与穷举结果不一致" + + def test_state_preparation_reproduces_the_distribution(self): + """态制备电路产生的测量概率应等于目标分布。""" + from pyqpanda3.core import CPUQVM, QProg + + dist = LossDistribution(SIMPLE_PROBABILITIES) + qubits = list(range(dist.num_qubits)) + prog = QProg() + prog << dist.state_preparation(qubits) + machine = CPUQVM() + machine.run(prog, 1) + amplitudes = np.asarray(machine.result().get_state_vector()) + measured = np.abs(amplitudes) ** 2 + assert np.max(np.abs(measured - dist.probabilities)) < 1e-9, "态制备结果与目标分布不符" + + @pytest.mark.parametrize("probabilities", [[], [-0.1, 0.5], [0.0, 0.0]]) + def test_rejects_invalid_probabilities(self, probabilities): + with pytest.raises(ValueError): + LossDistribution(probabilities) + + def test_rejects_mismatched_portfolio_lengths(self): + with pytest.raises(ValueError): + LossDistribution.from_credit_portfolio([0.1, 0.2], [1]) + + @pytest.mark.parametrize("confidence", [0.0, 1.0, -0.5, 1.5]) + def test_rejects_invalid_confidence(self, confidence): + dist = LossDistribution(SIMPLE_PROBABILITIES) + with pytest.raises(ValueError): + dist.value_at_risk(confidence) + + +class Test_QuantumRiskAnalyzer: + + def test_tail_probability_matches_exact_value(self): + """振幅估计得到的尾部概率应落在给定精度内。""" + dist = LossDistribution(SIMPLE_PROBABILITIES) + epsilon = 0.01 + analyzer = QuantumRiskAnalyzer(dist, epsilon=epsilon) + for threshold in range(1, dist.probabilities.size): + exact = dist.tail_probability(threshold) + estimate = analyzer.tail_probability(threshold) + assert abs(estimate - exact) < 10 * epsilon, \ + f"阈值 {threshold}: 估计值 {estimate:.5f} 与精确值 {exact:.5f} 偏差过大" + + def test_tail_probability_boundaries_need_no_estimation(self): + dist = LossDistribution(SIMPLE_PROBABILITIES) + analyzer = QuantumRiskAnalyzer(dist, epsilon=0.05) + assert analyzer.tail_probability(0) == 1.0, "阈值为 0 时尾部概率必为 1" + assert analyzer.tail_probability(dist.probabilities.size) == 0.0, \ + "阈值超出取值范围时尾部概率必为 0" + assert analyzer.estimate_calls == 0, "边界情形不应调用振幅估计" + + def test_value_at_risk_matches_exact_value(self): + dist = LossDistribution(SIMPLE_PROBABILITIES) + for confidence in (0.80, 0.90): + analyzer = QuantumRiskAnalyzer(dist, epsilon=recommended_epsilon(confidence)) + assert analyzer.value_at_risk(confidence) == dist.value_at_risk(confidence), \ + f"置信度 {confidence} 下的 VaR 与精确值不一致" + + def test_value_at_risk_uses_logarithmic_number_of_estimations(self): + """二分搜索的振幅估计次数应为 O(log2(取值数))，而非线性扫描。""" + dist = LossDistribution(SIMPLE_PROBABILITIES) + analyzer = QuantumRiskAnalyzer(dist, epsilon=0.02) + analyzer.estimate_calls = 0 + analyzer.value_at_risk(0.80) + levels = dist.probabilities.size + assert analyzer.estimate_calls <= int(np.ceil(np.log2(levels))) + 1, \ + f"{levels} 个取值上的二分搜索调用了 {analyzer.estimate_calls} 次振幅估计" + + def test_expected_shortfall_matches_exact_value(self): + dist = LossDistribution(SIMPLE_PROBABILITIES) + confidence = 0.80 + analyzer = QuantumRiskAnalyzer(dist, epsilon=recommended_epsilon(confidence)) + exact = dist.expected_shortfall(confidence) + estimate = analyzer.expected_shortfall(confidence) + assert abs(estimate - exact) < 0.5, \ + f"期望损失估计 {estimate:.4f} 与精确值 {exact:.4f} 偏差过大" + + def test_report_contains_quantum_and_exact_values(self): + dist = LossDistribution(SIMPLE_PROBABILITIES) + confidence = 0.80 + analyzer = QuantumRiskAnalyzer(dist, epsilon=recommended_epsilon(confidence)) + summary = analyzer.report(confidence) + for key in ('var_quantum', 'var_exact', 'expected_shortfall_quantum', + 'expected_shortfall_exact', 'amplitude_estimations', 'num_qubits'): + assert key in summary, f"报告缺少字段 {key}" + assert summary['num_qubits'] == 2 * dist.num_qubits, "寄存器宽度应为损失寄存器的两倍" + assert summary['amplitude_estimations'] > 0, "报告应记录振幅估计的调用次数" + + def test_warns_when_epsilon_cannot_resolve_the_confidence_level(self): + """精度不足以分辨 1 - alpha 时应给出警告。""" + dist = LossDistribution(SIMPLE_PROBABILITIES) + analyzer = QuantumRiskAnalyzer(dist, epsilon=0.05) + with pytest.warns(RuntimeWarning): + analyzer.value_at_risk(0.99) + + def test_rejects_invalid_arguments(self): + dist = LossDistribution(SIMPLE_PROBABILITIES) + with pytest.raises(TypeError): + QuantumRiskAnalyzer(SIMPLE_PROBABILITIES) + with pytest.raises(ValueError): + QuantumRiskAnalyzer(dist, epsilon=0.0) + with pytest.raises(ValueError): + QuantumRiskAnalyzer(dist, epsilon=1.5) + + +class Test_query_complexity: + + def test_recommended_epsilon_scales_with_the_tail(self): + assert recommended_epsilon(0.90) > recommended_epsilon(0.99), \ + "置信度越高，所需精度应越小" + assert abs(recommended_epsilon(0.99) - 0.2 * 0.01) < 1e-12 + + def test_classical_sample_count_grows_quadratically(self): + coarse = classical_monte_carlo_samples(0.01) + fine = classical_monte_carlo_samples(0.005) + ratio = fine / coarse + assert 3.5 < ratio < 4.5, \ + f"精度提高一倍时经典采样量应约变为 4 倍，实际比值为 {ratio:.2f}" + + @pytest.mark.parametrize("epsilon", [0.0, 1.0, -0.1]) + def test_classical_sample_count_rejects_invalid_accuracy(self, epsilon): + with pytest.raises(ValueError): + classical_monte_carlo_samples(epsilon) + + +if __name__ == "__main__": + pytest.main([__file__, "-v"]) diff --git a/contest2/OriginQCup_DU_Fanta/example_credit_risk.py b/contest2/OriginQCup_DU_Fanta/example_credit_risk.py new file mode 100644 index 00000000..71bcf532 --- /dev/null +++ b/contest2/OriginQCup_DU_Fanta/example_credit_risk.py @@ -0,0 +1,161 @@ +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Credit portfolio risk analysis with quantum amplitude estimation. + +Runs the full workflow on a small loan book and prints every quantum estimate +next to the exact value, so the accuracy of the method is visible rather than +asserted. + +Usage + python example_credit_risk.py +""" + +import numpy as np + +from quantum_risk import ( + LossDistribution, + QuantumRiskAnalyzer, + classical_monte_carlo_samples, + recommended_epsilon, +) + +# A small loan book: five obligors with individual default probabilities and +# integer exposures in units of 100k CNY. +DEFAULT_PROBABILITIES = [0.08, 0.15, 0.04, 0.22, 0.10] +EXPOSURES = [2, 1, 3, 1, 2] +UNIT = 1e5 +EPSILON = 0.01 +CONFIDENCE_LEVELS = [0.90, 0.95, 0.99] + + +def print_portfolio(): + print('=' * 72) + print('Credit portfolio') + print('=' * 72) + print(f"{'obligor':>8} {'default prob':>14} {'exposure':>10} {'exposure (CNY)':>16}") + for index, (probability, exposure) in enumerate(zip(DEFAULT_PROBABILITIES, EXPOSURES), 1): + print(f'{index:>8} {probability:>14.2%} {exposure:>10} {exposure * UNIT:>16,.0f}') + expected = sum(p * e for p, e in zip(DEFAULT_PROBABILITIES, EXPOSURES)) + print(f"\ntotal exposure : {sum(EXPOSURES) * UNIT:,.0f} CNY") + print(f'expected loss : {expected * UNIT:,.0f} CNY') + + +def print_distribution(distribution): + print() + print('=' * 72) + print(f'Loss distribution on {distribution.num_qubits} qubits ' + f'({distribution.probabilities.size} levels)') + print('=' * 72) + print(f"{'loss':>6} {'CNY':>12} {'P(L = l)':>12} {'P(L >= l)':>12}") + for level, probability in enumerate(distribution.probabilities): + if probability < 1e-6: + continue + print(f'{level:>6} {level * UNIT:>12,.0f} {probability:>12.5f} ' + f'{distribution.tail_probability(level):>12.5f}') + + +def print_tail_probabilities(distribution, analyzer): + print() + print('=' * 72) + print('Tail probabilities: quantum amplitude estimation vs exact') + print('=' * 72) + print(f"{'threshold':>10} {'exact':>10} {'quantum':>10} {'abs error':>11}") + errors = [] + for threshold in range(1, distribution.probabilities.size): + exact = distribution.tail_probability(threshold) + if exact < 1e-6: + continue + estimate = analyzer.tail_probability(threshold) + errors.append(abs(estimate - exact)) + print(f'{threshold:>10} {exact:>10.5f} {estimate:>10.5f} {abs(estimate - exact):>11.5f}') + print(f'\nmax absolute error {max(errors):.5f}, target accuracy {EPSILON}') + + +def print_risk_measures(analyzer): + print() + print('=' * 72) + print('Risk measures') + print('=' * 72) + # Value at Risk is decided by comparing tail probabilities against 1 - alpha, + # so each confidence level needs its own accuracy; a coarse epsilon cannot + # resolve a 1% tail. Each report is also estimated once and reused, because + # amplitude estimation is randomised. + summaries = [] + for confidence in CONFIDENCE_LEVELS: + epsilon = min(EPSILON, recommended_epsilon(confidence)) + summaries.append(QuantumRiskAnalyzer(analyzer.distribution, epsilon=epsilon) + .report(confidence)) + + print(f"{'alpha':>7} {'epsilon':>9} {'VaR (q)':>9} {'VaR (exact)':>12} {'ES (q)':>9} " + f"{'ES (exact)':>11} {'ES error':>10} {'AE calls':>9}") + for summary in summaries: + print(f"{summary['confidence']:>7.2f} {summary['epsilon']:>9.4f} " + f"{summary['var_quantum']:>9} " + f"{summary['var_exact']:>12} " + f"{summary['expected_shortfall_quantum']:>9.4f} " + f"{summary['expected_shortfall_exact']:>11.4f} " + f"{summary['expected_shortfall_abs_error']:>10.4f} " + f"{summary['amplitude_estimations']:>9}") + + print() + for summary in summaries: + print(f"alpha = {summary['confidence']:.0%}: " + f"VaR = {summary['var_monetary']:,.0f} CNY, " + f"ES = {summary['expected_shortfall_monetary']:,.0f} CNY") + + boundary = [s for s in summaries if s['var_quantum'] != s['var_exact']] + if boundary: + print('\nNote: Value at Risk is an integer level, so a confidence level that') + print('falls very close to a jump of the cumulative distribution can round to') + print('either side. The levels affected here are: ' + + ', '.join(f"{s['confidence']:.0%}" for s in boundary) + '.') + + +def print_query_comparison(distribution): + print() + print('=' * 72) + print('Oracle queries vs classical Monte Carlo samples') + print('=' * 72) + print(f"{'accuracy':>10} {'classical samples':>19} {'ratio to 1/eps':>16}") + for epsilon in (0.05, 0.01, 0.005, 0.001): + samples = classical_monte_carlo_samples(epsilon) + print(f'{epsilon:>10} {samples:>19,} {samples * epsilon:>16,.1f}') + print('\nClassical Monte Carlo needs O(1/eps^2) samples; iterative amplitude') + print('estimation needs O(1/eps) oracle queries. The advantage is in the query') + print('count of the tail-probability primitive, not in the whole workflow, and') + print('it is asymptotic: on a simulator the classical computation is faster.') + + levels = distribution.probabilities.size + print(f'\nValue at Risk uses bisection over {levels} loss levels, so it needs') + print(f'O(log2({levels})) = {int(np.ceil(np.log2(levels)))} amplitude estimations') + print(f'rather than the {levels} a linear scan would take.') + + +def main(): + print_portfolio() + + distribution = LossDistribution.from_credit_portfolio( + DEFAULT_PROBABILITIES, EXPOSURES, unit=UNIT) + print_distribution(distribution) + + analyzer = QuantumRiskAnalyzer(distribution, epsilon=EPSILON) + print(f'\nquantum register: {analyzer.num_qubits} qubits ' + f'({distribution.num_qubits} loss + {distribution.num_qubits} comparator)') + + print_tail_probabilities(distribution, analyzer) + print_risk_measures(analyzer) + print_query_comparison(distribution) + + +if __name__ == '__main__': + main() diff --git a/contest2/OriginQCup_DU_Fanta/quantum_risk.py b/contest2/OriginQCup_DU_Fanta/quantum_risk.py new file mode 100644 index 00000000..e5a74219 --- /dev/null +++ b/contest2/OriginQCup_DU_Fanta/quantum_risk.py @@ -0,0 +1,479 @@ +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Quantum risk analysis: Value at Risk and Expected Shortfall by amplitude estimation. + +This module composes two existing ``pyqpanda_alg`` components into a complete +risk-measurement workflow: + +* ``QCmp.int_comparator`` builds the reversible predicate :math:`L \\geq t` ; +* ``QAE.IQAE`` estimates the probability that the predicate holds. + +Together they evaluate the tail probability :math:`P(L \\geq t)` with +:math:`O(1 / \\epsilon)` oracle queries, against the :math:`O(1 / \\epsilon ^ 2)` +samples a classical Monte Carlo simulation needs for the same accuracy. + +Every risk figure produced here reduces to tail probabilities, which is what +makes a single quantum primitive sufficient for the whole workflow. For an +integer-valued loss :math:`L` : + +.. math:: + + E[(L - t) ^ +] = \\sum_{k > t} P(L \\geq k) + + \\mathrm{ES}_\\alpha = \\mathrm{VaR}_\\alpha + + \\frac{E[(L - \\mathrm{VaR}_\\alpha) ^ +]}{P(L \\geq \\mathrm{VaR}_\\alpha)} + +so Value at Risk follows from a bisection over tail probabilities and Expected +Shortfall from a short sum of them. No separate estimation primitive is needed. + +References + [1] S. Woerner, D. J. Egger, Quantum risk analysis. npj Quantum Information 5, 15 (2019). + https://doi.org/10.1038/s41534-019-0130-6 + [2] D. Grinko, J. Gacon, C. Zoufal, S. Woerner, Iterative quantum amplitude estimation. + npj Quantum Information 7, 52 (2021). https://doi.org/10.1038/s41534-021-00379-1 +""" + +import warnings + +import numpy as np +from pyqpanda3.core import QCircuit, Encode + +from pyqpanda_alg import QAE +from pyqpanda_alg import QCmp + +__all__ = ['LossDistribution', 'QuantumRiskAnalyzer', 'recommended_epsilon', + 'classical_monte_carlo_samples'] + +# Value at Risk is decided by comparing tail probabilities against 1 - alpha, so +# the estimation accuracy has to be a small fraction of that gap or the +# comparison is dominated by estimation error. A fifth of the gap keeps the +# probability of choosing the wrong level low without making the estimate +# needlessly expensive. +EPSILON_SAFETY_FACTOR = 0.2 + + +def recommended_epsilon(confidence: float, safety: float = EPSILON_SAFETY_FACTOR) -> float: + """Accuracy needed to resolve Value at Risk at a given confidence level. + + The bisection in :meth:`QuantumRiskAnalyzer.value_at_risk` tests whether + :math:`P(L \\geq t + 1) \\leq 1 - \\alpha` . If the estimation accuracy + :math:`\\epsilon` is comparable to :math:`1 - \\alpha` , that comparison is + decided by noise and the returned level can be off by one. This helper + returns :math:`\\text{safety} \\times (1 - \\alpha)` . + + Parameters + confidence : ``float``\n + Confidence level :math:`\\alpha` in :math:`(0, 1)` . + safety : ``float``\n + Fraction of the gap :math:`1 - \\alpha` to use. Default 0.2. + + Returns + epsilon : ``float``\n + The recommended estimation accuracy. + + Examples + >>> from quantum_risk import recommended_epsilon + >>> round(recommended_epsilon(0.99), 4) + 0.002 + + """ + _check_confidence(confidence) + if not 0 < safety < 1: + raise ValueError(f'safety must lie in (0, 1), got {safety}') + return safety * (1.0 - confidence) + + +class LossDistribution: + """A discrete loss distribution held on :math:`n` qubits. + + The loss takes the integer values :math:`0, 1, \\ldots, 2 ^ n - 1` , which is + the representation the quantum comparator works on directly. Monetary losses + are recovered by scaling with :attr:`unit`. + + Parameters + probabilities : ``list[float]``, ``np.ndarray``\n + Probability of each loss level. The length is padded up to the next + power of two and the result is normalised. + unit : ``float``\n + Monetary value of one loss level. Default 1.0. + + Examples + >>> from quantum_risk import LossDistribution + >>> dist = LossDistribution([0.5, 0.3, 0.15, 0.05]) + >>> dist.num_qubits + 2 + >>> round(dist.tail_probability(2), 4) + 0.2 + + """ + + def __init__(self, probabilities, unit: float = 1.0): + probabilities = np.asarray(probabilities, dtype=float).flatten() + if probabilities.size == 0: + raise ValueError('probabilities must contain at least one entry') + if np.any(probabilities < 0): + raise ValueError('probabilities must be non-negative') + total = probabilities.sum() + if total <= 0: + raise ValueError('probabilities must contain at least one positive entry') + if unit <= 0: + raise ValueError(f'unit must be positive, got {unit}') + + size = 1 << max(int(np.ceil(np.log2(probabilities.size))), 1) + padded = np.zeros(size) + padded[:probabilities.size] = probabilities + self.probabilities = padded / padded.sum() + self.unit = float(unit) + + @property + def num_qubits(self) -> int: + """Number of qubits holding the loss value.""" + return int(np.log2(self.probabilities.size)) + + @property + def levels(self) -> np.ndarray: + """The integer loss levels :math:`0, \\ldots, 2 ^ n - 1` .""" + return np.arange(self.probabilities.size) + + @property + def amplitudes(self) -> np.ndarray: + """Amplitudes :math:`\\sqrt{p_i}` loaded by the state-preparation circuit.""" + return np.sqrt(self.probabilities) + + def state_preparation(self, qubits) -> QCircuit: + """Circuit mapping :math:`|0\\rangle ^ {\\otimes n}` to :math:`\\sum_i \\sqrt{p_i} |i\\rangle` . + + Parameters + qubits : ``list[int]``\n + The :math:`n` qubits holding the loss value. + + Returns + circuit : ``QCircuit``\n + The state-preparation circuit. + + """ + if len(qubits) != self.num_qubits: + raise ValueError( + f'expected {self.num_qubits} qubits for the loss register, got {len(qubits)}') + encoder = Encode() + encoder.amplitude_encode(list(qubits), list(self.amplitudes)) + circuit = QCircuit() + circuit << encoder.get_circuit() + return circuit + + def tail_probability(self, threshold: int) -> float: + """Exact :math:`P(L \\geq \\text{threshold})` , used as the reference value.""" + threshold = int(np.clip(threshold, 0, self.probabilities.size)) + return float(self.probabilities[threshold:].sum()) + + def value_at_risk(self, confidence: float) -> int: + """Exact Value at Risk: the smallest level :math:`t` with :math:`P(L \\leq t) \\geq \\alpha` .""" + _check_confidence(confidence) + cumulative = np.cumsum(self.probabilities) + return int(np.searchsorted(cumulative, confidence)) + + def expected_shortfall(self, confidence: float) -> float: + """Exact Expected Shortfall :math:`E[L \\mid L \\geq \\mathrm{VaR}_\\alpha]` .""" + var = self.value_at_risk(confidence) + tail = self.tail_probability(var) + if tail <= 0: + return float(var) + mask = self.levels >= var + return float((self.levels[mask] * self.probabilities[mask]).sum() / tail) + + @classmethod + def from_credit_portfolio(cls, default_probabilities, exposures, unit: float = 1.0): + """Build the loss distribution of a portfolio of independent obligors. + + Each obligor :math:`i` defaults independently with probability + :math:`p_i` and contributes an integer exposure :math:`e_i` when it + does, so the portfolio loss is :math:`L = \\sum_i e_i X_i` with + :math:`X_i \\sim \\mathrm{Bernoulli}(p_i)` . The exact distribution is + obtained by convolving the per-obligor distributions. + + Parameters + default_probabilities : ``list[float]``\n + Default probability of each obligor, each in :math:`[0, 1]` . + exposures : ``list[int]``\n + Integer loss contributed by each obligor on default. + unit : ``float``\n + Monetary value of one exposure unit. Default 1.0. + + Returns + distribution : ``LossDistribution``\n + The exact portfolio loss distribution. + + Examples + >>> from quantum_risk import LossDistribution + >>> dist = LossDistribution.from_credit_portfolio([0.1, 0.2], [1, 2]) + >>> round(dist.tail_probability(3), 4) + 0.02 + + """ + default_probabilities = np.asarray(default_probabilities, dtype=float).flatten() + exposures = np.asarray(exposures, dtype=int).flatten() + if default_probabilities.size != exposures.size: + raise ValueError('default_probabilities and exposures must have the same length') + if default_probabilities.size == 0: + raise ValueError('the portfolio must contain at least one obligor') + if np.any(default_probabilities < 0) or np.any(default_probabilities > 1): + raise ValueError('every default probability must lie in [0, 1]') + if np.any(exposures < 0): + raise ValueError('exposures must be non-negative') + + distribution = np.zeros(int(exposures.sum()) + 1) + distribution[0] = 1.0 + for probability, exposure in zip(default_probabilities, exposures): + updated = distribution * (1.0 - probability) + if exposure > 0: + updated[exposure:] += distribution[:distribution.size - exposure] * probability + else: + updated += distribution * probability + distribution = updated + return cls(distribution, unit=unit) + + +class QuantumRiskAnalyzer: + """Estimate tail risk measures with iterative quantum amplitude estimation. + + The circuit acting on the loss register and the comparator ancillas is + + .. parsed-literal:: + + loss register ---- state preparation ----*---- + | + comparator ------------------- int_comparator ----> objective qubit + + and ``QAE.IQAE`` estimates the probability that the objective qubit is + :math:`|1\\rangle` , which is exactly :math:`P(L \\geq t)` . + + Parameters + distribution : ``LossDistribution``\n + The loss distribution to analyse. + epsilon : ``float``\n + Target accuracy of each amplitude estimate. Default 0.005. + + Examples + >>> from quantum_risk import LossDistribution, QuantumRiskAnalyzer + >>> dist = LossDistribution([0.5, 0.3, 0.15, 0.05]) + >>> analyzer = QuantumRiskAnalyzer(dist, epsilon=0.01) + >>> estimate = analyzer.tail_probability(2) + >>> abs(estimate - dist.tail_probability(2)) < 0.05 + True + + """ + + def __init__(self, distribution: LossDistribution, epsilon: float = 5e-3): + if not isinstance(distribution, LossDistribution): + raise TypeError('distribution must be a LossDistribution instance') + if not 0 < epsilon < 1: + raise ValueError(f'epsilon must lie in (0, 1), got {epsilon}') + + self.distribution = distribution + self.epsilon = float(epsilon) + + n = distribution.num_qubits + # int_comparator needs one ancilla register the same width as the loss + # register, and reports the outcome on its last qubit. + self.loss_qubits = list(range(n)) + self.comparator_qubits = list(range(n, 2 * n)) + self.objective_qubit = self.comparator_qubits[-1] + self.num_qubits = 2 * n + self.estimate_calls = 0 + + def _state_operator(self, threshold: int): + """The operator :math:`A` whose :math:`|1\\rangle` amplitude on the objective qubit is the tail probability.""" + distribution = self.distribution + loss_qubits = self.loss_qubits + comparator_qubits = self.comparator_qubits + + def operator(qlist): + circuit = QCircuit() + circuit << distribution.state_preparation([qlist[i] for i in loss_qubits]) + circuit << QCmp.int_comparator(int(threshold), + [qlist[i] for i in loss_qubits], + [qlist[i] for i in comparator_qubits], + function='geq') + return circuit + + return operator + + def tail_probability(self, threshold: int) -> float: + """Estimate :math:`P(L \\geq \\text{threshold})` by amplitude estimation. + + Parameters + threshold : ``int``\n + The loss level to compare against. + + Returns + probability : ``float``\n + The estimated tail probability. + + """ + size = self.distribution.probabilities.size + # The comparator is only defined on representable levels; outside that + # range the answer is known exactly and no quantum work is needed. + if threshold <= 0: + return 1.0 + if threshold >= size: + return 0.0 + + self.estimate_calls += 1 + return float(QAE.IQAE(operator_in=self._state_operator(threshold), + qnumber=self.num_qubits, + epsilon=self.epsilon, + res_index=self.objective_qubit).run()) + + def value_at_risk(self, confidence: float) -> int: + """Estimate Value at Risk by bisection on quantum tail probabilities. + + :math:`\\mathrm{VaR}_\\alpha` is the smallest level :math:`t` with + :math:`P(L \\geq t + 1) \\leq 1 - \\alpha` . Because :math:`P(L \\geq t)` + is monotone in :math:`t` , bisection finds it in + :math:`O(\\log 2 ^ n) = O(n)` amplitude estimations rather than the + :math:`O(2 ^ n)` a linear scan would need. + + The accuracy :attr:`epsilon` must be small compared with :math:`1 - \\alpha` ; + see :func:`recommended_epsilon`. A warning is issued otherwise, because + the comparison would then be decided by estimation noise. + + Parameters + confidence : ``float``\n + Confidence level :math:`\\alpha` in :math:`(0, 1)` . + + Returns + var : ``int``\n + The estimated Value at Risk, as a loss level. + + """ + _check_confidence(confidence) + advised = recommended_epsilon(confidence) + if self.epsilon > advised: + warnings.warn( + f'epsilon={self.epsilon:g} is too coarse to resolve Value at Risk at ' + f'confidence {confidence:g}: the bisection compares tail probabilities ' + f'against {1.0 - confidence:g}. Use epsilon <= {advised:g}.', + RuntimeWarning, stacklevel=2) + low, high = 0, self.distribution.probabilities.size - 1 + target = 1.0 - confidence + while low < high: + middle = (low + high) // 2 + if self.tail_probability(middle + 1) <= target: + high = middle + else: + low = middle + 1 + return int(low) + + def expected_shortfall(self, confidence: float, var: int = None) -> float: + """Estimate Expected Shortfall :math:`E[L \\mid L \\geq \\mathrm{VaR}_\\alpha]` . + + Uses :math:`E[(L - t) ^ +] = \\sum_{k > t} P(L \\geq k)` , so the whole + quantity is a sum of tail probabilities and needs no second estimation + primitive. + + Parameters + confidence : ``float``\n + Confidence level :math:`\\alpha` in :math:`(0, 1)` . + var : ``int``, optional\n + A previously estimated Value at Risk. Recomputed when omitted. + + Returns + shortfall : ``float``\n + The estimated Expected Shortfall, as a loss level. + + """ + _check_confidence(confidence) + if var is None: + var = self.value_at_risk(confidence) + + tail = self.tail_probability(var) + if tail <= 0: + return float(var) + + size = self.distribution.probabilities.size + excess = sum(self.tail_probability(k) for k in range(var + 1, size)) + return float(var + excess / tail) + + def report(self, confidence: float) -> dict: + """Estimate every risk measure at one confidence level and compare with the exact values. + + Parameters + confidence : ``float``\n + Confidence level :math:`\\alpha` in :math:`(0, 1)` . + + Returns + summary : ``dict``\n + Quantum estimates, exact references, absolute errors, the number + of amplitude estimations used, and the monetary values obtained + by scaling with the distribution's unit. + + """ + _check_confidence(confidence) + self.estimate_calls = 0 + + quantum_var = self.value_at_risk(confidence) + quantum_tail = self.tail_probability(quantum_var) + quantum_shortfall = self.expected_shortfall(confidence, var=quantum_var) + + exact_var = self.distribution.value_at_risk(confidence) + exact_tail = self.distribution.tail_probability(exact_var) + exact_shortfall = self.distribution.expected_shortfall(confidence) + unit = self.distribution.unit + + return { + 'confidence': confidence, + 'epsilon': self.epsilon, + 'num_qubits': self.num_qubits, + 'amplitude_estimations': self.estimate_calls, + 'var_quantum': quantum_var, + 'var_exact': exact_var, + 'var_abs_error': abs(quantum_var - exact_var), + 'var_monetary': quantum_var * unit, + 'tail_probability_quantum': quantum_tail, + 'tail_probability_exact': exact_tail, + 'tail_probability_abs_error': abs(quantum_tail - exact_tail), + 'expected_shortfall_quantum': quantum_shortfall, + 'expected_shortfall_exact': exact_shortfall, + 'expected_shortfall_abs_error': abs(quantum_shortfall - exact_shortfall), + 'expected_shortfall_monetary': quantum_shortfall * unit, + } + + +def classical_monte_carlo_samples(epsilon: float, alpha: float = 0.05) -> int: + """Samples a classical Monte Carlo estimate needs for accuracy ``epsilon``. + + From the Chernoff-Hoeffding bound, estimating a probability to within + :math:`\\epsilon` with confidence :math:`1 - \\alpha` needs + :math:`\\lceil \\log(2 / \\alpha) / (2 \\epsilon ^ 2) \\rceil` samples. It is + quoted here so the :math:`O(1 / \\epsilon)` versus :math:`O(1 / \\epsilon ^ 2)` + scaling can be stated with a concrete number rather than asymptotically. + + Parameters + epsilon : ``float``\n + Target accuracy. + alpha : ``float``\n + Failure probability. Default 0.05. + + Returns + samples : ``int``\n + The required number of classical samples. + + """ + if not 0 < epsilon < 1: + raise ValueError(f'epsilon must lie in (0, 1), got {epsilon}') + return int(np.ceil(np.log(2.0 / alpha) / (2.0 * epsilon ** 2))) + + +def _check_confidence(confidence: float) -> None: + if not 0 < confidence < 1: + raise ValueError(f'confidence must lie in (0, 1), got {confidence}')