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}')