Fix GroverOptimizer's rotation_count (backport #132)

qiskit_optimization/algorithms/grover_optimizer.py

@@ -128,7 +128,7 @@ def _get_oracle(self, qr_key_value):
         oracle = QuantumCircuit(qr_key_value, oracle_bit)
         oracle.z(self._num_key_qubits)  # recognize negative values.
 
-        def is_good_state(self, measurement):
+        def is_good_state(measurement):
             """Check whether ``measurement`` is a good state or not."""
             value = measurement[self._num_key_qubits:self._num_key_qubits + self._num_value_qubits]
             return value[0] == '1'
@@ -209,7 +209,7 @@ def solve(self, problem: QuadraticProgram) -> OptimizationResult:
             while not improvement_found:
                 # Determine the number of rotations.
                 loops_with_no_improvement += 1
-                rotation_count = int(np.ceil(algorithm_globals.random.uniform(0, m - 1)))
+                rotation_count = algorithm_globals.random.integers(0, m)
                 rotations += rotation_count
                 # Apply Grover's Algorithm to find values below the threshold.
                 # TODO: Utilize Grover's incremental feature - requires changes to Grover.
@@ -290,11 +290,9 @@ def solve(self, problem: QuadraticProgram) -> OptimizationResult:
     def _measure(self, circuit: QuantumCircuit) -> str:
         """Get probabilities from the given backend, and picks a random outcome."""
         probs = self._get_probs(circuit)
-        freq = sorted(probs.items(), key=lambda x: x[1], reverse=True)
+        logger.info("Frequencies: %s", probs)
         # Pick a random outcome.
-        idx = algorithm_globals.random.choice(len(freq), 1, p=[x[1] for x in freq])[0]
-        logger.info('Frequencies: %s', freq)
-        return freq[idx][0]
+        return algorithm_globals.random.choice(list(probs.keys()), 1, p=list(probs.values()))[0]
 
     def _get_probs(self, qc: QuantumCircuit) -> Dict[str, float]:
         """Gets probabilities from a given backend."""
@@ -312,8 +310,8 @@ def _get_probs(self, qc: QuantumCircuit) -> Dict[str, float]:
         else:
             state = result.get_counts(qc)
             shots = self.quantum_instance.run_config.shots
-            hist = {key[::-1]: val / shots for key, val in state.items() if val > 0}
-            self._circuit_results = {b[::-1]: np.sqrt(v / shots) for (b, v) in state.items()}
+            hist = {key[::-1]: val / shots for key, val in sorted(state.items()) if val > 0}
+            self._circuit_results = {b: (v / shots) ** 0.5 for (b, v) in state.items()}
         return hist
 
     @staticmethod

releasenotes/notes/fix-grover-rotation-count-853baefc6b7e8476.yaml

@@ -0,0 +1,12 @@
+---
+fixes:
+  - |
+    Fixes ``rotation_count`` in :class`qiskit_optimization.algorithms.GroverOptimizer`.
+    This fix uses ``algorithm_globals.random.integers(0, m)`` to generate a random integer
+    in a range 0..``m``-1.
+  - |
+    Sorts the order of ``result.get_counts(qc)`` by bitstring
+    in :class`qiskit_optimization.algorithms.GroverOptimizer` when ``qasm_simulator`` is used
+    so that the algorithm behaves deterministically.
+    The previous version sorts the counts by probabilities, but some bitstrings may have
+    the same probability and the algorithm could behave probabilistically.

test/algorithms/test_grover_optimizer.py

@@ -21,13 +21,17 @@
 from qiskit import Aer
 from qiskit.utils import QuantumInstance, algorithm_globals
 from qiskit.algorithms import NumPyMinimumEigensolver
-from qiskit_optimization.algorithms import (GroverOptimizer,
-                                            MinimumEigenOptimizer,
-                                            OptimizationResultStatus)
-from qiskit_optimization.converters import (InequalityToEquality,
-                                            IntegerToBinary,
-                                            LinearEqualityToPenalty,
-                                            QuadraticProgramToQubo)
+from qiskit_optimization.algorithms import (
+    GroverOptimizer,
+    MinimumEigenOptimizer,
+    OptimizationResultStatus,
+)
+from qiskit_optimization.converters import (
+    InequalityToEquality,
+    IntegerToBinary,
+    LinearEqualityToPenalty,
+    QuadraticProgramToQubo,
+)
 from qiskit_optimization.problems import QuadraticProgram
 
 
