Using QAOA to find solutions

Using QAOA to find solutions¶

We start by importing the version of Numpy provided by Pennylane.

from pennylane import numpy as np

The following line sets the parameter $N$ of the $N$ queens problem. This parameter only needs to be set here, everything downstream is written in terms of this $N$.

N = 4

The exact cover problem and its generalization¶

Given a (countable) set $X$ and a collection $S$ of subsets of $X$, a subcollection $S^\star$ of $S$ is called an exact cover of $X$ if any two sets in $S^\star$ are disjoint, and the union of all sets in $S^\star$ is $X$. We need to find such an exact cover if it exists. We can represent this in matrix form as follows. Let the matrix $M$ have as many rows as there are sets in $S$, and as many columns as there are elements in $X$. For each row of the matrix (corresponding to each set $s$ in $S$), let the $i$^th element be $1$ if the corresponding element is in $s$ and zero otherwise. Then, the objective is to find a set of rows such that their sum is the all-ones vector. The set $X$ can be thought of as a set of constraints, and the choice of rows as selections. The objective is then to select a set of rows so that each constraint is satisfied by exactly one selection.

The generalized exact cover problem divides $X$ up into two sets. In one set, the constraints must be satisfied by exactly one selection (these are called the primary constraints), while the secondary constraints may be satisfied by at most one (that is, either zero or one) selection. In matrix language, we need to find a set of rows that sum to $1$ in the columns corresponding to the primary constraints and either $0$ or $1$ in the columns corresponding to the secondary constraints.

We can frame the $N$ queens problem in this framework as follows. Let the matrix $M$ have $N^2$ columns and $6N - 6$ columns. The first $N$ columns corresponds to the $N$ files (columns) of the chessboard. The next $N$ corresponds to the $N$ ranks (rows). The next $2N - 3$ columns correspond to the diagonals, with the first and last diagonal omitted because they only consist of one square each. Similarly, the last $2N - 3$ columns correspond to the reverse diagonals. Each row $M$ corresponds to a particular position of a queen on the board; it has ones in the columns corresponding to the rank, file, diagonal and reverse diagonal that the square is in, and zeros everywhere else. The objective is to find a selection of rows (or, equivalently, positions of queens) such that they sum to $1$ for each column corresponding to ranks and files (because each rank and file must contain a queen), and either $0$ for $1$ for each column corresponding to a diagonal or a reverse diagonal (because it is not necessary for each diagonal/reverse diagonal to contain a queen).

This matrix is generated by the following code. For each row in the matrix (which goes from $1$ to $N^2$ corresponding to the $N^2$ choices of squares), it places a $1$ in the index for the rank, file, diagonal and reverse diagonal of the square, and zeros everywhere else.

M = np.zeros((N*N, 6*N - 6), requires_grad=False)
for m in range(np.shape(M)[0]):
    for n in range(np.shape(M)[1]):
        file = m // N
        rank = m % N
        diagonal = rank - file + (N-1) - 1
        rdiagonal = rank + file - 1
        if ((n == file 
            or n == N + rank 
            or n == 2*N + diagonal 
            or n == 4*N - 3 + rdiagonal
            )
            and diagonal >= 0 
            and diagonal < 2*N - 3 
            and rdiagonal >= 0 
            and rdiagonal < 2*N - 3
            ):
            M[m][n] = 1
        if n == file or n == N + rank:
            M[m][n] = 1
        if diagonal >= 0 and diagonal < 2*N - 3:
            if n == 2*N + diagonal:
                M[m][n] = 1
        if rdiagonal >= 0 and rdiagonal < 2*N - 3:
            if n == 4*N - 3 + rdiagonal:
                M[m][n] = 1
M

As shown in [1], a generalized exact cover problem can be reduced to an exact cover problem by adding a row for each secondary constraint with a $1$ in the corresponding column and zeros everywhere else. The solution to the generalized problem is obtained by taking the solution of the exact problem and picking the original rows that are selected. The following code takes the matrix above defining the generalized problem and creates the matrix for the corresponding exact problem.

concat = np.concatenate((np.zeros([4*N - 6, 2*N]), np.eye(4*N - 6)), axis=1)
M = np.concatenate((M, concat), axis=0)
M

Setting up the QAOA¶

To apply the QAOA, we need to turn the above problem into one of finding the ground state of an appropriate Hamiltonian. In [2], it is shown how to find the relevant Hamiltonian starting from the matrix defining an exact cover problem. The Hamiltonian is given by $$ H = \sum_{i < j} J_{ij} \sigma_i^z \sigma_j^z + \sum_{i} h_i \sigma_i^z, $$ where $$ J_{rr'} = \frac{1}{2} \sum_{j} M_{rj}M_{r'j} $$ and $$ h_r = \frac{1}{2} \sum_{i} M_{ri} \left(\sum_{r'} M_{r'i} - 2\right). $$

