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.