Fix recursive optimization (qiskit-community/qiskit-aqua#976)

* fix recursive optimization

* Update recursive_minimum_eigen_optimizer.py

Co-authored-by: Manoel Marques <[email protected]>
Co-authored-by: Julien Gacon <[email protected]>

Author: 3 people

qiskit/optimization/algorithms/minimum_eigen_optimizer.py

@@ -152,17 +152,29 @@ def solve(self, problem: QuadraticProgram) -> MinimumEigenOptimizerResult:
         operator_converter = QuadraticProgramToIsing()
         operator, offset = operator_converter.encode(problem_)
 
-        # approximate ground state of operator using min eigen solver
-        eigen_results = self._min_eigen_solver.compute_minimum_eigenvalue(operator)
-
-        # analyze results
-        samples = eigenvector_to_solutions(eigen_results.eigenstate, operator)
-        samples = [(res[0], problem_.objective.sense.value * (res[1] + offset), res[2])
-                   for res in samples]
-        samples.sort(key=lambda x: problem_.objective.sense.value * x[1])
+        # only try to solve non-empty Ising Hamiltonians
+        if operator.num_qubits > 0:
+
+            # approximate ground state of operator using min eigen solver
+            eigen_results = self._min_eigen_solver.compute_minimum_eigenvalue(operator)
+
+            # analyze results
+            samples = eigenvector_to_solutions(eigen_results.eigenstate, operator)
+            samples = [(res[0], problem_.objective.sense.value * (res[1] + offset), res[2])
+                       for res in samples]
+            samples.sort(key=lambda x: problem_.objective.sense.value * x[1])
+            x = samples[0][0]
+            fval = samples[0][1]
+
+        # if Hamiltonian is empty, then the objective function is constant to the offset
+        else:
+            x = [0]*problem_.get_num_binary_vars()
+            fval = offset
+            x_str = '0'*problem_.get_num_binary_vars()
+            samples = [(x_str, offset, 1.0)]
 
         # translate result back to integers
-        opt_res = MinimumEigenOptimizerResult(samples[0][0], samples[0][1], samples, qubo_converter)
+        opt_res = MinimumEigenOptimizerResult(x, fval, samples, qubo_converter)
         opt_res = qubo_converter.decode(opt_res)
 
         # translate results back to original problem

qiskit/optimization/algorithms/recursive_minimum_eigen_optimizer.py

@@ -146,7 +146,6 @@ def solve(self, problem: QuadraticProgram) -> OptimizationResult:
             x_j = problem_.variables[j].name
             if correlations[i, j] > 0:
                 # set x_i = x_j
-                problem_.substitute_variables()
                 problem_ = problem_.substitute_variables(variables={i: (j, 1)})
                 if problem_.status == QuadraticProgram.Status.INFEASIBLE:
                     raise QiskitOptimizationError('Infeasible due to variable substitution')
@@ -158,16 +157,23 @@ def solve(self, problem: QuadraticProgram) -> OptimizationResult:
 
                 # 1a. get additional offset
                 constant = problem_.objective.constant
-                constant += problem_.objective.quadratic[i, i]
                 constant += problem_.objective.linear[i]
+                constant += problem_.objective.quadratic[i, i]
                 problem_.objective.constant = constant
 
                 # 1b. get additional linear part
                 for k in range(problem_.get_num_vars()):
-                    coeff = problem_.objective.quadratic[i, k]
+                    coeff = problem_.objective.linear[k]
+                    if k == i:
+                        coeff += 2*problem_.objective.quadratic[i, k]
+                    else:
+                        coeff += problem_.objective.quadratic[i, k]
+
+                    # set new coefficient if not too small
                     if np.abs(coeff) > 1e-10:
-                        coeff += problem_.objective.linear[k]
                         problem_.objective.linear[k] = coeff
+                    else:
+                        problem_.objective.linear[k] = 0
 
                 # 2. replace x_i by -x_j
                 problem_ = problem_.substitute_variables(variables={i: (j, -1)})
@@ -211,7 +217,16 @@ def find_value(x, replacements, var_values):
         return results
 
     def _find_strongest_correlation(self, correlations):
-        m_max = np.argmax(np.abs(correlations.flatten()))
+
+        # get absolute values and set diagonal to -1 to make sure maximum is always on off-diagonal
+        abs_correlations = np.abs(correlations)
+        for i in range(len(correlations)):
+            abs_correlations[i, i] = -1
+
+        # get index of maximum (by construction on off-diagonal)
+        m_max = np.argmax(abs_correlations.flatten())
+
+        # translate back to indices
         i = int(m_max // len(correlations))
         j = int(m_max - i*len(correlations))
         return (i, j)