Bug fix for isolated nodes when running RQAOA

openqaoa/rqaoa/rqaoa.py

@@ -390,7 +390,14 @@ def redefine_problem(problem: QUBO, spin_map: dict):
     """
 
     # Define new QUBO problem as a dictionary
+    old_register = set([spin for term in problem.terms for spin in term])
     new_problem_dict = {}
+
+    # get a set of all the spins to be eliminated
+    eliminated_spins = set()
+    for spin in spin_map.keys():
+        if spin != spin_map[spin][1] or spin_map[spin][1] == None:
+            eliminated_spins.add(spin)
 
     # Scan all terms and weights
     for term, weight in zip(problem.terms, problem.weights):
@@ -478,10 +485,17 @@ def redefine_problem(problem: QUBO, spin_map: dict):
     # Define quadratic register after removing vanishing terms
     new_quadratic_register = set([spin for edge in new_problem_dict.keys() if len(edge) == 2 for spin in edge])
 
-    # If lengths do not match, there are isolated nodes
-    if len(new_register) != len(new_quadratic_register):
+    # Check for isolated nodes after constructing the new problem
+    there_is_isolated_nodes = False
+
+    if len(new_register) != len(new_quadratic_register):  # If lengths do not match, there are isolated nodes
         isolated_nodes = new_register.difference(new_quadratic_register)
-
+        there_is_isolated_nodes = True
+    elif old_register - new_register != eliminated_spins: # If too few eliminations, there are isolated nodes; only important for bias-free problems.
+        isolated_nodes = old_register.difference(new_register)
+        there_is_isolated_nodes = True
+
+    if there_is_isolated_nodes:
         # Fix isolated nodes
         for node in isolated_nodes:
             singlet = (node,)

openqaoa/workflows/optimizer.py

@@ -963,7 +963,7 @@ def optimize(self, verbose=False):
 
         # timestamp for the start of the optimization
         self.header['execution_time_start'] = int(time.time())
-
+        
         # If above cutoff, loop quantumly, else classically
         while n_qubits > n_cutoff:
 
@@ -981,7 +981,7 @@ def optimize(self, verbose=False):
             spin_map = rqaoa.spin_mapping(problem, max_terms_and_stats)
             # Eliminate spins and redefine problem
             new_problem, spin_map = rqaoa.redefine_problem(problem, spin_map)
-
+            
             # Extract final set of eliminations with correct dependencies and update tracker
             eliminations = [{'pair': (spin_map[spin][1],spin), 'correlation': spin_map[spin][0]} for spin in sorted(spin_map.keys()) if spin != spin_map[spin][1]]
             elimination_tracker.append(eliminations)
@@ -996,16 +996,16 @@ def optimize(self, verbose=False):
 
             # problem is updated
             problem = new_problem
-            
+
             # Compile qaoa with the problem
             q.compile(problem, verbose=False)
 
             # Add one step to the counter
             counter += 1
-
+        
         # Solve the new problem classically
         cl_energy, cl_ground_states = ground_state_hamiltonian(problem.hamiltonian)
-
+            
         # Retrieve full solutions including eliminated spins and their energies
         full_solutions = rqaoa.final_solution(
             elimination_tracker, cl_ground_states, self.problem.hamiltonian)

tests/test_rqaoa.py

@@ -14,9 +14,13 @@
 
 import numpy as np
 import unittest
+import networkx as nx
 
 from openqaoa.qaoa_parameters import Hamiltonian
 from openqaoa.rqaoa import *
+from openqaoa.problems.problem import MaximumCut
+from openqaoa.workflows.optimizer import RQAOA
+from openqaoa.devices import create_device
 
 """
 Unittest based testing of current implementation of the RQAOA Algorithm
@@ -289,6 +293,58 @@ def test_redefine_problem(self):
         assert np.allclose(hamiltonian.constant,comp_hamiltonian.constant), f'Constant in the computed Hamiltonian is incorrect'
 
 
+    def test_isolated_nodes_minimum_example(self):
+        """
+        Testing an edge case: solving MaxCut on a specific random unweighted graph must identify correctly the isolated nodes.
+
+        The test recreates the graph instance and MaxCut QUBO, runs standard RQAOA and compare the result to the expected one if the nodes has been correctly isolated.
+        """
+        test_qubo = QUBO(n=6, terms=[[0, 1], [2, 3], [2, 5], [3, 4], [4, 5]], weights=[2.0, 1.0, 1.0, 1.0, 1.0])
+        spin_map = {0: (1, 0), 1: (-1.0, 0), 2: (1, 2), 3: (1, 3), 4: (1, 4), 5: (1, 5)}
+        new_problem, new_spin_map = redefine_problem(test_qubo, spin_map)
+
+        assert new_problem.terms == [[0, 1], [0, 3], [1, 2], [2, 3]]
+
+
+    def test_isolated_nodes_whole_workflow(self):
+        """
+        Testing an edge case: solving MaxCut on a specific random unweighted graph must identify correctly the isolated nodes.
+
+        The test recreates the graph instance and MaxCut QUBO, runs standard RQAOA and compare the result to the expected one if the nodes has been correctly isolated.
+        """
+        # Generate the graph
+        g = nx.generators.gnp_random_graph(n=12, p=0.7, seed=7, directed=False)
+
+        # Define the problem and translate it into a binary QUBO.
+        maxcut_qubo = MaximumCut(g).get_qubo_problem()
+
+        # Define the RQAOA object
+        R = RQAOA()
+
+        # Set parameters for RQAOA: standard with cut off size 3 qubits
+        R.set_rqaoa_parameters(steps=1, n_cutoff=3)
+
+        # Set more parameters with a very specific starting point
+        R.set_circuit_properties(p=1, init_type='custom', variational_params_dict={"betas":[0.7044346364592513], "gammas":[8.926807699153894]}, mixer_hamiltonian='x')
+
+        # Define the device to be vectorized
+        device = create_device(location='local', name='vectorized')
+        R.set_device(device)
+
+        # Set the classical method used to optimiza over QAOA angles and its properties
+        R.set_classical_optimizer(method="cobyla", maxiter=200)
+
+        # Compile and optimize the problem instance on RQAOA
+        R.compile(maxcut_qubo)
+        R.optimize()
+
+        # Get results
+        opt_results = R.results
+
+        # Compare results to known behaviour:
+        assert opt_results['schedule'] == [1, 1, 1, 1, 1, 1, 2, 1]
+        assert opt_results['solution'] == {'011001101001': -6.0, '011000110011': -6.0}
+
 if __name__ == "__main__":
 	unittest.main()