Number Partitioning using QAOA - Variatioal Algorithms

In this notebook, I solve the Number Partitioning problem using the Quantum Approximate Optimization Algorithm (QAOA). The goal of the Number Partitioning problem is to divide a set of numbers into two subsets such that the sums of the numbers in each subset are as close as possible.

Problem¶

Given a set of numbers $S = \{a_1, a_2, ..., a_n\}$ , the objective is to find a partition of $S$ into two subsets $A$ and $B$ such that the absolute difference between the sums of the numbers in each subset is minimized. This can be framed as minimizing the cost function:

$C = (\sum_{i \in A} a_i - \sum_{i \in B} a_i)^2$

This can be further simplified to: $C = (\sum_{i=1}^{n} a_i s_i)^2$ , where $s_i = 1$ if $a_i \in A$ and $s_i = -1$ if $a_i \in B$ .

In the above definition, we can see the subtle hint of spin-up and spin-down states in the form of $s_i$ . More precisely, $\pm 1$ are the eigenvalues of the Pauli $Z$ operator, and we can use them to represent the spin states.

Algorithm¶

The QAOA algorithm is used to find an approximate solution to this optimization problem. The algorithm involves the following:

Define the problem, and prepare the initial state of the qubit system.
Construct a Hamiltonian whose ground state corresponds to the optimal partition.
Construct a mixing Hamiltonian to explore the solution space.
Create a quantum circuit consisting of alternating layers of the Cost and Mixing Hamiltonians.
Use a classical optimizer to find the optimal values for the variational parameters in the quantum circuit.
Measure the output of the quantum circuit to obtain a candidate solution.
Evaluate the quality of the solution and repeat steps 6-7 till we achieve the desired accuracy.

1. Setup and Dependencies¶

This code is meant to be run in Google Colab, so make sure to run it in that environment.

Checks for the cirq library. If not installed, it installs it using pip.
Imports necessary libraries: cirq, sympy, and numpy.

try:
    import cirq
except ImportError:
    print("installing cirq...")
    !pip install --quiet cirq
    print("installed cirq.")
    import cirq
import sympy
import numpy as np

2. Define the Problem¶

Generates a list of random integers (numbers) to be partitioned. The length of the list determines the number of qubits required.

Create Qubits¶

Creates a list of cirq.NamedQubit objects, one for each number in the numbers list.

# Here we use Random numbers, but you can replace this with any other set of numbers
np.seed(42)
numbers = np.random.randint(-10, 100, 6)
n_qubits = len(numbers)
qubits = [cirq.NamedQubit(f'q{i}') for i in range(n_qubits)]

3. Build Cost Hamiltonian¶

The cost Hamiltonian $H_C$ is constructed to penalize unequal partitions.
$H_C = (\sum a_i Z_i)^2 = \sum_i a_i^2 I + \sum_{i \ne j} a_i a_j Z_i Z_j$
The code calculates the constant term and the ZZ terms. The constant_term does not affect the evolution but contributes to the total energy. The cost_terms list stores the ZZ terms, each weighted by $2 * a_i * a_j$ .

cost_terms = []


constant_term = sum(a*a for a in numbers)  # This adds to the total energy but doesn't affect evolution

# Add ZZ terms (i ≠ j): Σ_i,j a_i a_j Z_i Z_j
for i in range(n_qubits):
    for j in range(i+1, n_qubits):
        weight = 2 * numbers[i] * numbers[j]  # Factor of 2 because we're only counting i < j
        cost_terms.append(cirq.ZZ(qubits[i], qubits[j]) ** weight)

Here we have to carefully construct the Hamiltonian to ensure that we are accurately representing the problem. The Hamiltonian is constructed as a sum of Pauli-Z operators, where each term corresponds to a pair of qubits. The coefficients of the terms are determined by the values of the numbers in the input list.

While the implemented cost Hamiltonian penalizes imbalance, alternative cost functions could be used, e.g., one that promotes solutions with similar cardinality of subsets A and B. The current implementation only focuses on minimizing the sum difference, regardless of subset size.

4. Build Mixing Hamiltonian¶

The mixing Hamiltonian $H_B$ is a sum of Pauli-X operators applied to each qubit. It helps to explore the solution space by flipping the qubit states.

# Mixing Hamiltonian: H_B = Σ X_i
mixer = [cirq.X(q) for q in qubits]

5. QAOA Circuit Construction¶

The QAOA circuit is constructed with alternating layers of the cost and mixing Hamiltonians.
Initial Superposition: Applies a Hadamard gate to each qubit to create an equal superposition of all possible states.
Variational Parameters: Defines symbolic variables gamma(γ) and beta(β) for the cost and mixing Hamiltonian evolution times, respectively. (Symbolic variables with sympy are the cirq way of representing parameters for a quantum circuit.)
QAOA Layers: Appends the exponentiated cost and mixing Hamiltonians to the circuit, controlled by the variational parameters.
Measurement: Measures all qubits at the end of the circuit.

