Problem with Hamiltonian from Qiskit to Pennylane

PennyLane Help
1

Hi everyone, I would like to make this code work but in PennyLane:

from docplex.mp.model import Model
from qiskit_optimization.problems import QuadraticProgram
from qiskit_optimization.translators import from_docplex_mp
import matplotlib.pyplot as plt 
import networkx as nx

#QAOA
from qiskit.algorithms.minimum_eigensolvers import QAOA
from qiskit.algorithms.optimizers import COBYLA
from qiskit.primitives import Sampler

#VQE
from qiskit.algorithms.optimizers import SPSA
from qiskit.circuit.library import TwoLocal
from qiskit.primitives import Sampler
from qiskit.algorithms.minimum_eigensolvers import SamplingVQE

edges = [(0, 1), (1, 2), (2, 0), (2, 3)]
G = nx.Graph(edges)

nx.draw(G, with_labels=True)

# Parameters of the problem.
a = 4.0
b = 1.0
N = G.number_of_nodes()

# Optimizer
optimizerQAOA = COBYLA()
optimizerVQE = SPSA(maxiter=300)

mdl = Model()

# Here, we define the binary variable X_i.
x = {i: mdl.binary_var(name=f'x_{i}') for i in range(N)}

# For the sake of readability, we have subdivided H_a into three distinct equations,
# with each equation mirroring the one previously explained, partitioned by the addition operator (+).
one = a * mdl.sum((1 - x[u]) * (1 - x[v]) for u in range(N) for v in range(N) if v in G[u])

two = b * mdl.sum(x[v] for v in range(N))

# Here we combine the two equations to form the entire Hamiltonian.
objective = one + two

# We decide to minimize the problem.
mdl.minimize(objective)

# We transform the Hamiltonian equation into a Qubo for Qiskit.
qp = from_docplex_mp(mdl)

# Finally, we create the Ising model with the Qubo for Qiskit.
H, offset = QuadraticProgram.to_ising(qp)
print('Offset:', offset)
print('Ising Hamiltonian:')
print(H)

The two problems are:

  1. If I extract the Pauli-String and give it to PennyLane, the result is not the same and not correct.
  2. The second one is if I use the method implemented by PennyLane for the vertex cover, I don’t obtain the same result. I also saw that if I change the coefficients before cost_h = 3 * edge_driver(graph, ["11", "10", "01"]) + bit_driver(graph_nodes, 0) for another one than 3, the result is completely false. However, with my comprehension of Hamiltonian, the behavior should not be like that.

So, I want to know if it’s me who doesn’t understand something specific or if the way Hamiltonian works is different between Qiskit and PennyLane or something else.

Furthermore, for additional information, I have the same number of operators, but the coefficients differ, and based on my understanding, PennyLane seems to be quite sensitive to coefficients.

Here is my code in PennyLane:

# G = nx.circular_ladder_graph(5)
edges = [(0, 1), (1, 2), (2, 0), (2, 3)]
G = nx.Graph(edges)

nx.draw(G, with_labels=True)

def test1():
    coeffs = [-1 for _ in G.nodes()]
    ops = [qml.PauliZ(w) for w in G.nodes()]
    return qml.Hamiltonian(coeffs, ops)

def test2():
    coeffs = []
    ops = []

    for e in G.edges:
        sign = 1
        coeffs.extend([0.25, 0.25, 0.25])
        ops.extend(
            [
                qml.PauliZ(e[0]) @ qml.PauliZ(e[1]),
                qml.PauliZ(e[0]),
                qml.PauliZ(e[1]),
            ]
        )
    return qml.Hamiltonian(coeffs, ops)

wires = range(len(G.nodes()))
depth = 2

# This is in case you want to try the Pauli string that I got from Qiskit
coeffs = [3.5, 3.5, 5.5, 1.5, 2.0, 2.0, 2.0, 2.0]
ops = [
    qml.PauliZ(0),
    qml.PauliZ(1),
    qml.PauliZ(2),
    qml.PauliZ(3),
    qml.PauliZ(0) @ qml.PauliZ(1),
    qml.PauliZ(0) @ qml.PauliZ(2),
    qml.PauliZ(1) @ qml.PauliZ(2),
    qml.PauliZ(2) @ qml.PauliZ(3),
]

cost_h = (3 * test2()) + test1()
mixer_h = qaoa.x_mixer(G.nodes)


def qaoa_layer(gamma, alpha):
    qaoa.cost_layer(gamma, cost_h)
    qaoa.mixer_layer(alpha, mixer_h)

def circuit(params, **kwargs):
    for w in wires:
        qml.Hadamard(wires=w)
    qml.layer(qaoa_layer, depth, params[0], params[1])

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

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

optimizer = qml.GradientDescentOptimizer()
steps = 70
params = np.array([[0.5, 0.5], [0.5, 0.5]], requires_grad=True)

for i in range(steps):
    params = optimizer.step(cost_function, params)

print("Optimal Parameters")
print(params)

@qml.qnode(dev)
def probability_circuit(gamma, alpha):
    circuit([gamma, alpha])
    return qml.probs(wires=wires)


