sirgeorgesawcon commented 5 months ago

Required prerequisites

[X] Consult the security policy. If reporting a security vulnerability, do not report the bug using this form. Use the process described in the policy to report the issue.
[X] Make sure you've read the documentation. Your issue may be addressed there.
[X] Search the issue tracker to verify that this hasn't already been reported. +1 or comment there if it has.
[X] If possible, make a PR with a failing test to give us a starting point to work on!

Describe the bug

I tried solving a simple $6$ node Max Cut problem , using both qiskit and cudaq. While qiskit usually solves it correctly, cudaq never solved to the right answer. It is always defaulting to 010101.

Steps to reproduce the bug

How to Recreate?

Step 1: Make a MaxCut Graph Problem

import networkx as nx

from qiskit_optimization.applications import Maxcut

seed = 1
num_nodes = 6

G = nx.random_regular_graph(d=3, n=num_nodes, seed=seed)
nx.draw(G, with_labels=True, pos=nx.spring_layout(G, seed=seed))

maxcut = Maxcut(G)
problem = maxcut.to_quadratic_program()
print(problem.prettyprint())

output This is the Quadratic Program for it:

Problem name: Max-cut

Maximize
  -2*x_0*x_1 - 2*x_0*x_3 - 2*x_0*x_4 - 2*x_1*x_2 - 2*x_1*x_5 - 2*x_2*x_3
  - 2*x_2*x_4 - 2*x_3*x_5 - 2*x_4*x_5 + 3*x_0 + 3*x_1 + 3*x_2 + 3*x_3 + 3*x_4
  + 3*x_5

Subject to
  No constraints

  Binary variables (6)
    x_0 x_1 x_2 x_3 x_4 x_5

Step 2: Making a Hamiltonian

qubitOp, offset = qp.to_ising()
print("Offset:", offset)
print("Ising Hamiltonian:")
print(str(qubitOp))

The resulting Hamiltonian is

Offset: -4.5
Ising Hamiltonian:
SparsePauliOp(['IIIIZZ', 'IIZIIZ', 'IZIIIZ', 'IIIZZI', 'ZIIIZI', 'IIZZII', 'IZIZII', 'ZIZIII', 'ZZIIII'],
              coeffs=[0.5+0.j, 0.5+0.j, 0.5+0.j, 0.5+0.j, 0.5+0.j, 0.5+0.j, 0.5+0.j, 0.5+0.j,
 0.5+0.j])

Step 3 : Make a QAOA and Solve;

# QAOA ansatz circuit
ansatz = QAOAAnsatz(hamiltonian, reps=3)
def cost_func(params, ansatz, hamiltonian, estimator):
    """Return estimate of energy from estimator

    Parameters:
        params (ndarray): Array of ansatz parameters
        ansatz (QuantumCircuit): Parameterized ansatz circuit
        hamiltonian (SparsePauliOp): Operator representation of Hamiltonian
        estimator (Estimator): Estimator primitive instance

    Returns:
        float: Energy estimate
    """
    cost = estimator.run(ansatz, hamiltonian, parameter_values=params).result().values[0]
    return cost

x1 = np.random.uniform(- np.pi / 8.0, np.pi/ 8.0 ,ansatz.num_parameters)
res = minimize(cost_func, x1, args=(ansatz, hamiltonian, estimator), method="COBYLA")

The result it gives is: 011010 which is exactly similar to the one I'm getting via Brute Force.

But the moment I make use of cudaq and use the code:

import cudaq
from cudaq import spin
import matplotlib.pyplot as plt
import numpy as np
import time
# Here we build up a kernel for QAOA with `p` layers, with each layer
# containing the alternating set of unitaries corresponding to the problem
# and the mixer Hamiltonians. The algorithm leverages the VQE algorithm
# to compute the Max-Cut of a rectangular graph illustrated below.

#       v0  0---------------------0 v1
#           |                     |
#           |                     |
#           |                     |
#           |                     |
#       v3  0---------------------0 v2
# The Max-Cut for this problem is 0101 or 1010.

# The problem Hamiltonian
#hamiltonian = 0.5 * spin.z(0) * spin.z(1) + 0.5 * spin.z(1) * spin.z(2) + 0.5 * spin.z(0) * spin.z(5) + 0.5 * spin.z(2) * spin.z(3) + 0.5 * spin.z(3) * spin.z(4) + 0.5 * spin.z(4) * spin.z(5) - 3.0
# 10 Node Max Cut
# Set the target to our density matrix simulator.
#cudaq.set_target('density-matrix-cpu')

