Expectation value of a Hamiltonian

[AdrianPerezSalinas]

Hi all,

I am working on a new algorithm that is to be implementable in an actual experiment. In order to get the feeling about what is going on, I would need to compute the expectation value for a given Hamiltonian with shots and not with the exact wavefunction of the state. In other words, it would be like measuring in the real physical world. Is it possible?

Thank you all

[scarrazza]

Do you have references about this feature?

[AdrianPerezSalinas]

I do not have any, but I don't think it is necessary.
Let us think in the Z hamiltonian for one qubit. We can perform the state * hamiltonian * state matrix multiplication and obtain a number. But we could also measure with shots the probability of obtaining 0 (p0) and 1 (p1) and say Z = p0 - p1

Since we measure in the Z basis, it would be immediate in this case. For the X and Y basis, I would apply a Hadamard gate or the equivalent in Y to measure in another basis

[stavros11]

Thanks for opening this, it is certainly not possible to do directly with Qibo. For now, I will try to find a temporary solution using custom code and will write here if I come up with something.

It would be great to add it as a feature. I agree that a reference is not needed, I am just not sure how easy it would be to implement this with the current representation we have for Hamiltonians, which is mainly based on numpy arrays. We may need redesign this a bit.

[stavros11]

I found a relatively simple way to do this using our existing TrotterHamiltonian functionality. It may be considered cheating because the state vector is used in the calculation, however it may be possible to use it to emulate the "measurement noise" in the expectation value. I am not 100% sure if it is a good approach though. The idea is the following:

# construct the circuit that produces the measurements in the Z basis
c = models.Circuit(nqubits)
c.add(...) # all your gates
# measure all qubits in the Z basis
c.add(gates.M(*range(nqubits)))
# execute
result = c(nshots=nshots)
# get the decimal samples as np.array
measurements = result.samples(binary=False).numpy()
# get the final state vector as np.array
# note this is the state before the measurement (not collapsed)
final_state = c.final_state.numpy()

# construct your Hamiltonian as `TrotterHamiltonian`, eg.
ham = hamiltonians.TFIM(nqubits, trotter=True)

# calculate the EV using measurements
ev = 0
psiloc0 = final_state[measurements]
# add Hamiltonian terms one-by-one
for term in ham.terms():
    psiloc = term(np.copy(final_state)).numpy()[measurements]
    # np.copy is important here, otherwise the custom operator will alter the state
    # and the remaining terms will be wrong!
    ev += (psiloc / psiloc0).mean()

# (optional) calculate exact EV to compare
exact_ev = ham.expectation(final_state)

If you try this you will see that ev converges to exact_ev as you increase nshots which is what is expected. I used TFIM in this example but you can create your own Hamiltonian as TrotterHamiltonians relatively easily using sympy symbols (see our example).

Perhaps it is not exactly what we should implement in Qibo but let me know if it helps.

[stavros11]

In order to make clearer what the code in the previous post does, the idea is based in the following equation for the expectation value of an arbitrary operator A:

Note that while \sum _\sigma denotes a sum over all possible bit-strings/spin configurations, the last sum is a sum over samples that follow the probability distribution defined by the wave function. This is the typical averaging approach used in Monte Carlo methods and in our case we can use the measurements as the MC samples (instead of doing the standard Metropolis). To be more explicit, in the code psiloc0 is the final state indexed by the measurement samples (that is <sigma | psi >) and psiloc the same after A has been applied (that is <sigma | A |psi >). So ev calculates the sum from the equation above using every local term of the Hamiltonian as your A.

As I said, this approach still uses the full state vector so it is not the one you would use in the quantum computer, however I think it can give an estimate of the statistical noise in the EV, since it is not an exact matrix multiplication calculation. It is also slightly easier to implement because it does not require rotating and re-measuring for every term like the pure quantum approach would require.

[AdrianPerezSalinas]

Thanks for your notes and your code Stavros, it is now clear to me. The fact that we need the state to compute the hamiltonian is not a problem for me, since we are simulating in the most efficient manner we are able to do.

I just have a question: you are using a hamiltonian TFIM which can be trotterized. What if we want to compute the expected value of a Hamiltonian that has not been defined in the same way? Can I just trotterize any hamiltonian (from a qibo perspective)

[stavros11]

Thanks for checking this.

I just have a question: you are using a hamiltonian TFIM which can be trotterized. What if we want to compute the expected value of a Hamiltonian that has not been defined in the same way? Can I just trotterize any hamiltonian (from a qibo perspective)

Yes, you can trotterize (from the qibo perspective) any Hamiltonian. The TrotterHamiltonian name may be a bit deceptive, I guess we named it this way because we created with Trotter evolution in mind. However, when it comes down to the actual implementation it is just a Hamiltonian that is defined using local terms instead of the full 2^Nx2^N Hamiltonian. It has some methods (such as the .circuit) that are related to Trotter evolution, but if you don't use them you never actually perform any trotterization. The above application is an example of this.

I think the easiest way to define your custom Hamiltonian as TrotterHamiltonian is to write it using sympy symbols and use TrotterHamiltonian.from_symbolic, like the example I referred above. The symbols method is easy to use because you essentially use the Hamiltonian with the form you would write it on paper.