probs = probability_circuit(params[0], params[1])


plt.style.use("seaborn")
plt.bar(range(2 ** len(wires)), probs)
plt.show()

Thank you in advance for your response.

Name: PennyLane Version: 0.33.0 Summary: PennyLane is a Python quantum machine learning library by Xanadu Inc. Home-page: GitHub - PennyLaneAI/pennylane: PennyLane is a cross-platform Python library for differentiable programming of quantum computers. Train a quantum computer the same way as a neural network. Author: Author-email: License: Apache License 2.0 Location: /opt/homebrew/lib/python3.11/site-packages Requires: appdirs, autograd, autoray, cachetools, networkx, numpy, pennylane-lightning, requests, rustworkx, scipy, semantic-version, toml, typing-extensions Required-by: PennyLane-Lightning Platform info: macOS-13.0-arm64-arm-64bit Python version: 3.11.2 Numpy version: 1.23.5 Scipy version: 1.10.1 Installed devices: - lightning.qubit (PennyLane-Lightning-0.33.1) - default.gaussian (PennyLane-0.33.0) - default.mixed (PennyLane-0.33.0) - default.qubit (PennyLane-0.33.0) - default.qubit.autograd (PennyLane-0.33.0) - default.qubit.jax (PennyLane-0.33.0) - default.qubit.legacy (PennyLane-0.33.0) - default.qubit.tf (PennyLane-0.33.0) - default.qubit.torch (PennyLane-0.33.0) - default.qutrit (PennyLane-0.33.0) - null.qubit (PennyLane-0.33.0)

1 Like
2

Hey @Antoine_Lemelin! Welcome to the forum :slight_smile:

At first glance, it seems like what you want to do should be doable in PennyLane. Looking at your code, I’m a little confused about the details. Can you explain the problem you want to solve in words?

3

Thank you for your response.

I am working on simulating and obtaining the ground state with QAOA of the minimum vertex cover problem. Initially, I implemented it in Qiskit and achieved some good results. However, when I attempted to adapt my code to Pennylane, I noticed that using the Pauli string obtained in Qiskit did not yield the correct results in Pennylane (I got the same operators but not the same coefficient (See exemple above)).

Upon reviewing your Pennylane code for the minimum vertex cover method, I observed that augmenting the constant before the edge_drive method did not produce the correct outcome. In my understanding of Hamiltonian, increasing the constant should only impact the weight of the constraint, not the overall result. This seems to be where my error lies.

4

Hi @Antoine_Lemelin,

Thank you for your question.

When I run print(qml.Hamiltonian(coeffs, ops)) with the coefficients and ops you obtained from Qiskit I get the correct answer, same as in Qiskit. What answer are you getting?

When you run the full example as shown here it does give you the correct answer to the minimum-vertex-cover problem too.

If you run your Hamiltonian obtain with Qiskit in a Qiskit minimum-vertex-cover algorithm, do you get the right states as a result? Can you please share the code for this?

I don’t understand your test 1 and 2, and your definition of cost_h. I would recommend instead using cost_h, mixer_h = qaoa.min_vertex_cover(G, constrained=False)

Please let me know if this answers your question!

5

Thank you for your response.

When I run print(qml.Hamiltonian(coeffs, ops)) with the coefficients and ops obtained from Qiskit, it works. However, when I use it with QAOA, it doesn’t give me the correct answer.

Additionally, if I use the example given above, it works. When I run my Hamiltonian in Qiskit for the minimum vertex cover problem, Qiskit provides the correct answer.

I am bringing up this issue because if I want to implement another problem, such as the Traveling Salesman Problem (TSP), using the paper https://www.frontiersin.org/articles/10.3389/fphy.2014.00005/full, I can’t use just any coefficients; I need to use specific ones. Based on my test yesterday, I observed that when my coefficient exceeds approximately 2, the QAOA method no longer performs well. I suspect this issue may arise because the variable of time in the aproxTimeEvolution function does not consider the ratio of the coefficients. I believe that if we normalized the Hamiltonian before passing to the QAOA method this should work. (This is an hypothesis)

Furthermore, test 1 and test 2 are functions that I use in my code, they are the same as qaoa.min_vertex_cover(G, constrained=False) which I copied to understand the behavior. (Test1 is bit_driver(graph_nodes, 0) and test2 is edge_driver(graph, ["11", "10", "01"]) ) I wanted to figure out how to build a proper Hamiltonian for Pennylane for combinatorial problems.

Here’s my code to run minimum vertex cover problem in Qiskit:
(If you want more exemples for more combinatory problems I can provide you some code without problem)

from docplex.mp.model import Model
from qiskit_optimization.problems import QuadraticProgram
from qiskit_optimization.translators import from_docplex_mp
import matplotlib.pyplot as plt 
import networkx as nx
from qiskit.visualization import plot_histogram

#QAOA
from qiskit.algorithms.minimum_eigensolvers import QAOA
from qiskit.algorithms.optimizers import COBYLA
from qiskit.primitives import Sampler

#VQE
from qiskit.algorithms.optimizers import SPSA
from qiskit.circuit.library import TwoLocal
from qiskit.primitives import Sampler
from qiskit.algorithms.minimum_eigensolvers import SamplingVQE