hamiltonian =  0.5 * spin.z(0) * spin.z(1) + 0.5 * spin.z(0) * spin.z(2) + 0.5 * spin.z(1) * spin.z(3) + 0.5 * spin.z(2) * spin.z(3) + 0.5 * spin.z(0) * spin.z(4) + 0.5 * spin.z(3) * spin.z(4) + 0.5 * spin.z(1) * spin.z(5)  + 0.5 * spin.z(2) * spin.z(5) + 0.5 * spin.z(4) * spin.z(5) - 4.5 
# Problem parameters.
qubit_count: int = 6
layer_count: int = 3
parameter_count: int = 2 * layer_count

def kernel_qaoa() -> cudaq.Kernel:
    """QAOA ansatz for Max-Cut"""
    kernel, thetas = cudaq.make_kernel(list)
    qvec = kernel.qalloc(qubit_count)

    # Create superposition
    kernel.h(qvec)

    # Loop over the layers
    for i in range(layer_count):
        # Loop over the qubits
        # Problem unitary
        for j in range(qubit_count):
            kernel.cx(qvec[j], qvec[(j + 1) % qubit_count])
            kernel.rz(2.0 * thetas[i], qvec[(j + 1) % qubit_count])
            kernel.cx(qvec[j], qvec[(j + 1) % qubit_count])

        # Mixer unitary
        for j in range(qubit_count):
            kernel.rx(2.0 * thetas[i + layer_count], qvec[j])

    return kernel

# Specify the optimizer and its initial parameters. Make it repeatable.
cudaq.set_random_seed(2)
optimizer = cudaq.optimizers.COBYLA()
optimizer.max_iterations = 1000
np.random.seed(20)
optimizer.initial_parameters = np.random.uniform(- np.pi / 8.0, np.pi / 8.0, parameter_count)
#optimizer.initial_parameters = 2.0 * np.pi * np.random.rand(parameter_count)
print("Initial parameters = ", optimizer.initial_parameters)

# Pass the kernel, spin operator, and optimizer to `cudaq.vqe`.
tic = time.time()
optimal_expectation, optimal_parameters = cudaq.vqe(
    kernel=kernel_qaoa(),
    spin_operator=hamiltonian,
    optimizer=optimizer,
    parameter_count=parameter_count)

# Print the optimized value and its parameters
print("Optimal value = ", optimal_expectation)
print("Optimal parameters = ", optimal_parameters)

# Sample the circuit using the optimized parameters
counts = cudaq.sample(kernel_qaoa(), optimal_parameters)
toc = time.time()
print("Time taken = ", toc-tic)
ny_dict = dict(sorted(counts.items(),key=lambda item: item[1], reverse=True))
print(dict(sorted(counts.items(),key=lambda item: item[1], reverse=True)))
#counts.dump()

# plot the histogram of my_dict, for only first 10 elements
# Extract first 10 key-value pairs
first_10_items = list(ny_dict.items())[:10]
x_values = [item[0] for item in first_10_items]
y_values = [item[1] for item in first_10_items]

# Plot the data
plt.figure(figsize=(10, 6))
plt.bar(x_values, y_values, color='skyblue')
plt.xlabel('Keys')
plt.ylabel('Values')
plt.title('Plot of First 10 Key-Value Pairs')
plt.xticks(rotation=90)  # Rotate x-axis labels for better readability
plt.show()

Expected behavior

It should give the correct result. At least some of the times.

Is this a regression? If it is, put the last known working version (or commit) here.

Not a regression

Environment

CUDA Quantum version: 0.6.0
Python version: 3.10
C++ compiler:
Operating system: Ubuntu 22.04.3 LTS

Suggestions

No response

bettinaheim commented 3 months ago

This may be an issue with the COBYLA optimizer that is used. It would be good to confirm if using the scipy cobyla is doing better.

mmvandieren commented 3 months ago