The following lines compute the matrix $J$ from the matrix $M$ and checks that $J$ is symmetric ($J$ must be Hermitian, and is real since $M$ only has real entries.).

rows = np.shape(M)[0]
cols = np.shape(M)[1]
J = np.zeros((rows, rows), requires_grad=False)
for i in range(rows):
    for j in range(rows):
        J[i][j] = (0.5)*np.sum([M[i][f] * M[j][f] for f in range(cols)])
np.allclose(J, np.transpose(J))

The following lines construct the vector $h$ from the matrix $M$.

h = np.zeros(rows, requires_grad=False)
for r in range(rows):
    h[r] = (0.5)*np.sum([M[r][f]*(np.sum([M[s][f] for s in range(rows)]) - 2) for f in range(cols)])
h

We now have everything in place for using QAOA. We need to create the cost and mixer Hamiltonians. We first begin by defining the cost Hamiltonian using the $J$ and $h$ we defined above.

import pennylane as qml

cost_coeffs = []
cost_observables = []
for j in range(np.shape(J)[0]):
    for i in range(j-1):
        cost_coeffs.append(J[i][j])
        cost_observables.append(qml.PauliZ(i) @ qml.PauliZ(j))
for j in range(np.shape(h)[0]):
    cost_coeffs.append(h[j])
    cost_observables.append(qml.PauliZ(j))
cost_hamiltonian = qml.Hamiltonian(cost_coeffs, cost_observables, simplify=True)
cost_hamiltonian

The mixer coefficients consist of Pauli $X$ gates acting on the qubits.

mixer_coeffs = []
mixer_observables = []
for r in range(rows):
    mixer_coeffs.append(1)
    mixer_observables.append(qml.PauliX(r))
mixer_hamiltonian = qml.Hamiltonian(mixer_coeffs, mixer_observables)
mixer_hamiltonian
    

Optimization¶

We shall use the qaoa module from Pennylane and define a layer of the QAOA circuit.

from pennylane import qaoa

def qaoa_layer(params):
    qaoa.cost_layer(params[0], cost_hamiltonian)
    qaoa.mixer_layer(params[1], mixer_hamiltonian)

Here we set the depth of the QAOA circuit. As with N above, everything downstream is written in terms of this parameter and so to control the number of depths, the only change to be made is here.

DEPTH = 1

For the circuit, we start with a uniform superposition over the starting qubits and then apply the cost and mixer circuits in succession, as usual for the QAOA.

wires = range(rows)
depth = DEPTH

def circuit(params):
    for w in wires:
        qml.Hadamard(wires=w)
    qml.layer(qaoa_layer, depth, params)

The cost function is simply the expectation value of the cost Hamiltonian defined above.

dev = qml.device("default.qubit", wires=wires)

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

The parameters are initialized to $0.5$ each (we have not investigated other starting parameter values). We then run the optimizer for $30$ steps using Pennylane.

optimizer = qml.GradientDescentOptimizer()
steps = 30
params = np.array([[0.5, 0.5] for i in range(depth)], requires_grad=True)
for i in range(steps):
    params = optimizer.step(cost_function, params)
    print(i, cost_function(params))

print("Optimal Parameters")
print(params)

Next, we use the optimal parameters and sample the qubits corresponding to the rows of the original generalized problem. This data is stored in the positions list.

run_dev = qml.device("default.qubit", wires=wires, shots=1)

@qml.qnode(run_dev)
def optimized_circuit(params):
    circuit(params)
    return qml.sample(wires=[i for i in range(N*N)])

positions = optimized_circuit(params)

Finally, we create the $N \times N$ chessboard with the queens in the computed positions.

for i in range(N):
    for j in range(N):
        if positions[N*i + j] == 1:
            print('🟥', end='')
        else:
            if (i+j) % 2 == 0:
                print('⬛', end='')
            else:
                print('⬜', end='')
    print('')
    

References¶

[1]: Knuth, Donald E. “Dancing links.” arXiv preprint cs/0011047 (2000).

[2]: Vikstål, Pontus, Mattias Grönkvist, Marika Svensson, Martin Andersson, Göran Johansson, and Giulia Ferrini. “Applying the quantum approximate optimization algorithm to the tail-assignment problem.” Physical Review Applied 14, no. 3 (2020): 034009.

[3]: Farhi, Edward, Jeffrey Goldstone, and Sam Gutmann. “A quantum approximate optimization algorithm.” arXiv preprint arXiv:1411.4028 (2014).

Introduction

N-queens using QAOA