QUBO as a QAOAProblem Instance#

In this tutorial you’ll be guided through the process of defining a new phase separator to be used within the scope of the Alternating Operator Ansatz focussed on solving various QUBO problems with only needing the QUBO matrix $Q$ as an input. QUBO, or Quadratic Unconstrained Binary Optimization, is a type of problem that involves optimizing a quadratic function of binary variables.

After first translating the QUBO Hamiltonian $H_{C}$ from the binary basis $x_{i}$ to the state basis with the translation $x_{i} \to \frac{I - Z_{i}}{2}$ , we’ll construct the phase separator unitary and run QAOA to solve the set partitioning problem.

Not to get too ahead of ourselves with how to do it, let’s first show how to solve a QUBO problem using Qrisp. To obtain the optimal solutions to a problem, we only need to know the QUBO matrix $Q$ . In order to get familiar with that, we propose having a look at the Tutorial on Formulating and Using QUBO Models, which explains how to derive such a matrix $Q$ .

One can then simply minimize the cost function $min C (x) = min x^{T} Q x$ for binary variables $x_{i} \in {0, 1}$ . This is done either by sending the QUBO matrix $Q$ to the annealer/quantum computer, which calculates and returns the bitstrings with corresponding values of $C (x)$ . The lower the $C (x)$ , the better the solution. In other words, we want to find the Hamiltonian $H_{C} | x ⟩ = C (x) | x ⟩$ .

Let’s borrow the QUBO matrix (explained in depth in the above mentioned tutorial) for the set partitioning problem. The QUBO matrix in that case is

$Q = (\begin{matrix} - 17 & 10 & 10 & 10 & 0 & 20 \\ 10 & - 18 & 10 & 10 & 10 & 20 \\ 10 & 10 & - 29 & 10 & 20 & 20 \\ 10 & 10 & 10 & - 19 & 10 & 10 \\ 0 & 10 & 20 & 10 & - 17 & 10 \\ 20 & 20 & 20 & 10 & 10 & - 28 \end{matrix})$

Usually, QUBO matrices are upper triangular (by convention) - this means that only elements above the diagonal (with the diagonal included) are not equal to zero. This is because the variables in QUBO problems are binary, which results in the identity $q_{i} q_{j} = q_{j} q_{i}$ for binary variables. Because of this, the elements above the diagonal are doubled compared to the symmetrized version of $Q$ we wrote down above.

Of course it’s easy to go from the conventional QUBO $Q_{up Δ}$ formulation to the symmetrized $Q_{sym}$ QUBO: $Q_{sym} = \frac{1}{2} (Q_{up Δ} + Q_{up Δ}^{T})$ , and back:

import numpy as np
Q_up_triang = 2 * np.triu(Q_sym) - np.diag(np.diag(Q_sym))

Our implementation of solving QUBO problems using QAOA works for both the upper triangular, as well as symmetrized matrix conventions. As promised, we can now immediately solve the QUBO by calling the solve_QUBO method:

from qrisp.qaoa import *
import numpy as np

Q = np.array(
    [
        [-17,  10,  10,  10,   0,  20],
        [ 10, -18,  10,  10,  10,  20],
        [ 10,  10, -29,  10,  20,  20],
        [ 10,  10,  10, -19,  10,  10],
        [ 0,   10,  20,  10, -17,  10],
        [ 20,  20,  20,  10,  10, -28],
    ]
)

solve_QUBO(Q, depth = 1, shots=5000)[:5]

That’s it! You can try running this block of code on the website using the Thebe interactivity integration, or run it on your own Qrispy environment. We see that we obtain the 5 best solutions with the corresponding bitsring of the binary variables. You can test the statements above, and convert the symmetrized matrix into the upper triangular one, and run solve_QUBO again.

We see, that we only needed to provide the matrix $Q$ , specify the depth of the QAOA circuit that is running in the background, and the backend you want to run the algorithm on. Clean and qrispy!

From QUBO to QAOA#

Of course it’s beneficial to not only know the “how”, but also understand the “why”. So let’s dig in!

To construct such a circuit with quantum gates and run it on a quantum computer, one has to translate between the basis of $x_{i}$ to Pauli gate operators.

$Z_{i} | x ⟩ = (- 1)^{x_{i}} \to \frac{1 - Z_{i}}{2} | x ⟩ = x_{i}$

To find $H_{C}$ one can calculate

$\begin{array}{r} \begin{aligned} H_{C} & = \sum_{i, j = 1}^{n} Q_{i j} \frac{1 - Z_{i}}{2} \frac{1 - Z_{j}}{2} \\ = \sum_{i, j = 1}^{n} \frac{Q_{i j}}{4} + \sum_{i, j}^{n} \frac{Q_{i j}}{4} Z_{i} Z_{j} \\ - \sum_{i = 1}^{n} (\sum_{j = 1}^{n} \frac{Q_{i j}}{4}) Z_{i} - \sum_{j = 1}^{n} (\sum_{i = 1}^{n} \frac{Q_{i j}}{4}) Z_{j} \end{aligned} \end{array}$

Swapping the indices $i$ and $j$ in the last sum, and using the identity $Z_{i} Z_{i} = 1$ , we get