The kernel_qaoa did not match the problem graph and the hamiltonian also did not match the problem graph (e.g. vertex 0 and 2 are not connected so the term 0.5 spin.z(0) spin.z(2) s should not be part of the hamiltonian.

The code below fixes these issues and generates the correct answer.

mmvandieren commented 3 months ago

# Define a function to generate the Hamiltonian for a max cut problem using the graph G

def hamiltonian_max_cut(G): 
    """Hamiltonian for finding the max cut for the graph G

    Parameters
    ----------
    G: networkX graph 
        Graph whose max cut we want to find

    Returns
    -------
    cudaq.SpinOperator
        Hamiltonian for finding the max cut of the graph G
    """

    hamiltonian = 0
    # Since our vertices may not be a list from 0 to n, or may not even be integers,
    # we need to associate  each vertex with its location in the list of nodes for mapping vertices to qubits.
    # To do this we sort the vertices, and then map the first vertex in this sorted list to the qubit index by 0
    # and the next vertex in the sorted list is associated with qubit 1, etc. 
    nodes = sorted(list(nx.nodes(G)))

    for u, v in nx.edges(G):
        # We can use the index() command to read out the qubits associated with the vertex u and v.
        qubitu = nodes.index(u)
        qubitv = nodes.index(v)

        # Add a term to the Hamiltonian for the edge (u,v)
        hamiltonian += 0.5*(spin.z(qubitu)*spin.z(qubitv)-spin.i(qubitu)*spin.i(qubitv))

    return hamiltonian

# Define a kernel function that takes as input two qubits and a parameter

qaoaProblem, qubit_0, qubit_1, alpha = cudaq.make_kernel(cudaq.qubit, cudaq.qubit, float)
qaoaProblem.cx(qubit_0, qubit_1)
qaoaProblem.rz(2*alpha,qubit_1)
qaoaProblem.cx(qubit_0, qubit_1)

# Define a function taking as input the problem graph and layer count
# and returning the QAOA kernel for the graph and given layer count
def kernel_qaoa(G,layer_count) -> cudaq.Kernel:
    """Build the QAOA circuit for max cut of the graph G
    Parameters
    ----------
    G: networkX graph 
        Graph whose max cut we want to find
    layer_count : int 
        Number of layers in the QAOA kernel

    Returns
    -------
    cudaq.Kernel
        QAOA circuit for Max-Cut for max cut of the graph G
    """
    kernel, thetas = cudaq.make_kernel(list)
    nodes = sorted(list(nx.nodes(G)))
    qubit_count = len(nodes)
    qreg = kernel.qalloc(qubit_count)

    # Create superposition
    kernel.h(qreg)

    # Loop over the layers
    for i in range(layer_count):
        # Problem unitary
        for u, v in nx.edges(G):
            qubitu = qreg[nodes.index(u)]
            qubitv = qreg[nodes.index(v)]
            kernel.apply_call(qaoaProblem, qubitu, qubitv, thetas[i])

    # Define a mixer function to apply to each qubit in each layer
    # Keeping this function inside the layer loop so that we can change the parameters
        def mixer(index:int):
            kernel.rx(2.0 * thetas[i + layer_count], qreg[index])

        kernel.for_loop(start =0, stop = qubit_count, function =mixer)

    return kernel

the code below calls on the corrected kernel and hamiltonian functions that match the problem graph

from cudaq import spin
import matplotlib.pyplot as plt
import numpy as np
import time
# Here we build up a kernel for QAOA with `p` layers, with each layer
# containing the alternating set of unitaries corresponding to the problem
# and the mixer Hamiltonians. The algorithm leverages the VQE algorithm
# to compute the Max-Cut of a rectangular graph illustrated below.

#       v0  0---------------------0 v1
#           |                     |
#           |                     |
#           |                     |
#           |                     |
#       v3  0---------------------0 v2
# The Max-Cut for this problem is 0101 or 1010.

# The problem Hamiltonian
#hamiltonian = 0.5 * spin.z(0) * spin.z(1) + 0.5 * spin.z(1) * spin.z(2) + 0.5 * spin.z(0) * spin.z(5) + 0.5 * spin.z(2) * spin.z(3) + 0.5 * spin.z(3) * spin.z(4) + 0.5 * spin.z(4) * spin.z(5) - 3.0
# 10 Node Max Cut
# Set the target to our density matrix simulator.
#cudaq.set_target('density-matrix-cpu')

hamiltonian = hamiltonian_max_cut(G)
# the hamiltonian below does not match the graph G
#hamiltonian =  0.5 * spin.z(0) * spin.z(1) + 0.5 * spin.z(0) * spin.z(2) + 0.5 * spin.z(1) * spin.z(3) + 0.5 * spin.z(2) * spin.z(3) + 0.5 * spin.z(0) * spin.z(4) + 0.5 * spin.z(3) * spin.z(4) + 0.5 * spin.z(1) * spin.z(5)  + 0.5 * spin.z(2) * spin.z(5) + 0.5 * spin.z(4) * spin.z(5) - 4.5 

# Problem parameters.
qubit_count: int = 6
layer_count: int = 3
parameter_count: int = 2 * layer_count

# this kernel_qaoa doesn't match the graph 
#def kernel_qaoa() -> cudaq.Kernel:
#    """QAOA ansatz for Max-Cut"""
#    kernel, thetas = cudaq.make_kernel(list)
#    qvec = kernel.qalloc(qubit_count)
#
#    # Create superposition
#    kernel.h(qvec)#
#
#    # Loop over the layers
#    for i in range(layer_count):
#        # Loop over the qubits
#        # Problem unitary
#        for j in range(qubit_count):
#            kernel.cx(qvec[j], qvec[(j + 1) % qubit_count])
#            kernel.rz(2.0 * thetas[i], qvec[(j + 1) % qubit_count])
#            kernel.cx(qvec[j], qvec[(j + 1) % qubit_count])
#
#        # Mixer unitary
#        for j in range(qubit_count):
#            kernel.rx(2.0 * thetas[i + layer_count], qvec[j])
#
#return kernel

# Specify the optimizer and its initial parameters. Make it repeatable.
cudaq.set_random_seed(2)
optimizer = cudaq.optimizers.COBYLA()
optimizer.max_iterations = 1000
np.random.seed(20)
optimizer.initial_parameters = np.random.uniform(- np.pi / 8.0, np.pi / 8.0, parameter_count)
#optimizer.initial_parameters = 2.0 * np.pi * np.random.rand(parameter_count)
print("Initial parameters = ", optimizer.initial_parameters)

# Pass the kernel, spin operator, and optimizer to `cudaq.vqe`.
tic = time.time()
optimal_expectation, optimal_parameters = cudaq.vqe(
    kernel=kernel_qaoa(G, layer_count),
    spin_operator=hamiltonian,
    optimizer=optimizer,
    parameter_count=parameter_count)

# Print the optimized value and its parameters
print("Optimal value = ", optimal_expectation)
print("Optimal parameters = ", optimal_parameters)

# Sample the circuit using the optimized parameters
counts = cudaq.sample(kernel_qaoa(G, layer_count), optimal_parameters)
toc = time.time()
print("Time taken = ", toc-tic)
ny_dict = dict(sorted(counts.items(),key=lambda item: item[1], reverse=True))
print(dict(sorted(counts.items(),key=lambda item: item[1], reverse=True)))
#counts.dump()

# plot the histogram of my_dict, for only first 10 elements
# Extract first 10 key-value pairs
first_10_items = list(ny_dict.items())[:10]
x_values = [item[0] for item in first_10_items]
y_values = [item[1] for item in first_10_items]

# Plot the data
plt.figure(figsize=(10, 6))
plt.bar(x_values, y_values, color='skyblue')
plt.xlabel('Keys')
plt.ylabel('Values')
plt.title('Plot of First 10 Key-Value Pairs')
plt.xticks(rotation=90)  # Rotate x-axis labels for better readability
plt.show()

It generates the correct answer 010110

mmvandieren commented 3 months ago

there was no issue running this example with the cudaq COBYLA optimizer

zohimchandani commented 3 months ago

@mmvandieren thanks a lot for taking a look at this

sirgeorgesawcon commented 3 months ago

Thank you very much, the code works.

NVIDIA / cuda-quantum

`cudaq` never giving correct result for maxcut QAOA problem #1271

Required prerequisites

Describe the bug

Steps to reproduce the bug

How to Recreate?

Step 1: Make a MaxCut Graph Problem

Step 2: Making a Hamiltonian

Step 3 : Make a QAOA and Solve;

Expected behavior

Is this a regression? If it is, put the last known working version (or commit) here.

Environment

Suggestions