# Qiskit - Approximation of Hamiltonian energy via QPE

I'm trying to study QPE with the motivation of obtaining eigenvalues of Hamiltonian, i.e. energies of a system. My problem is, that while np.linalg.eig and VQE are agreeing on the lowest energy, energy obtained via QPE is completely off.

My Hamiltonian looks like this $$\widehat{H} = -II + 0.3IZ -0.01ZY + 0.1YX$$

and I'd expect its eigenvalues to be (roughly) [-0.39328755, -0.89, -1.11, -1.60671245]. Using VQE I'm obtaining -1.6067123508302064, i.e. the number very close to the lowest eigenvalue.

On the other hand, when I run my QPE code, I'm getting number like -0.392699, not knowing, where is the problem.

### My code

import numpy as np
from qiskit import QuantumCircuit, execute
from qiskit.algorithms.minimum_eigensolvers import VQE
from qiskit.algorithms.optimizers import SLSQP
from qiskit.circuit.library import PhaseEstimation, TwoLocal
from qiskit.extensions import HamiltonianGate
from qiskit.primitives import Estimator
from qiskit.quantum_info import SparsePauliOp
from qiskit_aer import AerSimulator

H = SparsePauliOp.from_list([('II', -1), ('IZ', 0.3), ('XI', -0.3), ('ZY', -0.01), ('YX', 0.1)])
Hmat = H.to_matrix()

eig = np.linalg.eigvals(Hmat)
print(eig)

estimator = Estimator()
optimizer = SLSQP()
ansatz = TwoLocal(rotation_blocks=['ry', 'rz'], entanglement_blocks='cz')

vqe = VQE(estimator, ansatz, optimizer)
result = vqe.compute_minimum_eigenvalue(operator=H)
print(result.eigenvalue)

n_qpe_qbits = 10

U = HamiltonianGate(Hmat, 1, label='H')

# Obtain a solution via QPE
total_qbits = U.num_qubits + n_qpe_qbits
measure_circ = QuantumCircuit(total_qbits, n_qpe_qbits)
qpe = PhaseEstimation(n_qpe_qbits, U)

measure_circ = measure_circ.compose(qpe)
measure_circ.measure(range(n_qpe_qbits), range(n_qpe_qbits))
print(measure_circ.decompose())

backend = AerSimulator(method='statevector')
job = execute(measure_circ, backend, shots=2048)
counts = job.result().get_counts()
print(counts)

max_count = max(counts.items(), key=lambda x: x)
print(f'MAX count: {max_count}')

theta = int(max_count[::-1], 2) / 2**n_qpe_qbits
print(f'Theta value: {theta}')
print(f'QPE-approximated U-eigenvalue: {np.exp(2*1j*np.pi * theta)}')
print(f'QPE-approximated H-eigenvalue: {-2 * np.pi * theta}')
$$$$


Most likely because your initial state has very small overlap with your ground state. You can check by just looking at the eigenvector correspond to the lowest eigenvalue. So change the initial state. I changed it to uniformly superposition state by adding the following line of code right after you defined the measure_circ quantum circuit.

measure_circ.h([-1, -2])


This seems to give back the right result.

MAX count: ('0110000010', 4605)
Theta value: 0.255859375
QPE-approximated U-eigenvalue: (-0.036807222941358866+0.9993223845883495j)
QPE-approximated H-eigenvalue: -1.607611865704152


Also remember that you may have to run QPE for some amount of times before you see the ground state energy (lowest eigenvalue). This depends on the overlap between the ground state and the inital state you prepared.

Code:

import numpy as np
from qiskit import QuantumCircuit, execute
from qiskit.algorithms.minimum_eigensolvers import VQE
from qiskit.algorithms.optimizers import SLSQP
from qiskit.circuit.library import PhaseEstimation, TwoLocal
from qiskit.extensions import HamiltonianGate
from qiskit.primitives import Estimator
from qiskit.quantum_info import SparsePauliOp
from qiskit_aer import AerSimulator

H = SparsePauliOp.from_list([('II', -1), ('IZ', 0.3), ('XI', -0.3), ('ZY', -0.01), ('YX', 0.1)])
Hmat = H.to_matrix()

eig = np.linalg.eigvals(Hmat)
print('NUMPY eig',eig)

estimator = Estimator()
optimizer = SLSQP()
ansatz = TwoLocal(rotation_blocks=['ry', 'rz'], entanglement_blocks='cz')

vqe = VQE(estimator, ansatz, optimizer)
result = vqe.compute_minimum_eigenvalue(operator=H)
print('VQE eig', result.eigenvalue)

n_qpe_qbits = 10

U = HamiltonianGate(Hmat, 1, label='H')

# Obtain a solution via QPE
total_qbits = U.num_qubits + n_qpe_qbits
measure_circ = QuantumCircuit(total_qbits, n_qpe_qbits)
measure_circ.h([-1, -2])

qpe = PhaseEstimation(n_qpe_qbits, U)

measure_circ = measure_circ.compose(qpe)
measure_circ.measure(range(n_qpe_qbits), range(n_qpe_qbits))
print(measure_circ.decompose())

backend = AerSimulator(method='statevector')
job = execute(measure_circ, backend, shots=10000)
counts = job.result().get_counts()
print(counts)

max_count = max(counts.items(), key=lambda x: x)
print(f'MAX count: {max_count}')

theta = int(max_count[::-1], 2) / 2**n_qpe_qbits
print(f'Theta value: {theta}')
print(f'QPE-approximated U-eigenvalue: {np.exp(2*1j*np.pi * theta)}')
print(f'QPE-approximated H-eigenvalue: {-2 * np.pi * theta}')
`