Efficient Hamiltonian Sum for QAOA

PennyLane Help
1

Hello everyone,

I am trying to optimise the code to construct a Hamiltonian as a sum of Hamiltonians. I have notised that using qml.Hamiltonian with an array of coefficients and an array of observables is by far faster (running time in milliseconds even for 15 or more qubits) than constructing it as a sum of coefficient times observable (running time in minutes for 15 or more qubits). The code below constructs Hamiltonians using both methods.

import pennylane as qml

N = 2

s = [(1/2 * (qml.Identity(i) - qml.PauliZ(i))) for i in range(2*N)]
s_plus = s[:N]
s_minus = s[N:]

cost_h_from_coeffs_obs = qml.Hamiltonian(coeffs=[1 for i in range(4)],
                                         observables=[s_plus[0],
                                                      (s_plus[1] - 2*(s_plus[1] @ s_minus[1]) + s_minus[1]),
                                                      (s_plus[0] - s_minus[0]),
                                                      (s_plus[0] - s_minus[0]) @ (s_plus[1] - s_minus[1])])
print("Hamiltonian from coeffs. obs. construction")
print(cost_h_from_coeffs_obs)

cost_h_from_sum = sum([s_plus[0],
                       (s_plus[1] - 2*(s_plus[1] @ s_minus[1]) + s_minus[1]),
                       (s_plus[0] - s_minus[0]),
                       (s_plus[0] - s_minus[0]) @ (s_plus[1] - s_minus[1])])
print("Hamiltonian from sum")
print(cost_h_from_sum)

print("Equivalent Hamiltonians?", cost_h_from_coeffs_obs.compare(cost_h_from_sum))

However, when I run this code, it returns that both Hamiltonians are not equivalent. Moreover, when I print them, I see the PauliZ and Identity gates for the summed up Hamiltonian, but not for the one computed with qml.Hamiltonian.

Hamiltonian from coeffs. obs. construction
  (1) [Hamiltonian0]
+ (1) [Hamiltonian1,3]
+ (1) [Hamiltonian0,2]
+ (1) [Hamiltonian0,1,2,3]
Hamiltonian from sum
  (-1.0) [Z0]
+ (0.5) [Z2]
+ (1.0) [I0]
+ (-0.5) [Z1 Z3]
+ (-0.25) [Z0 Z3]
+ (-0.25) [Z2 Z1]
+ (0.25) [Z0 Z1]
+ (0.25) [Z2 Z3]
Equivalent Hamiltonians? False

Can someone explain me why these two do not match?

Finally, when I try to use the Hamiltonian computed with qml.Hamiltonian for QAOA, I get the following error message:

ValueError: hamiltonian must be written only in terms of PauliZ and Identity gates

How can I deal with it? Is there any other way to speed up the computation of a cost Hamiltonian from the sum of s_plus and s_minus elements times some fixed factors? Thank you very much in advance.

2

Hi @pormelrog, thanks for the question!

The two objects still represent the same mathematical Hamiltonian, but since one is formed by qml.Hamiltonian, and the other by operator arithmetic directly, they will have slightly different representations as objects in Python.

If you examine your cost_h_from_sum closely, it contains all the same information as cost_h_from_coeffs_obs:

cost_h_from_coeffs_obs.terms()
# should give the following:
# (array([1, 1, 1, 1]),
# [<Hamiltonian: terms=2, wires=[0]>,
#  <Hamiltonian: terms=2, wires=[1, 3]>,
#  <Hamiltonian: terms=2, wires=[0, 2]>,
#  <Hamiltonian: terms=4, wires=[0, 1, 2, 3]>])

# examine each of these sub-terms and see what they contain:

cost_h_from_coeffs_obs.terms()[1][0]
# gives
#  (-0.5) [Z0]
# + (0.5) [I0]

cost_h_from_coeffs_obs.terms()[1][1]
# gives: 
#  (0.5) [I3]
# + (-0.5) [Z1 Z3]