$\begin{array}{r} \begin{aligned} H_{C} & = \frac{1}{4} \sum_{i \neq j} Q_{i j} Z_{i} Z_{j} - \frac{1}{4} \sum_{i = 1}^{n} \sum_{j = 1}^{n} (Q_{i j} + Q_{j i}) Z_{i} & + \frac{1}{4} \sum_{i, j = 1}^{n} Q_{i j} + \frac{1}{4} \sum_{i = 1}^{n} Q_{i i} \end{aligned} \end{array}$

Note that for each single $Z_{i}$ we sum the $i$ -th row and the $i$ -th column of the matrix $Q$ .

For the cost operator $U_{C}$ , which we feed into QAOAProblem, we get

$\begin{array}{r} \begin{aligned} U_{C} = e^{- i γ H_{C}} = & \prod_{i, j = 1}^{n} R_{Z_{i} Z_{j}} (\frac{γ}{2} Q_{i j}) \times \prod_{i = 1}^{n} R_{Z_{i}} (- \frac{γ}{2} (\sum_{j = 1}^{n} (Q_{i j} + Q_{j i}))) \\ \times \exp (- \frac{i γ}{4} \sum_{i, j = 1}^{n} Q_{i j} - \frac{i γ}{4} \sum_{i = 1}^{n} Q_{i i}) \end{aligned} \end{array}$

Here, the last factor correspods to a global phase.

Let’s translate this into a function:

def create_QUBO_cost_operator(Q):

    def QUBO_cost_operator(qv, gamma):

        # Rescaling for enhancing the performance of the QAOA
        gamma = gamma/np.sqrt(np.linalg.norm(Q))

        gphase(-gamma/4*(np.sum(Q)+np.trace(Q)),qv[0])
        for i in range(len(Q)):
            rz(-gamma/2*(sum(Q[i])+sum(Q[:,i])), qv[i])
            for j in range(len(Q)):
                if i != j and Q[i][j] != 0:
                    rzz(gamma/2*Q[i][j], qv[i], qv[j])
    return QUBO_cost_operator

Like we did for MaxCut and M $κ$ CS we also define the general QUBO objective function, the classical cost function, as well as construct the QUBOProblem blueprint bringing everything together.

from qrisp import rzz, rz, gphase
import numpy as np

def QUBO_obj(bitstring, Q):
    x = np.array(list(bitstring), dtype=int)
    cost = x.T @ Q @ x
    return cost

def create_QUBO_cl_cost_function(Q):

    def cl_cost_function(counts):

        def QUBO_obj(bitstring, Q):
            x = np.array(list(bitstring), dtype=int)
            cost = x.T @ Q @ x
            return cost

        energy = 0
        for meas, meas_count in counts.items():
            obj_for_meas = QUBO_obj(meas,Q)
            energy += obj_for_meas * meas_count
        return energy

    return cl_cost_function

def QUBO_problem(Q):

    from qrisp.qaoa import QAOAProblem, RX_mixer

    return QAOAProblem(create_QUBO_cost_operator(Q), RX_mixer, create_QUBO_cl_cost_function(Q))

That’s it for the necessary ingredients you learned about in the QAOA theory 101 section! Let’s solve the set partitioning problem from above using this newly acquired information, and combine with how we already ran the QAOA algorithm using the run method:

define the QUBO matrix $Q$ ,
define the quantum argument qarg as a QuantumArray of QuantumVariables,
create the QUBO instance using QUBO_problem we defined above,
run the algorithm using the run method, and last but not least,
examine the QAOA solutions and perform for classical post processing: compute the cost functions, sort the solutions by their cost in ascending order, and print the solutions with their costs.

Warning

For small QUBO instance the number of shots typically exceeds the number of possible solutions. In this case, even QAOA with depth=0, i.e., sampling from a uniform superposition, may yield the optimal solution as the classical post-processing amounts to brute force search! Performance of solve_QUBO for small instance may not be indicative of performance for large instances.

These are exactly the pieces in the mosaic of code that solve_QUBO consists of and performs:

from qrisp.default_backend import def_backend
from qrisp import QuantumVariable, QuantumArray
from operator import itemgetter

Q = np.array(
    [
        [-17,  20,  20,  20,   0,  40],
        [  0, -18,  20,  20,  20,  40],
        [  0,   0, -29,  20,  40,  40],
        [  0,   0,   0, -19,  20,  20],
        [  0,   0,   0,   0, -17,  20],
        [  0,   0,   0,   0,   0, -28],
    ]
)

qarg = QuantumArray(qtype=QuantumVariable(1), shape=len(Q))

QUBO_instance = QUBO_problem(Q)

depth = 1
res = QUBO_instance.run(qarg, depth, mes_kwargs={"backend" : def_backend}, max_iter=50)

costs_and_solutions = [(QUBO_obj(bitstring, Q), bitstring) for bitstring in res.keys()]

sorted_costs_and_solutions = sorted(costs_and_solutions, key=itemgetter(0))

for i in range(5):
    print(f"Solution {i+1}: {sorted_costs_and_solutions[i][1]} with cost: {sorted_costs_and_solutions[i][0]} and probability: {res[sorted_costs_and_solutions[i][1]]}")

Now you are prepared to solve all QUBOs you derive and want to solve. On the other hand, if you would just like to play around instead, try out some QUBOs from this list of QUBO formulations.