Quantum Phase Estimation (QPE) with ProjectQ

by Fernando de la Iglesia

Quantum Phase Estimation is one of the most relevant algorithms in quantum computing, whose importance resides in that it is used as part of other more complex algorithms. The objective of this algorithm is to estimate the phase, and that means the value, of the eigenvalue corresponding to the eigenvector of a unitary operator.

In short, with this algorithm we can evaluate θ in the expression

For a quick reference you can visit the corresponding Wikipedia article, were you can find more references therein. Another important reference is the great paper by Daniel S. Abrams and Seth Lloyd: arXiv:quant-ph/9807070.

Recently this algorithm was added as a native operation in ProjectQ. The aim of this post is to show how to use this native operation.

Please refer to the ProjectQ documentation for setting up an environment to run the framework and to be able to run programs that make use of the framework. This post assumes that you have already installed all the required libraries. In addition, here you can find a link to a jupyter notebook including the examples presented below that include all the required installations and library imports.

To use the Quantum Phase Estimation native operation you have simply to import it in your program.

from projectq.ops import QPE

Provided you import as well all the required libraries typically you need to run ProjectQ programs.

The QPE operation acts on the qubits that form the system (system_qubits in the examples below) and on the ancillary qubits that at the end of the process store the estimated phase (qpe_ancillas in the examples). As it is described in the references provided above, the amount of ancillary qubits determine the precision of the phase estimation. In the first example, we will use a simple configuration with a single system qubit and three ancillas

n_qpe_ancillas = 3
qpe_ancillas = eng.allocate_qureg(n_qpe_ancillas)                          system_qubits = eng.allocate_qureg(1)

The other important component is the Unitary operation (U in the expression above) for which we want to obtain the eigenvalues corresponding to the eigenvectors. For this first example, we will use a simple Phase gate that adds a phase (0.125) to the initial state. The Phase gate is included as a native operation in ProjectQ

from projectq.ops import Ph
angle = math.pi*2.*0.125
U = Ph(angle)

With this we can now estimate the phase for the Unitary Ph and the state initialized as | ψ> = |0>

# Apply Quantum Phase Estimation
QPE(U) | (qpe_ancillas, system_qubits)                                             All(Measure) | qpe_ancillas
# Compute the phase from the ancilla measurement                          #(https://en.wikipedia.org/wiki/Quantum_phase_estimation_algorithm)                          phasebinlist = [int(q) for q in qpe_ancillas]                          phase_in_bin = ''.join(str(j) for j in phasebinlist)                          phase_int = int(phase_in_bin,2)
phase = phase_int / (2 ** n_qpe_ancillas)                          print (phase)
All(Measure) | system_qubits
eng.flush()

The result we find (as expected because we configured the Unitary in that way) is that the phase is 0.125.

A more elaborated example would use in one side a complex Unitary in the form of a class that inherits the Basic Gate class that allows being more general in terms of the unitary operator for which we would like to estimate the phase. And in the other side a more general state that is not in the form of an eigenvector for the unitary operation, but that at the end of the day can be expressed as a superposition of eigenvectors of the Unitary operator (I recommend to review the paper form Abrams and Lloyd at this point).

Let us clarify the ideas using a simple example.

First, we define a unitary operator whose eigenvectors will be described below

This operator has eigenvectors

With eigenvalues 1 and -1 respectively, or equivalently, phases 0.0 and 0.5.

(Please be aware of how ProjectQ manage the order of qubits in multiple qubits states, the operator and states above are taking this ordering into account to be consistent with the notation of the code below, that is (00,10,01,11)).

We could have a state that can be expressed as a superposition of these eigenvectors

That in the standard base looks as

Therefore, we can use our QPE operation in ProjectQ to estimate the phase for the eigenvectors of the defined operator included in a superposition of them. Let us review the required code.

First, we define the operator as an extension of BasicGate as stated before.

from projectq.ops import BasicGate
# Definition of the elaborated unitary operation
class unitary_operation(BasicGate):
  def __init__(self):
    BasicGate.__init__(self)
  
  @property
  def matrix(self):
     return np.matrix([[0, 1, 0, 0],
                       [1, 0, 0, 0],
                       [0, 0, 0, 1],
                       [0, 0, 1, 0]])