cost_h_from_coeffs_obs.terms()[1][2]
# gives: 
#   (-0.5) [Z0]
#  + (0.5) [Z2]


cost_h_from_coeffs_obs.terms()[1][3]
# gives: 
#  (-0.25) [Z0 Z3]
# + (-0.25) [Z2 Z1]
# + (0.25) [Z0 Z1]
# + (0.25) [Z2 Z3]

If you compare these terms to the other format, you’ll see that everything is accounted for (notice that there are two separate Z0 contributions, and the identity operators I0 and I3 can be safely summed together (there are implicit identities on all other wires).

So using these both in a circuit, you should expect the same numerical answers (possibly with somewhat different execution efficiencies :))

Now as for your QAOA question, QAOA cost Hamiltonians need to be diagonal in the Z basis (this is a unique structure of QAOA), so they can only consist of identities or PauliZ ops.

1 Like
3

For those looking for a solution to the problem I described in my original post: you can take advantage of the qml.Hamiltonian() method to construct Hamiltonians as the sum of other Hamiltonians times some factors using the following code.

import pennylane as qml
import numpy as np

# Example of Hamiltonians and their factors to sum up.
N = 4
# Replace the obs_to_sum array with the list of Hamiltonians which you would like to sum-up.
obs_to_sum = [(1/2 * (qml.Identity(i) - qml.PauliZ(i))) for i in range(N)]
# Replace the coeffs_to_sum array with a list of the factors with which each Hamiltonian in obs_to_array will be multiplied in the sum.
coeffs_for_sum = [i for i in range(N)]
# Check that there is one multiplication factor for each Hamiltonian and one Hamiltonian for each multiplication factor.
assert(len(obs_to_sum) == len(coeffs_for_sum))

# The methods below compute the Hamiltian sum([coeffs_to_sum[i]*obs_to_sum[i] for i in range(len(obs_to_sum))]).

def SimplePythonSum(obs_to_sum, coeffs_for_sum):
    # Check that there is one multiplication factor for each Hamiltonian and one Hamiltonian for each multiplication factor.
    assert(len(obs_to_sum) == len(coeffs_for_sum))

    # Sum up all Hamiltonians in obst_to_sum times their respective coefficient from coeffs_for_sum.
    cost_h = sum([coeffs_for_sum[i]*obs_to_sum[i] for i in range(len(obs_to_sum))])
    # Return summed Hamiltonian.
    return cost_h

def SumWithQMLHamiltonian(obs_to_sum, coeffs_for_sum):
    # Check that there is one multiplication factor for each Hamiltonian and one Hamiltonian for each multiplication factor.
    assert(len(obs_to_sum) == len(coeffs_for_sum))

    # Create a new list to save the observables.
    obs = []
    # Create a new numpy array to save the coefficients.
    coeffs = np.array([])

    # Extract the coefficients and observables from each Hamiltonian in obst_to_sum.
    for i in range(len(obs_to_sum)):
        # Save the extracted Hamiltonians from the ith obs_to_sum Hamiltonian.
        obs = obs + obs_to_sum[i].terms()[1]
        # Save the extracted coefficients from the ith obs_to_sum Hamiltonian and multiply them by the corresponding factor for the ithobs_to_sum Hamiltonian.
        coeffs = np.concatenate((coeffs, coeffs_for_sum[i] * obs_to_sum[i].terms()[0]))

    # Check that there is one multiplication factor for each Hamiltonian and one Hamiltonian for each multiplication factor.
    assert(len(obs) == len(coeffs))

    # Create a new Hamiltonian from the extracted observables and the adjusted coefficients.
    cost_h = qml.Hamiltonian(coeffs=coeffs,
                             observables=obs)
    # Return summed Hamiltonian.
    return cost_h