G = nx.circular_ladder_graph(2)
edges = [(0, 1), (1, 2), (2, 0), (2, 3)]
G = nx.Graph(edges)

nx.draw(G, with_labels=True)

# Parameters of the problem.
a = 3.0
b = 1.0
N = G.number_of_nodes()

#Optimizer
optimizerQAOA = COBYLA()
optimizerVQE = SPSA(maxiter=300)

mdl = Model()

# Here we define the binary variable so the X_i
x = {i: mdl.binary_var(name=f'x_{i}') for i in range(N)}

# For the sake of readability, we have subdivided H_a into three distinct equations, 
# with each equation mirroring the one previously explained, partitioned by the addition operator (+)
one = a* mdl.sum((1 - x[u])*(1 - x[v]) for u in range(N) for v in range(N) if v in G[u])

two = b * mdl.sum(x[v] for v in range(N))

# Here we combined the two equation to form the entire Hamiltonian.
objective = one + two

# We decided to minimize the problem.
mdl.minimize(objective)

# We transform the Hamiltonian euqation into a Qubo for Qiskit.
qp = from_docplex_mp(mdl)

# Finaly we create the ising model with the Qubo for Qiskit.
H, offset = QuadraticProgram.to_ising(qp)
print('Offset:', offset)
print('Ising Hamiltonian:')
print(H)


def printColoredNodes(result):
    colored_nodes = []
    bestresult = result.best_measurement.get('bitstring')[::-1]
    for i in range(N):
        if bestresult[i] == "1":
            colored_nodes.append(i)
    print(colored_nodes)

sampler = Sampler()
qaoa = QAOA(sampler, optimizerQAOA, reps=2)

result = qaoa.compute_minimum_eigenvalue(H)

x = result

print(x)
printColoredNodes(x)
6

Hi @Antoine_Lemelin,

Thanks for adding these additional details!

I’ve been trying to troubleshoot the issue with no success yet. If it’s any help I’ve created some code to build a PennyLane Hamiltonian from the Qiskit one. Below is the code. Remember that you need to map the wires in the reverse order since PennyLane uses big-endian and Qiskit uses little-endian. I’ll get someone else to take a look at this next week to try to figure out what’s going on. Please let me know in case you figure something out on your side.

# Extracting the Qiskit coeffs and Paulis from the ising Hamiltonian
qiskit_coeffs = H.coeffs
ls = H.to_list()
qiskit_paulis = list(zip(*ls))[0]
print(ls)
print(qiskit_paulis)

# Building the Hamiltonian in PL using the Paulis and coeffs obtained from Qiskit
wire_map = {3 : 0, 2 : 1, 1 : 2, 0 : 3} # change the order of the wires since Qiskit uses little endian and PL uses big endian
qiskit_pauli_ops = [qml.pauli.string_to_pauli_word(word,wire_map=wire_map) for word in qiskit_paulis]
H_from_qiskit = qml.Hamiltonian(qiskit_coeffs, qiskit_pauli_ops)
print(H_from_qiskit)
7

Hi @Antoine_Lemelin ,

I have some updates. On one side, I found that even though PennyLane and Qiskit produce different Hamiltonians, the minimum eigenvalues for both are in the same location. This means that they should both converge to the same solution. It also confirms that there’s not a unique Hamiltonian that can represent the problem.

I’ve updated the code I sent you on Friday so that it’s a function:

def qiskit_H_to_PL(H):
  """
  Args: 
      H (SparsePauliOp): Qiskit Ising Hamiltonian
  
  Returns:
      (.Hamiltonian): The Hamiltonian in PennyLane
  """
  # Extract the Qiskit coeffs and Paulis from the Ising Hamiltonian
  qiskit_coeffs = H.coeffs.astype(np.float64).tolist()
  ls = H.to_list()
  qiskit_paulis = list(zip(*ls))[0]

  # Build the Hamiltonian in PL using the Paulis and coeffs obtained from Qiskit
  wire_map = {}
  n = H.num_qubits
  for i in range(n):
    wire_map[n-i-1] = i # change the order of the wires since Qiskit uses little endian and PL uses big endian
  qiskit_pauli_ops = [qml.pauli.string_to_pauli_word(word,wire_map=wire_map) for word in qiskit_paulis]
  H_from_qiskit = qml.Hamiltonian(qiskit_coeffs, qiskit_pauli_ops)
  return H_from_qiskit


H_from_qiskit = qiskit_H_to_PL(H)

I used np.diag of the matrix representing the Hamiltonian to find the classical solution to the problem.

# The smallest eigenvalues show the solution to the problem.
np.diag(qml.matrix(H_from_qiskit))

In this case the minima are located in 6 and 10 (starting from zero), representing the two correct solutions.

The issue is that when you use the variational circuit in the PennyLane intro to QAOA demo the solution is wrong, meaning that it’s not converging to the right value. QAOA can be very tricky and I’m not sure why it’s failing so badly in this case.

I hope this can help you in your project and please let me know if you have any further comments or questions.

1 Like