# Qiskit: Taking a QUBO matrix into qubit_op'

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

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

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.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)

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

• 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 ? Feb 17 at 0:07
• @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.
– Lena
Feb 17 at 8:26
• 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.
– Lena
Feb 17 at 8:32
• 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! ;)
– Lena
Feb 17 at 8:32
• 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 ! Feb 23 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)

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

"""
number_of_variables = len(qubo_array[1]) # gets the number of variables from the length of the square qubo matrix
#mdl = Model('model_name')