cost_h_from_coeffs_obs = SumWithQMLHamiltonian(obs_to_sum, coeffs_for_sum)
print("Hamiltonian from coeffs. obs. construction")
print(cost_h_from_coeffs_obs)
print("Simplified Hamiltonian from coeffs. obs. construction")
print(cost_h_from_coeffs_obs.simplify())

cost_h_from_sum = SimplePythonSum(obs_to_sum, coeffs_for_sum)
print("Hamiltonian from sum")
print(cost_h_from_sum)

print("Equivalent Hamiltonians?", cost_h_from_coeffs_obs.compare(cost_h_from_sum))

This code snippet returns the following:

Hamiltonian from coeffs. obs. construction
  (-1.5) [Z3]
+ (-1.0) [Z2]
+ (-0.5) [Z1]
+ (-0.0) [Z0]
+ (0.0) [I0]
+ (0.5) [I1]
+ (1.0) [I2]
+ (1.5) [I3]
Simplified Hamiltonian from coeffs. obs. construction
  (-1.5) [Z3]
+ (-1.0) [Z2]
+ (-0.5) [Z1]
+ (-0.0) [Z0]
+ (3.0) [I0]
Hamiltonian from sum
  (-1.5) [Z3]
+ (-1.0) [Z2]
+ (-0.5) [Z1]
+ (-0.0) [Z0]
+ (3.0) [I0]
Equivalent Hamiltonians? True

The main advantage of the SumWithQMLHamiltonian method defined above is that it runs substantially faster than the SimplePythonSum method, which computes the Hamiltonian using Python’s sum method. I have tested this on my Asus UX410UAK laptop with an i7-7500U CPU and 8GB of RAM.

N = 100
obs_to_sum = [(1/2 * (qml.Identity(i) - qml.PauliZ(i))) for i in range(N)]
coeffs_for_sum = [i for i in range(N)]
assert(len(obs_to_sum) == len(coeffs_for_sum))

%timeit SumWithQMLHamiltonian(obs_to_sum, coeffs_for_sum)
%timeit SimplePythonSum(obs_to_sum, coeffs_for_sum)

The result:

1.56 ms ± 82.9 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
3.02 s ± 43.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

As you can see then SumWithQMLHamiltonian method is clearly the preferred method to compute the sum of Hamiltonians. An explanation why the code from my original post did not work for QAOA and why the observables have to be extracted is provided in the answer provided by @nathan.

1 Like
4

Thank you very much for the answer, @nathan. Is there any reason why the qml.Hamiltonian() method does not extract the coefficients and observables of the input Hamiltonians automatically/internally? As you can see in the code I have posted as answer to this question, doing that is pretty useful when computing a Hamiltonian as the sum of other Hamiltonians times some factors, specially in terms of computational time. Do you think this “feature” could be added to the qml.Hamiltonian() method?

5

Hi @pormelrog!

Thanks for pointing this out! This is something on our radar at the moment - we’re working on getting it so that the “simple python sum” style is both the recommended approach in PennyLane and is as fast as the nested-Hamiltonians approach :+1: Unfortunately, this might not be available for one or two releases, so the SumWithQMLHamiltonian function you created should be a good short-term solution for what you need.

2 Likes
6

@Tom_Bromley
It appears that this has not been fixed yet. Would you be willing to accept a PR? If yes, can you point me in the right direction to start working on it?

7

Hi @patelvyom, welcome to the Forum!

I’m not sure whether work on this is already in progress. Let me check and get back to you with an answer in the next couple of days.

1 Like
8

Welcome to the community, @patelvyom.

9

Hi @patelvyom,

I wanted to let you know that this is part of a big internal project that is still a few releases away. Normally we encourage PRs but this doesn’t seem like the right one. If this is blocking work for you then it’s best if you create a Fork or try to circumvent the problem. If you were looking to make an interesting contribution and you have other ideas I can let you know if they’re a better fit at the moment! You can also take on one of the Good First Issues which are already well-defined projects that you can work on.

Let me know if you have any questions and thank you once again for proposing to work on a PR for this.