@@ -38,17 +42,21 @@ class TestGroverOptimizer(QiskitOptimizationTestCase):
     def setUp(self):
         super().setUp()
         algorithm_globals.random_seed = 1
-        self.sv_simulator = QuantumInstance(Aer.get_backend('statevector_simulator'),
-                                            seed_simulator=921, seed_transpiler=200)
-        self.qasm_simulator = QuantumInstance(Aer.get_backend('qasm_simulator'),
-                                              seed_simulator=123, seed_transpiler=123)
+        self.sv_simulator = QuantumInstance(
+            Aer.get_backend("statevector_simulator"),
+            seed_simulator=921,
+            seed_transpiler=200,
+        )
+        self.qasm_simulator = QuantumInstance(
+            Aer.get_backend("qasm_simulator"), seed_simulator=123, seed_transpiler=123
+        )
+        self.n_iter = 8
 
     def validate_results(self, problem, results):
         """Validate the results object returned by GroverOptimizer."""
         # Get expected value.
         solver = MinimumEigenOptimizer(NumPyMinimumEigensolver())
         comp_result = solver.solve(problem)
-
         # Validate results.
         np.testing.assert_array_almost_equal(comp_result.x, results.x)
         self.assertEqual(comp_result.fval, results.fval)
@@ -59,9 +67,9 @@ def test_qubo_gas_int_zero(self):
 
         # Input.
         model = Model()
-        x_0 = model.binary_var(name='x0')
-        x_1 = model.binary_var(name='x1')
-        model.minimize(0*x_0+0*x_1)
+        x_0 = model.binary_var(name="x0")
+        x_1 = model.binary_var(name="x1")
+        model.minimize(0 * x_0 + 0 * x_1)
         op = QuadraticProgram()
         op.from_docplex(model)
 
@@ -77,15 +85,14 @@ def test_qubo_gas_int_simple(self):
 
         # Input.
         model = Model()
-        x_0 = model.binary_var(name='x0')
-        x_1 = model.binary_var(name='x1')
-        model.minimize(-x_0+2*x_1)
+        x_0 = model.binary_var(name="x0")
+        x_1 = model.binary_var(name="x1")
+        model.minimize(-x_0 + 2 * x_1)
         op = QuadraticProgram()
         op.from_docplex(model)
 
         # Get the optimum key and value.
-        n_iter = 8
-        gmf = GroverOptimizer(4, num_iterations=n_iter, quantum_instance=self.sv_simulator)
+        gmf = GroverOptimizer(4, num_iterations=self.n_iter, quantum_instance=self.sv_simulator)
         results = gmf.solve(op)
         self.validate_results(op, results)
 
@@ -98,19 +105,18 @@ def test_qubo_gas_int_simple_maximize(self):
 
         # Input.
         model = Model()
-        x_0 = model.binary_var(name='x0')
-        x_1 = model.binary_var(name='x1')
-        model.maximize(-x_0+2*x_1)
+        x_0 = model.binary_var(name="x0")
+        x_1 = model.binary_var(name="x1")
+        model.maximize(-x_0 + 2 * x_1)
         op = QuadraticProgram()
         op.from_docplex(model)
 
         # Get the optimum key and value.
-        n_iter = 8
-        gmf = GroverOptimizer(4, num_iterations=n_iter, quantum_instance=self.sv_simulator)
+        gmf = GroverOptimizer(4, num_iterations=self.n_iter, quantum_instance=self.sv_simulator)
         results = gmf.solve(op)
         self.validate_results(op, results)
 
