I have followed this tutorial and based on it, I've written the following function in qiskit, which can explicitly calculate the ground states of a transverse-field Ising Hamiltonian.
from qiskit import *
import numpy as np
def Hamiltonian(n,h):
pow_n=2**n
qc = np.empty(2*n-1, dtype=object)
#Creating the quantum circuits that are used in the calculation of the Hamiltonian based on the number of qubits
for i in range(0, 2*n-1): #2n-1 is the number of factors on the n-site Hamiltonian
qr = QuantumRegister(n)
qc[i] = QuantumCircuit(qr) #create quantum circuits for each factor of the Hamiltonian
#print(i)
if (i<=n-2): #for the first sum of the Hamiltonian
qc[i].z(i) #value of current spin
qc[i].z(i+1) #and value of its neighboring spin
else: #for the second sum of the Hamiltonian
qc[i].x(2*n-2-i) #2*n-2 gives the proper index since counting starts at 0
#Run each circuit in the simulator
simulator = Aer.get_backend('unitary_simulator')
result = np.empty(2*n-1, dtype=object)
unitary = np.empty(2*n-1, dtype=object)
Hamiltonian_Matrix=0
#Get the results for each circuit in unitary form
for i in range(0, 2*n-1):
result[i] = execute(qc[i], backend=simulator).result()
unitary[i] = result[i].get_unitary()
#print(unitary[i])
#And calculate the Hamiltonian matrix according to the formula
if (i<=n-2):
Hamiltonian_Matrix=np.add(Hamiltonian_Matrix,-unitary[i])
else:
Hamiltonian_Matrix=np.add(Hamiltonian_Matrix,-h*unitary[i])
print("The",pow_n,"x",pow_n, "Hamiltonian Matrix is:")
print(Hamiltonian_Matrix)
#Now that we have the Hamiltonian
#find the eigenvalues and eigenvectors
w, v = np.linalg.eig(Hamiltonian_Matrix)
print("Eigenvectors")
print(v)
print("Eigenvalues")
print(w)
minimum=w[0]
min_spot=0
for i in range(1, pow_n):
if w[i]<minimum:
min_spot=i
minimum=w[i]
print(min_spot)
groundstate = v[:,min_spot]
#the probability to measure each basic state of n qubits
probability = np.square(groundstate).real
print("The probability for each of the",pow_n,"base states is:")
print(probability)
print("The probabilities for each of the",pow_n,"base states add up to:")
print ("%.2f" % np.sum(probability))
My problem with this piece of code I've written is that it can only run on a unitary simulator. To my understanding (which may lack some of the underlying physics), the Hamiltonian itself is not a "purely" quantum calculation, since there are additions to be made which cannot be expressed with a quantum (unitary) gate, and this is why the resulting Hamiltonian matrix is also not unitary. For example, if you run Hamiltonian(3, 1), the Hamiltonian matrix is:
[[-2.+0.j -1.+0.j -1.+0.j 0.+0.j -1.+0.j 0.+0.j 0.+0.j 0.+0.j]
[-1.+0.j 0.+0.j 0.+0.j -1.+0.j 0.+0.j -1.+0.j 0.+0.j 0.+0.j]
[-1.+0.j 0.+0.j 2.+0.j -1.+0.j 0.+0.j 0.+0.j -1.+0.j 0.+0.j]
[ 0.+0.j -1.+0.j -1.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j -1.+0.j]
[-1.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j -1.+0.j -1.+0.j 0.+0.j]
[ 0.+0.j -1.+0.j 0.+0.j 0.+0.j -1.+0.j 2.+0.j 0.+0.j -1.+0.j]
[ 0.+0.j 0.+0.j -1.+0.j 0.+0.j -1.+0.j 0.+0.j 0.+0.j -1.+0.j]
[ 0.+0.j 0.+0.j 0.+0.j -1.+0.j 0.+0.j -1.+0.j -1.+0.j -2.+0.j]]
Does this mean that there is no way for this approach to run on a real quantum computer where all you can do is measurements on the qubits? I've seen different approaches online such as QAOA or the use of transformations, but I thought if it's so easy to do it with unitaries and some additions, there should be a way to do it with measurements as well.