# QAOA circuit construction
circuit = cirq.Circuit()

# Initial superposition
circuit.append(cirq.H.on_each(qubits))

# Variational parameters
gamma = sympy.Symbol('γ')
beta = sympy.Symbol('β')

# Add QAOA layers
for term in cost_terms:
    circuit.append(term ** gamma)
for term in mixer:
    circuit.append(term ** beta)

# Measurement
circuit.append(cirq.measure(*qubits, key='result'))

print(circuit)

                                           ┌──────────────────────────┐   ┌──────────────────────────┐   ┌───────────────────────────────────────┐   ┌──────────────────────────┐   ┌───────────────────────────┐
q0: ───H───ZZ──────────────ZZ───────────────ZZ─────────────────────────────ZZ─────────────────────────────ZZ──────────────────────────────────────────X^(β)─────────────────────────────────────────────────────────────────────────────────────────────────M('result')───
           │               │                │                              │                              │                                                                                                                                                 │
q1: ───H───ZZ^(1456.0*γ)───┼────────────────┼────────────ZZ────────────────┼────────────ZZ────────────────┼────────────ZZ─────────────────────────────ZZ─────────────────────────────X^(β)──────────────────────────────────────────────────────────────────M─────────────
                           │                │            │                 │            │                 │            │                              │                                                                                                     │
q2: ───H───────────────────ZZ^(1768.0*γ)────┼────────────ZZ^(1904.0*γ)─────┼────────────┼─────────────────┼────────────┼────────────ZZ────────────────┼────────────ZZ────────────────ZZ─────────────────────────────X^(β)───────────────────────────────────M─────────────
                                            │                              │            │                 │            │            │                 │            │                 │                                                                      │
q3: ───H────────────────────────────────────ZZ^(3588.0*γ)──────────────────┼────────────ZZ^(3864.0*γ)─────┼────────────┼────────────ZZ^(4692.0*γ)─────┼────────────┼─────────────────┼────────────ZZ────────────────ZZ──────────────X^(β)───────────────────M─────────────
                                                                           │                              │            │                              │            │                 │            │                 │                                       │
q4: ───H───────────────────────────────────────────────────────────────────ZZ^(4212.0*γ)──────────────────┼────────────ZZ^(4536.0*γ)──────────────────┼────────────ZZ^(5508.0*γ)─────┼────────────ZZ^(11178.0*γ)────┼───────────────ZZ──────────────X^(β)───M─────────────
                                                                                                          │                                           │                              │                              │               │                       │
q5: ───H──────────────────────────────────────────────────────────────────────────────────────────────────ZZ^(1144.0*γ)───────────────────────────────ZZ^(1232.0*γ)──────────────────ZZ^(1496.0*γ)──────────────────ZZ^(3036.0*γ)───ZZ^(3564.0*γ)───X^(β)───M─────────────
                                           └──────────────────────────┘   └──────────────────────────┘   └───────────────────────────────────────┘   └──────────────────────────┘   └───────────────────────────┘

Here I used a shallow QAOA circuit (one layer). In this problem a single layer is sufficient to find a good approximate solution. In general, the number of layers can be increased to improve the solution quality but at the cost of increased circuit complexity and longer optimization times.

The number of layers is a hyperparameter that can be tuned based on the problem size and complexity.

Improvment: Adaptive QAOA strategies could be employed to optimize the number of layers dynamically based on the problem characteristics.

6. Optimization Setup¶

Sets up the optimization process using cirq.Simulator() and scipy.optimize.minimize().
Energy Calculation: The energy_from_measurements function calculates the cost (energy) of a given set of measurements (bitstrings). It computes $(\sum a_i s_i)^2$ for each bitstring and returns the mean across all bitstrings.
Objective Function: The objective function takes the variational parameters (gamma and beta) as input, using them in the circuit, runs the simulation, and returns the energy.
Optimization: The scipy.optimize.minimize function is used to find the optimal values for gamma and beta that minimize the energy. Here I used COBYLA method optimization.

# Example optimization setup
simulator = cirq.Simulator()
params = [gamma, beta]

def energy_from_measurements(measurements):
    # Calculate (Σ a_i z_i)^2 from bitstrings
    return np.mean([(sum(numbers[i]*(-1)**b for i, b in enumerate(bs)))**2
                    for bs in measurements])

def objective(params):
    resolved = cirq.ParamResolver({'γ': params[0], 'β': params[1]})
    result = simulator.run(circuit, resolved, repetitions=500)
    return energy_from_measurements(result.measurements['result'])

from scipy.optimize import minimize

initial = np.random.uniform(0, np.pi/2, 2)
result = minimize(objective, initial, method='COBYLA')
print("Optimal parameters:", result.x)

Optimal parameters: [ 2.19935745 -0.17078835]

7. Final Result and Analysis¶

