Variational Quantum Eigensolver (VQE) for Number Partitioning Problem

# Import necessary libraries
try:
    import cirq
except ImportError:
    print("installing cirq...")
    !pip install --quiet cirq
    print("installed cirq.")
import cirq
import numpy as np
from scipy.optimize import minimize

# Generate random numbers to partition
numbers = np.random.randint(-10, 100, 6)
n_qubits = len(numbers)
print("Numbers to partition:", numbers)

# Create qubits
qubits = [cirq.NamedQubit(f'q{i}') for i in range(n_qubits)]

# Define a parameterized circuit (ansatz) for VQE
def create_ansatz(params, qubits):
    """
    Create a parameterized quantum circuit for VQE.
    This is a hardware-efficient ansatz with layers of single qubit rotations
    and two-qubit entangling operations.
    """
    circuit = cirq.Circuit()

    # Initial state - Start with all qubits in superposition
    circuit.append(cirq.H.on_each(qubits))

    # Number of layers in the ansatz
    n_layers = len(params) // (n_qubits + n_qubits*(n_qubits-1)//2)

    param_idx = 0
    for _ in range(n_layers):
        # Single qubit rotations (X-rotations)
        for i in range(n_qubits):
            circuit.append(cirq.rx(params[param_idx])(qubits[i]))
            param_idx += 1

        # Two-qubit entangling operations (ZZ interactions)
        for i in range(n_qubits):
            for j in range(i+1, n_qubits):
                if param_idx < len(params):  # Check to avoid index out of range
                    circuit.append(cirq.ZZPowGate(exponent=params[param_idx]/np.pi)(qubits[i], qubits[j]))
                    param_idx += 1

    return circuit

# Set up the simulator
simulator = cirq.Simulator()

# Function to calculate the expectation value of our Hamiltonian
def expectation_value(params):
    """
    Calculate the expectation value of the number partitioning Hamiltonian
    for the state prepared by the parameterized circuit.
    """
    # Create the ansatz circuit
    circuit = create_ansatz(params, qubits)

    # Simulate the circuit to get the state vector
    result = simulator.simulate(circuit)
    state_vector = result.final_state_vector

    # Calculate the expectation value of the Hamiltonian
    # For number partitioning, the Hamiltonian is (sum(a_i * Z_i))^2
    # which expands to sum(a_i^2) + sum_{i<j} 2*a_i*a_j*Z_i*Z_j

    energy = sum(a*a for a in numbers)  # Constant term

    # Calculate the contribution from each ZZ term
    for i in range(n_qubits):
        for j in range(i+1, n_qubits):
            weight = 2 * numbers[i] * numbers[j]

            # Create the Z_i * Z_j operator
            zz_op = cirq.Z(qubits[i]) * cirq.Z(qubits[j])

            # Calculate its expectation value
            zz_val = zz_op.expectation_from_state_vector(
                        state_vector,
                        qubit_map={q: i for i, q in enumerate(qubits)}
                    ).real

            energy += weight * zz_val

    return energy

# Set the number of layers and initialize parameters
n_layers = 2
total_params = n_layers * (n_qubits + n_qubits*(n_qubits-1)//2)
initial_params = np.random.uniform(0, 2*np.pi, total_params)

# Run the VQE optimization
print("Starting VQE optimization...")
result = minimize(
    expectation_value,
    initial_params,
    method='COBYLA',  # Classical optimizer for VQE
    options={'maxiter': 200}
)
optimal_params = result.x
print("Optimization complete!")
print("Final energy:", result.fun)

# Run the optimized circuit and measure
optimized_circuit = create_ansatz(optimal_params, qubits)
optimized_circuit.append(cirq.measure(*qubits, key='result'))
final_result = simulator.run(optimized_circuit, repetitions=1000)
measurements = final_result.measurements['result']

# Analyze the results
def bitstring_to_partition(bits):
    """Convert a bitstring to a partition (1 for subset A, -1 for subset B)."""
    return [1 if bit == 0 else -1 for bit in bits]

def compute_cost(spin_assignment):
    """Compute the cost (squared sum difference) of a partition."""
    S = sum(a * spin for a, spin in zip(numbers, spin_assignment))
    return S * S

# Count the occurrences of each partition and calculate their costs
partition_counts = {}
for bits in measurements:
    partition = tuple(bitstring_to_partition(bits))
    cost = compute_cost(partition)
    if partition in partition_counts:
        partition_counts[partition]['count'] += 1
    else:
        partition_counts[partition] = {'cost': cost, 'count': 1}

# Find the minimum cost partitioning
min_cost = min(info['cost'] for info in partition_counts.values())
optimal_partitions = [p for p, info in partition_counts.items() if info['cost'] == min_cost]

# Display the results
try:
    from tabulate import tabulate
except ImportError:
    !pip install tabulate
    from tabulate import tabulate

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_counts[partition]['count']])

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