-    @data('sv', 'qasm')
+    @data("sv", "qasm")
     def test_qubo_gas_int_paper_example(self, simulator):
         """
         Test the example from https://arxiv.org/abs/1912.04088 using the state vector simulator
@@ -119,81 +125,127 @@ def test_qubo_gas_int_paper_example(self, simulator):
 
         # Input.
         model = Model()
-        x_0 = model.binary_var(name='x0')
-        x_1 = model.binary_var(name='x1')
-        x_2 = model.binary_var(name='x2')
-        model.minimize(-x_0+2*x_1-3*x_2-2*x_0*x_2-1*x_1*x_2)
+        x_0 = model.binary_var(name="x0")
+        x_1 = model.binary_var(name="x1")
+        x_2 = model.binary_var(name="x2")
+        model.minimize(-x_0 + 2 * x_1 - 3 * x_2 - 2 * x_0 * x_2 - 1 * x_1 * x_2)
         op = QuadraticProgram()
         op.from_docplex(model)
 
         # Get the optimum key and value.
-        n_iter = 10
-
-        q_instance = self.sv_simulator if simulator == 'sv' else self.qasm_simulator
-        gmf = GroverOptimizer(6, num_iterations=n_iter, quantum_instance=q_instance)
+        q_instance = self.sv_simulator if simulator == "sv" else self.qasm_simulator
+        gmf = GroverOptimizer(6, num_iterations=self.n_iter, quantum_instance=q_instance)
         results = gmf.solve(op)
         self.validate_results(op, results)
 
     def test_converter_list(self):
         """Test converters list"""
         # Input.
+
         model = Model()
-        x_0 = model.binary_var(name='x0')
-        x_1 = model.binary_var(name='x1')
-        model.maximize(-x_0+2*x_1)
+        x_0 = model.binary_var(name="x0")
+        x_1 = model.binary_var(name="x1")
+        model.maximize(-x_0 + 2 * x_1)
         op = QuadraticProgram()
         op.from_docplex(model)
 
         # Get the optimum key and value.
-        n_iter = 8
         # a single converter.
         qp2qubo = QuadraticProgramToQubo()
-        gmf = GroverOptimizer(4, num_iterations=n_iter, quantum_instance=self.sv_simulator,
-                              converters=qp2qubo)
+        gmf = GroverOptimizer(
+            4,
+            num_iterations=self.n_iter,
+            quantum_instance=self.sv_simulator,
+            converters=qp2qubo,
+        )
         results = gmf.solve(op)
         self.validate_results(op, results)
+
         # a list of converters
         ineq2eq = InequalityToEquality()
         int2bin = IntegerToBinary()
         penalize = LinearEqualityToPenalty()
         converters = [ineq2eq, int2bin, penalize]
-        gmf = GroverOptimizer(4, num_iterations=n_iter, quantum_instance=self.sv_simulator,
-                              converters=converters)
+        gmf = GroverOptimizer(
+            4,
+            num_iterations=self.n_iter,
+            quantum_instance=self.sv_simulator,
+            converters=converters,
+        )
         results = gmf.solve(op)
         self.validate_results(op, results)
         # invalid converters
         with self.assertRaises(TypeError):
             invalid = [qp2qubo, "invalid converter"]
-            GroverOptimizer(4, num_iterations=n_iter,
-                            quantum_instance=self.sv_simulator,
-                            converters=invalid)
-
-    def test_samples_and_raw_samples(self):
+            GroverOptimizer(
+                4,
+                num_iterations=self.n_iter,
+                quantum_instance=self.sv_simulator,
+                converters=invalid,
+            )
+
+    @data("sv", "qasm")
+    def test_samples_and_raw_samples(self, simulator):
         """Test samples and raw_samples"""
+        algorithm_globals.random_seed = 2
         op = QuadraticProgram()
-        op.integer_var(0, 3, 'x')
-        op.binary_var('y')
-        op.minimize(linear={'x': 1, 'y': 2})
-        op.linear_constraint(linear={'x': 1, 'y': 1}, sense='>=', rhs=1, name='xy')
+        op.integer_var(0, 3, "x")
+        op.binary_var("y")
+        op.minimize(linear={"x": 1, "y": 2})
+        op.linear_constraint(linear={"x": 1, "y": 1}, sense=">=", rhs=1, name="xy")
+        q_instance = self.sv_simulator if simulator == "sv" else self.qasm_simulator
+        grover_optimizer = GroverOptimizer(
+            8, num_iterations=self.n_iter, quantum_instance=q_instance
+        )
         opt_sol = 1
         success = OptimizationResultStatus.SUCCESS
-        algorithm_globals.random_seed = 1
+        results = grover_optimizer.solve(op)
+        self.assertEqual(len(results.samples), 8)
+        self.assertEqual(len(results.raw_samples), 32)
+        self.assertAlmostEqual(sum(s.probability for s in results.samples), 1)
+        self.assertAlmostEqual(sum(s.probability for s in results.raw_samples), 1)
+        self.assertAlmostEqual(min(s.fval for s in results.samples), 0)
+        self.assertAlmostEqual(min(s.fval for s in results.samples if s.status == success), opt_sol)
+        self.assertAlmostEqual(min(s.fval for s in results.raw_samples), opt_sol)
+        for sample in results.raw_samples:
+            self.assertEqual(sample.status, success)
+        np.testing.assert_array_almost_equal(results.x, results.samples[0].x)
+        self.assertAlmostEqual(results.fval, results.samples[0].fval)
+        self.assertEqual(results.status, results.samples[0].status)
+        self.assertAlmostEqual(results.fval, results.raw_samples[0].fval)
+        self.assertEqual(results.status, results.raw_samples[0].status)
+        np.testing.assert_array_almost_equal([1, 0, 0, 0, 0], results.raw_samples[0].x)
+
+    @data("sv", "qasm")
+    def test_bit_ordering(self, simulator):
+        """Test bit ordering"""
+        # test minimize
+        algorithm_globals.random_seed = 2
+        q_instance = self.sv_simulator if simulator == "sv" else self.qasm_simulator
+        mdl = Model("docplex model")
+        x = mdl.binary_var("x")
+        y = mdl.binary_var("y")
+        mdl.minimize(x - 2 * y)
+        op = QuadraticProgram()
+        op.from_docplex(mdl)
+        opt_sol = -2
+        success = OptimizationResultStatus.SUCCESS
         grover_optimizer = GroverOptimizer(
-            8, num_iterations=5, quantum_instance=self.qasm_simulator)
-        result = grover_optimizer.solve(op)
-        self.assertEqual(len(result.samples), 8)
-        self.assertEqual(len(result.raw_samples), 32)
-        self.assertAlmostEqual(sum(s.probability for s in result.samples), 1)
-        self.assertAlmostEqual(sum(s.probability for s in result.raw_samples), 1)
-        self.assertAlmostEqual(min(s.fval for s in result.samples), 0)
-        self.assertAlmostEqual(min(s.fval for s in result.samples if s.status == success), opt_sol)
-        self.assertAlmostEqual(min(s.fval for s in result.raw_samples), opt_sol)
-        for sample in result.raw_samples:
+            3, num_iterations=self.n_iter, quantum_instance=q_instance
+        )
+        results = grover_optimizer.solve(op)
+        self.assertEqual(results.fval, opt_sol)
+        np.testing.assert_array_almost_equal(results.x, [0, 1])
+        self.assertEqual(results.status, success)
+        results.raw_samples.sort(key=lambda x: x.probability, reverse=True)
+        self.assertAlmostEqual(sum(s.probability for s in results.samples), 1, delta=1e-5)
+        self.assertAlmostEqual(sum(s.probability for s in results.raw_samples), 1, delta=1e-5)
+        self.assertAlmostEqual(min(s.fval for s in results.samples), -2)
+        self.assertAlmostEqual(min(s.fval for s in results.samples if s.status == success), opt_sol)
+        self.assertAlmostEqual(min(s.fval for s in results.raw_samples), opt_sol)
+        for sample in results.raw_samples:
             self.assertEqual(sample.status, success)
-        np.testing.assert_array_almost_equal(result.x, result.raw_samples[0].x[0:2])
-        self.assertAlmostEqual(result.fval, result.raw_samples[0].fval)
-        self.assertEqual(result.status, result.raw_samples[0].status)
 
 
-if __name__ == '__main__':
+if __name__ == "__main__":
     unittest.main()