Runs the circuit with the optimized parameters to obtain the final result.
The bitstring_to_partition function converts a bitstring to a partition, where 0 represents +1 (subset A) and 1 represents -1 (subset B).
The compute_cost function calculates the cost of a given partition.
The code iterates through the measurement results, converts each bitstring to a partition, computes the cost, and counts the frequency of each partition.
Finally, it identifies the optimal partitions (those with the minimum cost) and prints them along with their frequencies, the elements that belong to each partition and the difference between their sums, using the tabulate library for better readability.

best_params = result.x
final_result = simulator.run(
    circuit,
    cirq.ParamResolver({'γ': best_params[0], 'β': best_params[1]}),
    repetitions=1000
)

import collections


measurements = final_result.measurements['result']
n_qubits = len(numbers)

def bitstring_to_partition(bits):
    """
    Convention: 0 -> +1 (Subset A), 1 -> -1 (Subset B)
    """
    return [1 if bit == 0 else -1 for bit in bits]

def compute_cost(spin_assignment):
    """
    Computes the cost as the square of the imbalance.
    Imbalance S is given by sum(a_i * spin_i).
    """
    S = sum(a * spin for a, spin in zip(numbers, spin_assignment))
    return S * S

# Flatten the array to a list of bitstrings (each is a 1D array)
bitstring_list = [measurement for measurement in measurements]

# Dictionary to count frequency of each unique partition outcome and store its cost.
partition_costs = {}
for bits in bitstring_list:
    # Convert the bitstring into a tuple for dictionary keys.
    partition = tuple(bitstring_to_partition(bits))
    cost = compute_cost(partition)
    if partition in partition_costs:
        partition_costs[partition]['count'] += 1
    else:
        partition_costs[partition] = {'cost': cost, 'count': 1}

# Determine the partition(s) with the minimal cost.
min_cost = min(info['cost'] for info in partition_costs.values())
optimal_partitions = [p for p, info in partition_costs.items() if info['cost'] == min_cost]

try:
    from tabulate import tabulate
except ImportError:
    !pip install tabulate
    from tabulate import tabulate

# print("Minimum imbalance (cost):", min_cost)
# print("Optimal partition(s) and corresponding frequencies:")
# print(numbers)
# for partition in optimal_partitions:
#     sub_a = [numbers[i] for i in range(n_qubits) if partition[i] == 1]
#     sub_b = [numbers[i] for i in range(n_qubits) if partition[i] == -1]
#     print("Partition:", partition,
#           "First Partition", sub_a,
#           "Second Partition", sub_b,
#           "Sum Difference", sum(sub_a) - sum(sub_b),
#           "Frequency:", partition_costs[partition]['count'])

table_data = []
for partition in optimal_partitions:
    sub_a = [numbers[i] for i in range(n_qubits) if partition[i] == 1]
    sub_b = [numbers[i] for i in range(n_qubits) if partition[i] == -1]
    table_data.append([partition, sub_a, sub_b, sum(sub_a) - sum(sub_b), partition_costs[partition]['count']])

headers = ["Partition", "First Partition", "Second Partition", "Sum Difference", "Frequency"]
print(tabulate(table_data, headers=headers, tablefmt="grid"))

+-----------------------+-------------------+--------------------+------------------+-------------+
| Partition             | First Partition   | Second Partition   |   Sum Difference |   Frequency |
+=======================+===================+====================+==================+=============+
| (-1, 1, 1, 1, -1, -1) | [28, 34, 69]      | [26, 81, 22]       |                2 |          29 |
+-----------------------+-------------------+--------------------+------------------+-------------+
| (1, -1, 1, 1, -1, -1) | [26, 34, 69]      | [28, 81, 22]       |               -2 |           7 |
+-----------------------+-------------------+--------------------+------------------+-------------+
| (-1, 1, -1, -1, 1, 1) | [28, 81, 22]      | [26, 34, 69]       |                2 |           3 |
+-----------------------+-------------------+--------------------+------------------+-------------+
| (1, -1, -1, -1, 1, 1) | [26, 81, 22]      | [28, 34, 69]       |               -2 |          19 |
+-----------------------+-------------------+--------------------+------------------+-------------+

We have some assumptions¶

The optimal partitions obtained from the QAOA simulation are good approximations of the true optimal partitions. Since QAOA is a variational algorithm, it may not always find the exact optimal solution, especially for larger problem sizes. However, it can provide near-optimal solutions in a reasonable amount of time for the near-term quantum computers.

Improvements¶

Deeper circuits can potentially lead to better solutions.
Techniques such as variational quantum eigensolver (VQE) could be used to improve the accuracy of the results.
Here I considered a list of numbers with only 6 elements. The code can be easily modified to handle larger lists, but the optimization time will increase significantly. For larger lists, we can use more advanced techniques such as QAOA with adaptive depth or hybrid quantum-classical algorithms to improve the performance.

Variatioal Algorithms

Variational Quantum Eigensolver (VQE) for Number Partitioning Problem