【代码贡献】修复并加速 QAOA 的 CVaR/Gibbs 风险目标函数

pyqpanda-algorithm/pyqpanda_alg/QAOA/qaoa.py

@@ -10,11 +10,14 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
+from numbers import Real
+
 from pyqpanda3.core import CPUQVM, QCircuit, QProg, I, H, RZ, RX, CNOT, measure
 from pyqpanda3.hamiltonian import PauliOperator, Hamiltonian
 import numpy as np
 from scipy.optimize import minimize
 from scipy.interpolate import barycentric_interpolate as b_interp
+from scipy.special import logsumexp
 import sympy as sp
 from . import spsa
 from .default_circuits import *
@@ -309,9 +312,13 @@ def __init__(self, problem, init_circuit=None,
         problem_dimension = 0
         self.problem = problem
         self.operator = None
+        self._problem_symbols = ()
+        self._problem_function = None
         if isinstance(problem, sp.Basic):
             self.problem = sp.simplify(problem)
-            problem_dimension = len(problem.atoms(sp.Symbol))
+            self._problem_symbols = tuple(sorted(self.problem.free_symbols, key=lambda symbol: symbol.name))
+            self._problem_function = sp.lambdify(self._problem_symbols, self.problem, 'numpy')
+            problem_dimension = len(self._problem_symbols)
             self.operator = problem_to_z_operator(self.problem, norm)
 
         elif isinstance(problem, Hamiltonian):
@@ -347,7 +354,6 @@ def __init__(self, problem, init_circuit=None,
     def calculate_energy(self, x):
         """
         Calculate the function value for one solution.
-        TODO: using new method to acccelrate the calculation.
 
         Parameter
             x : ``array-like``\n
@@ -382,13 +388,7 @@ def calculate_energy(self, x):
         bit_form = x
         result = 0
         if isinstance(self.problem, sp.Basic):
-            symbols = sorted(self.problem.free_symbols, key=lambda symbol: symbol.name)
-            # problem = sp.Poly(self.problem)
-            value_dict = {}
-            for i in range(len(x)):
-                value_dict[symbols[i]] = bit_form[i]
-            f = sp.lambdify(symbols, self.problem, 'numpy')
-            raw_result = f(*bit_form)
+            raw_result = self._problem_function(*bit_form)
             result = raw_result.real if isinstance(raw_result, complex) else raw_result
 
         if isinstance(self.problem, PauliOperator):
@@ -559,7 +559,7 @@ def _loss_function_default(self, measure_result):
 
     def _loss_function_cvar(self, measure_result):
         """
-        Given a result, calculate the CVaR energy expectation.
+        Given a result, calculate the normalized lower-tail CVaR energy expectation.
 
         Parameter
             measure_result : ``dict``\n
@@ -569,26 +569,27 @@ def _loss_function_cvar(self, measure_result):
             lost : ``float``\n
                 CVaR energy expectation
         """
-        if not any(isinstance(self.alpha, t) for t in [int, float]):
-            raise ValueError('CVaR method needs parameter alpha to be a number between 0~1')
-        if self.alpha > 1 or self.alpha < 0:
-            raise ValueError('CVaR method needs parameter alpha to be a number between 0~1')
+        if (isinstance(self.alpha, bool) or not isinstance(self.alpha, Real)
+                or not np.isfinite(self.alpha) or not 0 < self.alpha <= 1):
+            raise ValueError('CVaR method needs parameter alpha to be a number in (0, 1]')
         cdf = 0.
         loss = 0.
 
-        measure_result = sorted(measure_result.items(), key=lambda k: k[1], reverse=True)
-        for solution, hits in measure_result:
+        energy_distribution = []
+        for solution, hits in measure_result.items():
             if solution not in self.energy_dict:
                 solution_list = [int(i) for i in solution[::-1]]
                 self.energy_dict[solution] = self.calculate_energy(solution_list)
-            prob = hits
-            if cdf < self.alpha:
-                if cdf + prob < self.alpha:
-                    loss += self.energy_dict[solution] * prob
-                else:
-                    loss += self.energy_dict[solution] * (self.alpha - cdf)
-                cdf += prob
-        return loss
+            energy_distribution.append((self.energy_dict[solution], hits))
+
+        for energy, prob in sorted(energy_distribution, key=lambda item: item[0]):
+            if cdf >= self.alpha:
+                break
+            included_probability = min(prob, self.alpha - cdf)
+            loss += energy * included_probability
+            cdf += included_probability
+
+        return loss / self.alpha
 
     def _loss_function_Gibbs(self, measure_result):
         """
@@ -602,17 +603,19 @@ def _loss_function_Gibbs(self, measure_result):
             lost : ``float``\n
                 Gibbs energy expectation
         """
-        if not any(isinstance(self.temperature, t) for t in [int, float]):
-            raise ValueError('Gibbs free energy method needs parameter temperature to be a number between 0~1')
-        if self.temperature > 1 or self.temperature < 0:
-            raise ValueError('Gibbs free energy method needs parameter temperature to be a number between 0~1')
-        lost = 0.
+        if (isinstance(self.temperature, bool) or not isinstance(self.temperature, Real)
+                or not np.isfinite(self.temperature) or not 0 < self.temperature <= 1):
+            raise ValueError('Gibbs free energy method needs parameter temperature to be a number in (0, 1]')
+        log_terms = []
         for solution, hits in measure_result.items():
             if solution not in self.energy_dict:
                 solution_list = [int(i) for i in solution[::-1]]
                 self.energy_dict[solution] = self.calculate_energy(solution_list)
-            lost += hits * np.exp(-self.energy_dict[solution] / self.temperature)
-        return - np.log(lost)
+            if hits > 0:
+                log_terms.append(np.log(hits) - self.energy_dict[solution] / self.temperature)
+        if not log_terms:
+            raise ValueError('Gibbs free energy method needs at least one positive probability')
+        return -float(logsumexp(log_terms))
 
     def _loss_function(self, paras):
         """
@@ -844,10 +847,12 @@ def run(self, layer=1, initial_para=None, shots=-1, loss_type=None, optimize_typ
             loss_option :\n
 
                 temperature : ``float``, ``optional``\n
-                    parameter calculated in _loss_function_Gibbs. Default is 1. See Note ``Gibbs energy``.
+                    Parameter in :math:`(0, 1]` calculated in _loss_function_Gibbs. Default is 1.
+                    See Note ``Gibbs energy``.
 
                 alpha : ``float``, ``optional``\n
-                    parameter calculated in _loss_function_cvar. Default is 1. See Note ``Gibbs energy``.
+                    Confidence level in :math:`(0, 1]` calculated in _loss_function_cvar. Default is 1.
+                    See Note ``CVaR loss function``.
 
         Return
             qaoa_result : ``dict``\n
@@ -926,16 +931,19 @@ def run(self, layer=1, initial_para=None, shots=-1, loss_type=None, optimize_typ
                     :math:`CVaR_\\alpha(X) = \mathbb{E}[X|X\leq F_X^{-1}(alpha)]`
 
                 Here :math:`\\alpha` is the confidence level. CVaR is the expected value of the lower α-tail of the
-                distribution of X. :math:`\\alpha=0` corresponds to the minimum, and :math:`\\alpha=1` corresponds to the
-                expectation value.
+                distribution of X. As :math:`\\alpha` approaches zero it approaches the minimum, while
+                :math:`\\alpha=1` corresponds to the expectation value.
 
                 If measure type is sample, it is calculated by
 
-                    :math:`E=\\frac{1}{\\alpha N}(\sum_{i=0}^{k} n_iE_i + (\\alpha N - n_{k+1})E_{k+1}),\sum_{i=0}^k n_i < \\alpha N`
+                    :math:`E=\\frac{1}{\\alpha N}(\sum_{i=0}^{k-1} n_iE_i + (\\alpha N - \sum_{i=0}^{k-1}n_i)E_k)`.
 
                 If measure type is theoretical, it is calculated by
 
-                    :math:`E=\sum_{i=0}^{k} p_iE_i + (\\alpha - p_{k+1})E_{k+1}, \sum_{i=0}^k p_i < \\alpha`
+                    :math:`E=\\frac{1}{\\alpha}(\sum_{i=0}^{k-1} p_iE_i + (\\alpha - \sum_{i=0}^{k-1}p_i)E_k)`.
+
+                In both formulas the states are ordered by ascending energy and :math:`k` is the state where the
+                cumulative probability first reaches :math:`\\alpha`.
 
             - Interpolate method:\n
                 Inspired by Ref[2].

test/QAOA/Test_qaoa_QAOA_risk_objectives.py

@@ -0,0 +1,92 @@
+import numpy as np
+import pytest
+import sympy as sp
+
+from pyqpanda_alg.QAOA.qaoa import QAOA
+
+
+@pytest.fixture
+def risk_qaoa():
+    qaoa = QAOA.__new__(QAOA)
+    qaoa.energy_dict = {
+        '00': 3.0,
+        '01': -1.0,
+        '10': -4.0,
+    }
+    return qaoa
+
+
+def test_cvar_uses_lowest_energy_tail_and_normalizes(risk_qaoa):
+    risk_qaoa.alpha = 0.2
+    distribution = {'00': 0.7, '01': 0.2, '10': 0.1}
+
+    result = risk_qaoa._loss_function_cvar(distribution)
+
+    # The lowest 20% consists of all p=0.1 at E=-4 and p=0.1 at E=-1.
+    assert result == pytest.approx((-4.0 * 0.1 - 1.0 * 0.1) / 0.2)
+
+
+def test_cvar_alpha_one_matches_energy_expectation(risk_qaoa):
+    risk_qaoa.alpha = 1.0
+    distribution = {'00': 0.7, '01': 0.2, '10': 0.1}
+
+    assert risk_qaoa._loss_function_cvar(distribution) == pytest.approx(
+        risk_qaoa._loss_function_default(distribution)
+    )
+
+
+@pytest.mark.parametrize('alpha', [0, -0.1, 1.1, np.nan, np.inf, True, '0.5'])
+def test_cvar_rejects_invalid_alpha(risk_qaoa, alpha):
+    risk_qaoa.alpha = alpha
+
+    with pytest.raises(ValueError, match=r'\(0, 1\]'):
+        risk_qaoa._loss_function_cvar({'00': 1.0})
+
+
+def test_gibbs_matches_direct_formula_for_regular_values(risk_qaoa):
+    risk_qaoa.temperature = 0.5
+    distribution = {'00': 0.7, '01': 0.2, '10': 0.1}
+    expected = -np.log(sum(
+        probability * np.exp(-risk_qaoa.energy_dict[solution] / risk_qaoa.temperature)
+        for solution, probability in distribution.items()
+    ))
+
+    assert risk_qaoa._loss_function_Gibbs(distribution) == pytest.approx(expected)
+
+
+def test_gibbs_remains_finite_for_large_energy_magnitudes(risk_qaoa):
+    risk_qaoa.temperature = 0.1
+    risk_qaoa.energy_dict = {'0': -1000.0, '1': 1000.0}
+
+    result = risk_qaoa._loss_function_Gibbs({'0': 0.5, '1': 0.5})
+
+    assert np.isfinite(result)
+    assert result == pytest.approx(-10000.0 + np.log(2.0))
+
+
+@pytest.mark.parametrize('temperature', [0, -0.1, 1.1, np.nan, np.inf, True, '0.5'])
+def test_gibbs_rejects_invalid_temperature(risk_qaoa, temperature):
+    risk_qaoa.temperature = temperature
+
+    with pytest.raises(ValueError, match=r'\(0, 1\]'):
+        risk_qaoa._loss_function_Gibbs({'00': 1.0})
+
+
+def test_gibbs_rejects_distribution_without_positive_probability(risk_qaoa):
+    risk_qaoa.temperature = 0.5
+
+    with pytest.raises(ValueError, match='positive probability'):
+        risk_qaoa._loss_function_Gibbs({'00': 0.0})
+
+
+def test_symbolic_energy_function_is_compiled_once(monkeypatch):
+    x0, x1 = sp.symbols('x0:2')
+    qaoa = QAOA(2 * x0 * x1 - x0)
+
+    def fail_if_recompiled(*args, **kwargs):
+        raise AssertionError('calculate_energy should reuse the compiled function')
+
+    monkeypatch.setattr(sp, 'lambdify', fail_if_recompiled)
+
+    assert qaoa.calculate_energy([1, 0]) == -1
+    assert qaoa.calculate_energy([1, 1]) == 1