Tag: QAOA

QHack 2023: 4 week prep-challenge, day 18, QAOA with custom Hamiltonian

The final exercise on the QML track, included some bug-hunting for me, but this made it much more rewarding, to eventually solve it.

The basic task was conceivably easy: given some data, create a hamiltonian and use a QAOA optimization scheme to obtain a ground state approximation. The problem described some special graph covering problem. The mistakes that took some efforts to resolve where

  • Mistake 1: PennyLane is exceptionally greedy, when arrays appear along the optimization scheme. The way I constructed the hamiltonian included numpy arrays for indexing qubits, PennyLane tried to differentiate the problem by them and vary them subsequently – that lead to exceptions, and it was a real pain to find the source of the error.
  • Mistake 2: I tried to construct operators, instead of decomposing the occupation operator into Pauli matrices. This might work somehow, but I had an easier time decomposing the operators eventually.
  • Mistake 3: I omitted the constant terms, the optimal configuration is not affected by the constant term, but my tests relied on the final energy.

The code for the Hamiltonian

def hamiltonian_coeffs_and_obs(graph):
    """Creates an ordered list of coefficients and observables used to construct
    the UDMIS Hamiltonian.

    Args:
        - graph (list((float, float))): A list of x,y coordinates. e.g. graph = [(1.0, 1.1), (4.5, 3.1)]

    Returns:
        - coeffs (list): List of coefficients for elementary parts of the UDMIS Hamiltonian
        - obs (list(qml.ops)): List of qml.ops
    """

    num_vertices = len(graph)
    E, num_edges = edges(graph)
    u = 1.35
    obs = []
    coeffs = []

    # single terms
    for i in range(num_vertices):
        obs.append(qml.PauliZ(wires=i))
        coeffs.append(-.5)
    
    # constant term
    obs.append(qml.Identity(wires=0))
    coeffs.append(-.5*num_vertices)
    
    # interaction terms
    pairs = np.argwhere(E>0)
    for pair in pairs:
        i = int(pair[0])
        j = int(pair[1])
        obs.append(qml.PauliZ(wires=i) @ qml.PauliZ(wires=j))
        coeffs.append(.25*u)
        coeffs[i] += .25*u
        coeffs[j] += .25*u
        coeffs[num_vertices] += .25*u
    
    return coeffs, obs

The tests check out, success!

QHack 2023: 4 week prep challenge, day 14, First steps in QAOA

Today, I’m not doing an exercise from the repository, but I looking into one of the basic goto algorithms in the world of QC, the Quantum Approximate Optimization Algorithm (QAOA). The code snippets are taken from the very instructive tutorial published here.

A common application of this algorithm is to solve so called Quadratic Unbound Binary Optimization (QUBO) problems. The goal of these problems is to find the binary feature vector x, that minimizes an at most quadratic function f:

\[ f_Q(x) = x^TQx = \sum^n_{i=1}\sum^n_{j=i}Q_{ij}x_ix_j \]

Therefore, the solution to the problem is a series of zeros and ones. As a matter of fact, a large class of problems can be formulated as such a QUBO problem.

QAOA algorithm

The QAOA algorithm employs a very simple idea. A quantum system is first initialized in a promising state (ideally close to the expected global minimum). Then a series of time evolutions is performed, the time evolutions apply alternating the Hamiltonian encoding the problem and a mixer Hamiltonian, to leave meta-stable local optima states.

\[ U(\gamma, \alpha) = \exp(-\alpha_nH_M)\exp(-i\gamma_nH_C)\cdots\exp(-\alpha_1H_M)\exp(-i\gamma_1H_C) \]

In the tutorial, the parameters Alpha and Gamma were optimized using a VQE scheme.

Implementation in PennyLane

First, you’ve got to define a Hamiltonian that characterizes your problem, here’s a toy Hamiltonian

import pennylane as qml

H_P = qml.Hamiltonian(
    [1, 1, 0.5],
    [qml.PauliX(0), qml.PauliZ(1), qml.PauliX(0) @ qml.PauliX(1)]
)
\[ H_P = \sigma_x^{(0)} + \sigma_z^{(1)} + \frac{1}{2} \sigma_x^{(0)} \otimes \sigma_x^{(1)} \]

and your mixer Hamiltonian

\[ H_M = \frac{1}{2} \left( \sigma_x^{(0)} + \sigma_x^{(1)} \right) \]

Then define a time evolution step

def qaoa_layer(params):
    qaoa.cost_layer(params[0], H_P)
    qaoa.mixer_layer(params[1], H_C)

Now, you can combine the initialization and the time evolution sequence to a single circuit

wires = range(4)
depth = 2

def circuit(params, **kwargs):
    for w in wires:
        qml.Hadamard(wires=w)
    qml.layer(qaoa_layer, depth, params[0], params[1])

And having a cost function in place, that evaluates the problem function …

dev = qml.device("qulacs.simulator", wires=wires)

@qml.qnode(dev)
def cost_function(params):
    circuit(params)
    return qml.expval(cost_h)

you can run an optimization procedure, to find a good parameterization for the time evolutions.

optimizer = qml.GradientDescentOptimizer()
steps = 70
params = np.array([[0.5, 0.5], [0.5, 0.5]], requires_grad=True)

for i in range(steps):
    params = optimizer.step(cost_function, params)

print("Optimal Parameters")
print(params)

Interesting observation

The procedure outlined in the tutorial yield good results, providing the highest probabilities for sensible configurations. However, after I went on to modify the algorithm, and increased the number of evolutions from 2 to 5, a less favorable configuration yielded the highest probability. I’ll have to do some research, into why this happens.