unitary_op = unitary_operation()

(Edited in June 30, 2019): I found that the correct (and simplest) way to create a gate or operator as a matrix is the following:

from projectq.ops import MatrixGate
unitary_op = MatrixGate([[0, 1, 0, 0],
                         [1, 0, 0, 0],
                         [0, 0, 0, 1],
                         [0, 0, 1, 0]])

Extension of BasicGate should not be used)

It follows our qubits and ancilla, the creation of the initial state, as well as the instantiation of the operator.

(As stated above, please be aware of how ProjectQ manage the order of qubits in multiple qubits states, the creation of the state with StatePreparation is taking this ordering into account to be consistent with the notation of the code below and the definition of the operator).

eng2 = MainEngine()
n_qpe_ancillas = 1
system_qubits = eng2.allocate_qureg(2)
qpe_ancillas = eng2.allocate_qureg(n_qpe_ancillas)
# Create the state
amplitude00 = (np.sqrt(3) + 1.)/(np.sqrt(2) * 2.)
amplitude10 = (np.sqrt(3) - 1.)/(np.sqrt(2) * 2.)
StatePreparation([amplitude00, amplitude10, 0., 0.]) | system_qubits

unit = unitary_op

And now the application of the QPE algorithm.

# Apply Quantum Phase Estimation
QPE(unit) | (qpe_ancillas, system_qubits)
All(Measure) | qpe_ancillas
fasebinlist = [int(q) for q in qpe_ancillas]
fasebin = ''.join(str(j) for j in fasebinlist)
faseint = int(fasebin, 2)
phase = faseint / (2. ** (len(qpe_ancillas)))
results = phase

The estimated phase stored in results could belong to either eigenvector |+> or |-> (see above). In order to be able to identify the corresponding eigenvector we need first to transform the final state to the eigenvectors’ base using in this case a simple Hadamard gate applied to the final state.

H | system_qubits[0]

Now, depending on the estimated phase we can verify to which eigenvector corresponds

if np.allclose(phase, .0, rtol=1e-1):
  results_plus = phase
  All(Measure) | system_qubits
  autovector_result = [int(q) for q in system_qubits]
  print (autovector_result, phase)
elif np.allclose(phase, .5, rtol=1e-1):
  results_minus = phase
  All(Measure) | system_qubits
  autovector_result = [int(q) for q in system_qubits]
  print (autovector_result, phase)
eng2.flush()

In accordance with the probabilities of obtaining each of the eigenstates when running this scripts, we will find 0.0 as the phase for state 00 (remember that we have applied H to the resulting state) and 0.5 for state 10. These probabilities depend on superposition in the original state and on the probability of finding the correct phase as per the QPE algorithm.

In addition QPE can be used to estimate the phase for a time parameter defined function so that

Let us clarify the idea with an example. First, we define the function.

def two_qubit_gate(system_q, time):
  CNOT | (system_q[0], system_q[1])
  Ph(2.0*cmath.pi*(time * 0.125)) | system_q[1]
  CNOT | (system_q[0], system_q[1])

As you can easily verify, when we exponentiate this function the result is equivalent to set the parameter time to the exponent. ProjectQ QPE is able to do phase estimation with functions that comply with this requirement.

We can estimate the phase for the eigenvector |10> of the previous function in the same way we did previously.

eng3 = MainEngine()
n_qpe_ancillas = 3
qpe_ancillas = eng3.allocate_qureg(n_qpe_ancillas)
system_qubits = eng3.allocate_qureg(2)
X | system_qubits[0]
# Apply Quantum Phase Estimation
QPE(two_qubit_gate) | (qpe_ancillas, system_qubits)
All(Measure) | qpe_ancillas
phasebinlist = [int(q) for q in qpe_ancillas]
phase_in_bin = ''.join(str(j) for j in phasebinlist)
phase_int = int(phase_in_bin,2)
phase = phase_int / (2 ** n_qpe_ancillas)
print (phase)
All(Measure) | system_qubits
eng3.flush()

The estimated phase as expected is 0.125 with high probability.