`cudaq` never giving correct result for `maxcut` QAOA problem.

[sirgeorgesawcon]

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.

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())

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()

I'm never getting the correct answer, I'm always getting 010101. Now, all the parameters used in both Qiskit and Cuda are the same, the optimizer is the same, the circuit is effectively the same and so is the Hamiltonian. But why one is giving the correct result, while cudaq is never giving the correct result.

The same happens when I scale the problem to more number of qubits, Why is this happening? Are there some other parameters that I need to optimize or hyperparameters to set before?

[sirgeorgesawcon]

For other graphs as well, like for a $8$ node graph like this:

whose Hamiltonian is

Offset: -6.0
Ising Hamiltonian:
SparsePauliOp(['IIIIIIZZ', 'IIIIZIIZ', 'IIIZIIIZ', 'IIIZIIZI', 'ZIIIIIZI', 'IIIIZZII', 'IIZIIZII', 'IZIIIZII', 'ZIIIZIII', 'IZIZIIII', 'IZZIIIII', 'ZIZIIIII'],
              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, 0.5+0.j, 0.5+0.j, 0.5+0.j])

I run the following cudaq 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(6)*spin.z(7) + 0.5*spin.z(4)*spin.z(7) + 0.5*spin.z(3)*spin.z(7) + 0.5*spin.z(3)*spin.z(6) + 0.5*spin.z(0)*spin.z(6) + 0.5*spin.z(4)*spin.z(5) + 0.5*spin.z(2)*spin.z(5) + 0.5*spin.z(1)*spin.z(5) + 0.5*spin.z(0)*spin.z(4) + 0.5*spin.z(1)*spin.z(3) + 0.5*spin.z(1)*spin.z(2) + 0.5*spin.z(0)*spin.z(2) 

# Problem parameters.
qubit_count: int = 8
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(13)
optimizer = cudaq.optimizers.COBYLA()
#optimizer.max_iterations = 1000
np.random.seed(13)
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)
# print(dict(sorted(counts.items(),key=lambda item: item[1]))) # print the counts in ascending order
# print the count in descending order
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()

The desired answer is [0, 1, 0, 1, 1, 1, 0, 0] or [0, 1, 0, 1, 0, 1, 1, 0], but using cudaq I'm always getting '01010101': 115, '10101010': 98

Why is it always giving me alternate bits like 010101010 as results, no matter what the graph is?

[leofang]

Hi @sirgeorgesawcon I think you filed in the wrong repo, cudaq issues should go to NVIDIA/cuda-quantum, not NVIDIA/cuQuantum. I don't think any of us have write access there, otherwise we could just transfer this issue over, so could you file the issue there please?

[sirgeorgesawcon]

Thanks @leofang for letting me know, I'll close this issue here, and will open it in NVIDIA/cuda-quantum.