Qiskit: Taking a QUBO matrix into `qubit_op'

Question

I'm trying to solve the maximum independent set problem using Qiskit and the QAOA.

I've a nice QUBO Matrix for this simple path graph:

as so:

My question is, how do I convert this into a general DOCPLEX method which can be put into the QuadracticProgram method of Qiskit? How can this be done for a general QUBO matrix?

I realise that for this particular problem I've found the matrix $Q$ as described in this [guide].(https://arxiv.org/ftp/arxiv/papers/1811/1811.11538.pdf)

So by doing the matrix multiplication:

$$y=x^{t} Q x$$

where $$ x = \begin{pmatrix} x_{0} \\ x_{1} \\ x_{2} \end{pmatrix} $$

gives the expression:

$$ y = -x_{0}^{2} -x_{1}^{2} -x_{2}^{2} + 2x_{0}x_{1} + 2x_{1}x_{2} $$

which I want to minimise (EDIT: Changed from maximise - whoops !). As the nodes are either in the maximum set or not, the variables are binary and so $ x_{0}= x_{0}^{2}$ (NB am I correct in thinking for solving the Maximum weighted independent set problem, say where $ x_{0} $ was twice as important, the weight would manifest as a different constant for the quadractic term?)

Is a function that can take a general matrix into DOCPLEX form qiskit can recognise and operate on?

EDIT 2:

I don't think the stable_set.get_operator() method works at all. For the simplest example

path = nx.to_numpy_array(nx.path_graph(3))


qubitOp, offset = stable_set.get_operator(path)
print('Offset:', offset)
print('Ising Hamiltonian:')
print(qubitOp.print_details())

gives:

which is not the correct function - the weights of the different nodes are not equal - and so has no hope of finding the right result.

EDIT 3: After a little bug hunt, the problem is nearly solved ! https://github.com/Qiskit/qiskit-aqua/issues/1553

Answer 1

There is actually a nice way in Qiskit to transform a matrix of an optimization problem into an qubit operator that can be translated into a quadratic program. I'll put here the example, note this is possible for many optimization problems, find every one here in case you want to test something else! For example here it is done for MaxCut and TSP, it could show you with other examples how to do what I did here !

import numpy as np

from qiskit import Aer
from qiskit.optimization.applications.ising import stable_set
from qiskit.aqua.algorithms import VQE, NumPyMinimumEigensolver, QAOA
from qiskit.aqua import aqua_globals
from qiskit.aqua import QuantumInstance
from qiskit.optimization.applications.ising.common import sample_most_likely
from qiskit.optimization.algorithms import MinimumEigenOptimizer
from qiskit.optimization.problems import QuadraticProgram


w = np.array([[-1., 2., 0.],
 [0., -1., 2.],
 [0., 0., -1.]])


qubitOp, offset = stable_set.get_operator(w)
print('Offset:', offset)
print('Ising Hamiltonian:')
print(qubitOp.print_details())

# mapping Ising Hamiltonian to Quadratic Program
qp = QuadraticProgram()
qp.from_ising(qubitOp, offset)
qp.to_docplex().prettyprint()

aqua_globals.random_seed = np.random.default_rng(123)
seed = 10598
backend = Aer.get_backend('statevector_simulator')
quantum_instance = QuantumInstance(backend, seed_simulator=seed, seed_transpiler=seed)

qaoa = QAOA(quantum_instance=quantum_instance, p = 3)

# create minimum eigen optimizer based on qaoa
qaoa_optimizer = MinimumEigenOptimizer(qaoa)

# solve quadratic program
result = qaoa_optimizer.solve(qp)
print(result)

And it will give you this result

Offset: 0.5
Ising Hamiltonian:
IZZ (1+0j)
ZZI (1+0j)
IIZ (0.5+0j)
IZI (0.5+0j)
ZII (-1.5+0j)

// This file has been generated by DOcplex
// model name is: AnonymousModel
// single vars section
dvar bool x_0;
dvar bool x_1;
dvar bool x_2;

minimize
 [ - 3 x_0^2 + 4 x_0*x_1 - 5 x_1^2 + 4 x_1*x_2 + x_2^2 ] + 2;
 
subject to {

}
optimal function value: -3.0
optimal value: [0. 1. 0.]
status: SUCCESS

Hope this helps, please tell me if something is not clear and I'll explain better! :)

Answer 2

Lena's answer doesn't work as we've noted in the comments underneath: the stable_set.get_operator() takes in a numpy adjacency matrix which is not what I'm inputting.

Instead I've written a function which takes the qubo_array I've got and creates a QuadracticProgram exactly.

To do the particular program I was after manually is done as so:

mdl = Model('docplex model')

x_0 = mdl.binary_var('x_0')
x_1 = mdl.binary_var('x_1')
x_2= mdl.binary_var('x_2')

mdl.minimize( -x_0 - x_1  - x_2 + 2*x_0*x_1 + 2*x_2*x_1)


print(mdl.export_as_lp_string())

The way I figured to do this for larger matrices automatically was the use of this function.

def the_auto_doco_mod(qubo_array,model_name,constant):

    """

    Function that takes the   QUBO array created for a graphing problem and converts it to a docplex model
    ready for qiskit

    Directly consrtructs the quadractic program with reference to this page
    """
    number_of_variables = len(qubo_array[1]) # gets the number of variables from the length of the square qubo matrix
    #mdl = Model('model_name')
    mod = QuadraticProgram()

    for variable in range(0,number_of_variables): # creates the binary variables from the size of the matrix 
        var_name = "x_" +str(variable)
        mod.binary_var(name =var_name)

    mod.minimize(constant = 2,quadratic =qubo_array)  # can put in all constraints as quadractic as the binary variables mean that x_0 ^ 2 = x_0 in both cases 
                                                    #  not sure of the impact of this on performance however 

    print(mod.export_as_lp_string())

Does anyone know whether constructing the QuadracticProgram this way will mean that the ising model taken in by qiskit will evaluate slower or be affected in any other way?