# A question about qiskit's WeightedPauliOperator.evolve method for IPEA (IQPE) implementation

I'm doing an implementation/tutorial of the iterative phase estimation algorithm (IPEA) and have a problem. My Hamiltonian is a diagonal matrix $$H = \begin{pmatrix} E_1 & 0 \\ 0 & E_2\end{pmatrix}$$, I want to estimate both eigenvalues. When I do it instead of $$E_1$$ the estimated value is 0, and instead of $$E_2$$ the estimated value is $$E_2 - E_1$$. Is it possible that when I call H.evolve(t=1) method it doesn't create a circuit for $$\mathrm{e}^{-iHt}$$, but instead, creates a circuit for $$\mathrm{e}^{-i\hat{H}t}$$ , where $$\hat{H} = \begin{pmatrix}0 & 0 \\ 0 & E_2 - E_1\end{pmatrix}$$?

In other words, if I have $$H = \begin{pmatrix} E_1 & 0 \\ 0 & E_2\end{pmatrix}$$, what should do the H.evolve(t=1) method?

a) Create a circuit for $$\mathrm{e}^{-iHt}$$ unitary operator.

b) Create a circuit for $$\mathrm{e}^{-i\hat{H}t}$$ unitary operator.

This is the link to my tutorial and the short version of it is presented here:

In the tutorial I am writing

E_1, E_2 = (0, random())


and it finds $$E_2$$, because when $$E_1 = 0$$, $$E_2 = E_2 - E_1$$.

Here in the following code, this line is changed:

E_1, E_2 = (random(), random())


and the algorithm estimates $$E_2 - E_1$$ value instead of $$E_2$$.

import numpy as np
from random import random
from qiskit import *
from qiskit.aqua.operators import WeightedPauliOperator, MatrixOperator
from qiskit.aqua.operators.op_converter import to_matrix_operator
from qiskit.aqua.utils.controlled_circuit import get_controlled_circuit

backend = BasicAer.get_backend('qasm_simulator')
q = QuantumRegister(1)
a = QuantumRegister(1)
c = ClassicalRegister(1)

def simple_hamiltonian_operator(E_1, E_2):
"""
Creates 0.5*(E_1 + E_2)*I + 0.5*(E_1 - E_2)*Z pauli sum
that will be our "simple" Hamiltonian operator. The corresponding
matrix for Hamiltonian is [[E_1, 0], [0, E_2]].

"""
pauli_dict = {
'paulis': [{"coeff": {"imag": 0.0, "real": 0.5 * (E_1 - E_2)}, "label": "Z"},
{"coeff": {"imag": 0.0, "real": 0.5 * (E_1 + E_2)}, "label": "I"}
]
}
return WeightedPauliOperator.from_dict(pauli_dict)

E_1, E_2 = (random(), random())
print("We want to estimate E_2 = {}".format(E_2))

H = simple_hamiltonian_operator(E_1, E_2)
print("The Hamiltonian in matrix form:")
print(to_matrix_operator(H).dense_matrix)

t = 1
H_circuit = H.evolve(evo_time=t, quantum_registers=q)
H_circuit.data.__delitem__(-1) # deleting a barrier at the end of circuit

# control version of the circuit
control_H = get_controlled_circuit(H_circuit, a[0])

num_bits_estimate = 10
phase = 0
for k_precision in reversed(range(num_bits_estimate)):
# Create a Quantum Circuit acting on the q register
k_circ = QuantumCircuit(q, a, c)

# (1) |1> eigenstate initialization
k_circ.x(q[0])

# (2) Initial Hadamard gate applied on ancillary qubit.
k_circ.h(a[0])

# (3) The control Hamiltonian applied on the qubits where control qubit is the ancillary qubit.
for order in range(2 ** k_precision):
k_circ += control_H

# (4) The phase gate and final Hadamard gate on ancillary qubit.
phase_shift = 2 * np.pi * phase * 2 ** k_precision
k_circ.u1(-phase_shift, a[0])
k_circ.h(a[0])

# (5) Measurement of ancillary qubit (findig the bit)
k_circ.measure(a[0], c[0])

# (6) executing on Quantum Computer and finding a bit from the phase
job = execute(k_circ, backend, shots=1000)
result = job.result()
counts = result.get_counts()
value = int(max(counts, key=counts.get))

# (7) phase after this iteration
phase += value / 2 ** (k_precision + 1)

# obtaining E_2 from the estimated phase
if E_2 > E_1:
eigenvalue = 2 * np.pi * (1 - phase) / t
else:
eigenvalue = -2 * np.pi * phase / t

print("Eigenvalue of Hamiltonian that we want to estimate: E_2 = " + str(E_2))
print("Meanwhile, E_2 - E_1 = {}".format(E_2 - E_1))
print("Estimated eigenvalue of Hamiltonian: " + str(eigenvalue))