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:

Nice Graph

as so:

enter image description here

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?


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:')


enter image description here

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


2 Answers 2


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:')

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

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)

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;

 [ - 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! :)

  • $\begingroup$ Hi Lena! Thanks so much - that's exactly what I was after. It has raised a couple more questions though, if you don't mind! First of all in the Ising Hamiltonian, what does the IZZ (1+0j) represent ? Secondly, how come the function to be minimised is different? I know it has the +2 offset but, even without this, it isn't equivalent to the expression. Finally, I didn't know about the QuantumInstance class, I thought you just put it an operator into the QAOA method. Having a look at the documentation, am I right in thinking QuantumInstance just provides simulation settings ? $\endgroup$
    – jolene
    Commented Feb 17, 2021 at 0:07
  • $\begingroup$ @bubba I'm glad this is what you were looking for! So first, about the Ising hamiltonian, in optimization you usually start with a qubo and deal with bits. To make it understandable for the backend, you need to translate this with something it'll understand, here the Pauli gate Z, and you do this translation via the following bijective transformation ($x_i$ is the bit I, $Z_i$ is the gate applied to the Qubit i) : $x_i = \frac{1-Z_i}{2}$. By applying this to the qubo form, you'll get something like the sum of terms IZZ, ZZI, etc. we have here. Just a way to write the problem for the backend. $\endgroup$
    – Lena
    Commented Feb 17, 2021 at 8:26
  • $\begingroup$ Next, the QuantumInstance is indeed a way to simulate your algorithm in the settings you want. It is quite useful if you want to simulate your Aqua algorithm in "real" conditions without having to use real devices, or play with the transpilation and execution via lots of parameters. $\endgroup$
    – Lena
    Commented Feb 17, 2021 at 8:32
  • $\begingroup$ Finally about the form of the function to minimize, I'll admit I don't know, I'll try to look at it and get back to you here if I find anything useful about that! ;) $\endgroup$
    – Lena
    Commented Feb 17, 2021 at 8:32
  • $\begingroup$ Hi @Lena, looking closer I don't think the use of the stable_set.get_operator() function works in this case. From the documentation, the function takes a numpy adjacency matrix which the matrix `w' sure ain't: it's part of the QUBO formulation. I'm gonna have a look now, but I suspect it might be easier to write a function that takes this form of matrix into a docplex formulation and then use the qiskit converters ..... lots to do ! $\endgroup$
    – jolene
    Commented Feb 23, 2021 at 1:17

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)


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 


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?


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.