Iterative Quantum Amplitude Estimation

What Are We Looking For?

IQAE, like QAE, is specifically focused on analyzing data sets composed of “good” and “bad” states. We can assign each of the subspaces a marker that indicates which of the two categories it falls in. Taking \(|0\rangle\) to be a “bad” marker and \(|1\rangle\) to be a “good” marker, the general IQAE state can be defined as

where \(a\) is the probability of measuring a “good” state, \(|\psi_0\rangle\) is a “bad” state, and \(|\psi_1\rangle\) is a “good” state.

In this implementation, the goal of the IQAE algorithm will be to estimate the fraction of states that are multiples of 8. When encoded in binary, multiples of 8 will always have 0 in the three least significant positions, which serves as a simple success criterion for the algorithm (i.e. a “good” state will always have 0s in the 3 least significant positions).

To carry this search out, we define an operator \(\mathcal{A}\) satisfying \(\mathcal{A}|0\rangle^{\otimes n+1}=|\Psi_{IQAE}\rangle\), where \(n+1\) is the number of qubits in the system, including the evaluation qubit. This operator is determined by the application.

In this case, \(\mathcal{A}\) prepares a random state and flips an auxiliary qubit via a multi-controlled CNOT gate if the 3 least significant bits adhere to the success criterion (i.e. are they all zero?). If the logic gate is triggered, this single-qubit evaluation register will be flipped to \(|1\rangle\), indicating a “good” result.

In IQAE, the goal is not to estimate the number of “good” results in one sweep. Instead, several iterations will be carried out in which the Grover operator is applied \(k\) times with the goal of extracting a probability amplitude by refining the interval within which the solution is likely to lie. In doing this, we avoid the expensive QPE procedure used in QAE. To do this, each iteration will yield the general result

where \(\mathcal{Q}\) is the Grover operator and \(\theta_a\) is related to the angle by which the state is rotated by Grover’s operator (note \(a=\sin^2(\theta_a)\)), and \(k\) is the number of times that the Grover operator is applied to the state in a single IQAE iteration. The specifications of this equation are covered thoroughly in [2], but the important result is that the probability of measuring a “good” state at the end of an iteration is given by

From this, it is clear that the probability is correlated to the angle imposed by the Grover operator. This relationship can be made clearer by taking into consideration the fact that \(\mathcal{Q}\) invokes a \(2\theta_a\) rotation in a 2D space defined by the “good” and “bad” state axes each time it is applied. So, if we can figure out the angle that has been imposed on the state as a result of the operator application, we can obtain the probability of extracting a “good” state. Since we do not aspire to use QPE, our best bet is to use our iterative \(k\) value combined with a measurement of the quantum circuit that outputs an intermediate amplitude guess.

The above equation also shows that the size of \(k\) determines the resolution of the search, with a large \(k\) corresponding to a high probability oscillation frequency and, therefore, a high resolution. If the frequency becomes too high, however, the measurement outcomes can become ambiguous, once again justifying the iterative approach to ensure an ideal \(k\) value is found using sequential information collected from the system [1].

Defining The Input State and Operators

First, we can define the circuit specifications. For this toy example, we will generate a random list of probabilities that will be assigned as weights in the input state.

import pennylane as qp
import numpy as np
import math
import jax.numpy as jnp
import time

import catalyst
from catalyst import for_loop

#Define system parameters
shots = 10000 #Number of shots
num_qubits = 5
N = 2**num_qubits #Number of possible states

#Generate a random list of probabilities to be assigned to the initial state
random_vector = np.sqrt(np.random.rand(N))
distribution = random_vector/np.linalg.norm(random_vector)

#Define which indices should be checked for success criteria
MCX_wires = [num_qubits-3,num_qubits-2,num_qubits-1,num_qubits]

As mentioned, the backbone of the quantum portion of the IQAE algorithm is the Grover operator \(\mathcal{Q}\), which aims to identify “good” states and introduce an identifiable phase flip and amplitude amplification. The basic structure of \(\mathcal{Q}\) is

In which \(\mathcal{S}_{\psi_1}\) acts as the oracle and flips the phase of (marks) a “good” state and \(\mathcal{S}_0\) should flip the phase of the \(|0\rangle\) state (see figures 1 through 3 in the amplitude amplification demo for an intuitive visualization). Since this is a non-uniform superposition, the operator that facilitates this process needs to be defined rather than using PennyLane’s built-in GroverOperator function, which assumes a uniform superposition.

First, \(\mathcal{A}\) can be defined according to the following procedure:

1. Generate \(n\) qubits with amplitudes according to the previously generated random probability distribution using StatePrep().

2. Flip the state of the 3 final qubits in the string so that MultiControlledX() is triggered by a \(|111\rangle\) state.

3. Implement MultiControlledX() such that wire \(n+1\) takes on the \(|1\rangle\) state if the success criterion is met.

4. Flip the state of the 3 final qubits back to the original.

#Define A operator
@qp.prod
def A(state):
  qp.StatePrep(state,wires=range(num_qubits)) #Randomly weighted superposition

  #Flip monitored qubits so that MCX is triggered by |111>
  qp.PauliX(wires=num_qubits-3)
  qp.PauliX(wires=num_qubits-2)
  qp.PauliX(wires=num_qubits-1)

  qp.MultiControlledX(MCX_wires) #State is marked with |1> iff number is a multiple of 8

  #Flip back monitored bits to original state
  qp.PauliX(wires=num_qubits-3)
  qp.PauliX(wires=num_qubits-2)
  qp.PauliX(wires=num_qubits-1)

Since the “good” state is marked by a \(|1\rangle_{eval}\), the \(\mathcal{S}_{\psi_1}\) operator (also known as the oracle) can be constructed simply using a PauliZ() gate, which will flip the phase of any state that has this marker and allow any state marked by \(|0\rangle_{eval}\) to pass unchanged. The \(\mathcal{S}_0\) operator is analogous to a simple FlipSign() operation defined to act on the \(|0\rangle\) state. When designing the operator, the global phase term can be neglected since it does not impact the final measurement in IQAE. It should be noted that this is not the case in canonical QAE, where global phase does make a difference in the final result [2].

These operators can be used to build the final, iterative circuit in which the number of Grover operator applications will vary per iteration.

Due to the repeated applications of the Grover operator (the exact number of which, as will be explored soon, grows quickly between iterations), the computational demand of this algorithm can become quickly unmanageable. To mitigate this, PennyLane’s Catalyst compiler can be used to compile the \(\mathcal{Q}\) loop and reduce the demand. For small systems (like, for example, 5 qubits), the benefit is negligible but becomes more apparent as the system grows.

k_i = 0
dev = qp.device("lightning.qubit", wires=num_qubits+1)

#Declare if Catalyst compiler should be used
catalyst_bool = True

#Build the circuit Q^kA|0>n|0>
def circuit_builder(catalyst_bool):
  if catalyst_bool:
    @qp.qnode(dev, shots=shots)
    def circuit(state, k_i):
      A(state)
      k_int = jnp.int64(k_i)
      @catalyst.for_loop(0, k_int, 1)
      def apply_Q(k_int):
        #Build the circuit using the Grover operator form Q=AS0A*Spsi
        qp.PauliZ(wires=num_qubits)
        qp.adjoint(A(state))
        qp.FlipSign(0, wires=range(num_qubits + 1))
        A(state)
      apply_Q()
      return qp.probs(wires=[num_qubits])
    circuit = catalyst.qjit(circuit)

  else:
    @qp.qnode(dev, shots=shots, interface=None)
    def circuit(state, k_i):
      A(state)
      for i in range(int(k_i)):
        #Build the circuit using the Grover operator form Q=AS0A*Spsi
        qp.PauliZ(wires=num_qubits)
        qp.adjoint(A(state))
        qp.FlipSign(0, wires=range(num_qubits + 1))
        A(state)
      return qp.probs(wires=[num_qubits])

  return circuit

circuit = circuit_builder(catalyst_bool)

Digesting the FindNextK Function

As shown, the iteration variable \(k\) is directly tied to the total angle of the state (which can be achieved via measurement of the quantum circuit) since the Grover operator invokes a deterministic rotation each time it is applied. The \(\sin^2(x)\) function adds complexity to the probability calculations, so standard trigonometric identities can be employed to achieve

Letting \(K_i=4k+2\) be the frequency.

In [1], the authors define FindNextK() to determine the number of times \(\mathcal{Q}\) is implemented per iteration. The goal of FindNextK() is to identify the largest possible \(k\) that adheres to what is called the half-plane condition. The core principle of IQAE is the narrowing of a range of potential amplitudes to, eventually, home in on an accurate estimate of the “good” state probability. To do this, each iteration of the algorithm must operate between an upper and lower bound that make up the function’s confidence interval, which corresponds to the range of probabilities within which the final probability amplitude will exist. Since the calculation involves a \(\arcsin(x)\) calculation, it is possible that this confidence interval could yield uninterpretable results if one angle falls in the upper half of the unit circle (i.e. between 0 and \(\pi\)) and the other falls in the lower half of the unit circle (i.e. between \(\pi\) and \(2\pi\)) since information about the measurement’s position on the probability curve would be lost. So, having both bounds of the confidence interval on the same half-plane of the unit circle will result in unambiguous knowledge on which branch of \(\arcsin(x)\) should be used. Thus, valid results should always fall in either the upper or lower half-plane of the unit circle.

The FindNextK function validates this condition. The logic is as follows:

1. For an initial guess \(k_i\) yielding confidence interval \([\theta_{min}^i,\theta_{max}^i]=[\theta_{lower}K_i,\theta_{upper}K_i]\), the function will return the current guess of \(k\) if either both the upper and lower bounds are less than \(\pi\) (i.e., they fall in the upper half of the unit circle) or both the upper and lower bounds are greater than \(\pi\) (i.e. they fall in the lower half of the unit circle).

2. If neither of these conditions are met (i.e., the two bounds fall in different half-planes), the magnitude of the guess needs to be reduced.

To carry out the actual comparison logic, however, some translation is required. First, the maximum possible value of \(k\) must be defined in relation to the available angles. [1] defines this value as

Where \(K_{max}\) can be interpreted as the maximum number of rotations that can be carried out before aliasing becomes an issue.

So, our goal is to find the largest integer \(K \leq K_{max}\) that both satisfies \(K = 4k+2\) and adheres to the half-plane condition. This can be carried out using a modulo 4 calculation, which enforces the required condition introduced by \(K=4k+2\). Once this is found, the final step is to compute scaling factor \(q\), the ratio between the current \(K\) guess and the previous, which shifts the values relative to the previous step to prevent backsliding.

The FindNextK function will achieve one of two outcomes.

1. Both bounds of the confidence interval fall in the same half-plane, causing the function to return the current \(k\) guess and a Boolean indicating which half-plane the interval falls in.

or

2. The bounds of the confidence interval fall in different half-planes, indicating the current guess is too large. If an adequate guess is not reached as the while loop runs, the previous guess is returned.

def FindNextK(k_i,theta_min, theta_max, HalfPlane_bool):
    K_i = 4*k_i+2 #Define coefficient
    theta_min_i = theta_min*K_i
    theta_max_i = theta_max*K_i
    Kmax = math.floor(math.pi/(theta_max-theta_min)) #Maximum K value
    K = Kmax-(Kmax-2)%4

    while K>=2*K_i:
      q = K/K_i
      if (q*theta_max_i)%(2*math.pi)<=math.pi and (q*theta_min_i)%(2*math.pi)<=math.pi:
        k_i_it = ((K-2)/4)
        HalfPlane_Bool_it = True
        return k_i_it, HalfPlane_Bool_it

      elif (q*theta_max_i)%(2*math.pi)>=math.pi and (q*theta_min_i)%(2*math.pi)>=math.pi:
        k_i_it = ((K-2)/4)
        HalfPlane_Bool_it = False
        return k_i_it, HalfPlane_Bool_it

      else:
        K-=4 #Decrease guess if the range spans two halves
    return (k_i,HalfPlane_bool)

Implementing the IQAE Algorithm

With FindNextK() defined, the IQAE algorithm can now be implemented! The main objective of this function is to apply the \(k\) value returned by FindNextK() to the previously defined quantum circuit, obtain a measurement, and determine if this measurement is adequately accurate with respect to a defined tolerance or if the confidence interval should be updated and passed back into the classical function for another iteration. The logic is as follows:

1. Call circuit() after FindNextK() outputs a guess for \(k\) and take a probability measurement.

2. Use this value to update the confidence interval, in which both the upper and lower bounds on the angles and probabilities are computed from the measured amplitude.

3. Compute the overlap between the previous confidence interval and the new confidence interval, taking this to be your final upper and lower bound definition.

4. Check to see if the difference between the new upper and lower bounds is smaller than \(\epsilon\), which represents a chosen accuracy parameter. If not, pass the final upper and lower bounds back into FindNextK() and repeat. If yes, return the probability amplitudes associated with the upper and lower amplitudes.

The confidence intervals can be determined using the Chernoff-Hoeffding bound, where the margin of error is given by \(\epsilon_{a_i}\) from [1].

Where \(\epsilon_{a_i}\) is the margin of error and \(T\) defines the maximum number of iterations required to achieve a precision of \(\epsilon_{a_i}\) and

We can then estimate the upper and lower bounds of the probability interval estimate.

In which \(a_i\) is the outcome of the quantum circuit measurement for iteration \(i\).

#Define IQAE parameters
eps = 0.0001 #Precision
alpha = 0.01 #Confidence

#Actually implement IQAE!
#Pre-selecting the use of Chernoff-Hoeffding to determine confidence interval
def IQAE(eps, alpha, shots):
  k_current = 0
  theta_lower = 0
  theta_upper = math.pi/2 #Begin search in the first quadrant
  HalfPlane_Bool = True
  T = math.ceil(math.log2(math.pi/(8*eps)))

  while (theta_upper-theta_lower)>2*eps:
    k_i, HalfPlane_Bool = FindNextK(k_current, theta_lower, theta_upper, HalfPlane_Bool) #determine current guess of k_i
    K_i = int(4*k_i+2)

    #Call circuit
    a_estimate = (circuit(distribution,k_i)[1])

    eps_ai = ((1/(2*shots))*math.log(2*T/alpha))**0.5

    p_max = (np.clip(a_estimate+eps_ai, 0, 1))
    p_min = (np.clip(a_estimate-eps_ai, 0, 1))

    #Update confidence interval based on the quadrant we are in
    if HalfPlane_Bool: #Upper quadrant
      theta_upper_est = 2*math.asin((p_max)**0.5)
      theta_lower_est = 2*math.asin((p_min)**0.5)
    else: #Lower quadrant - ensure the value stays within [0,2pi]
      theta_upper_est = 2*math.pi - 2*math.asin((p_min)**0.5)
      theta_lower_est = 2*math.pi - 2*math.asin((p_max)**0.5)

    #Compute total rotation
    new_lower = (2*math.pi*(math.floor(K_i*theta_lower/(2*math.pi))) + theta_lower_est)/K_i
    new_upper = (2*math.pi*(math.floor(K_i*theta_upper/(2*math.pi))) + theta_upper_est)/K_i

    #Intersection calculation
    theta_lower = max(theta_lower, new_lower)
    theta_upper = min(theta_upper, new_upper)

    k_current = k_i

  a_lower = math.sin(theta_lower)**2
  a_upper = math.sin(theta_upper)**2

  return a_lower, a_upper

For the purposes of this demonstration and to emphasize the bare bones of the IQAE algorithm, a few aspects of the full implementation presented in [1] were omitted. For example, the explicit overshooting condition was not translated here since it is not necessary to achieve the search outcome but required to obtain the performance guarantees (e.g., analytical bounds) derived by Grinko et al. Full exploration of these additional components can be found in the source paper.

Upon calling IQAE(), the output will consist of the upper and lower bounds between which the true amplitude lies.

dev_exact = qp.device("default.qubit", wires=num_qubits+1)

@qp.qnode(dev_exact)
def circuit_exact():
    A(distribution)
    return qp.probs(wires=[num_qubits])

#Define parameters
eps = 0.0001 #Precision
alpha = 0.01 #Confidence
HalfPlane_Bool = True

t0 = time.time()
a_lower, a_upper = IQAE(eps, alpha, shots)
t1 = time.time()
true_a = circuit_exact()[1]

#Read results
if catalyst_bool:
  print("Compiled with Catalyst")
print("Execution time:", t1-t0,"s")
print("Analytic probability:", true_a)
print("Lower prediction bound:", a_lower)
print("Upper prediction bound:", a_upper)
print("Contains true value?", a_lower <= true_a <= a_upper) #Boolean indicating whether the analytical value falls within the final confidence interval
print("-------------------------")
if a_lower <= true_a <= a_upper:
  print("The probability of measuring a multiple of 8 falls between", (100*a_lower),"%", "and", (100*a_upper),"%")
else:
  print("No valid answer found")

Compiled with Catalyst
Execution time: 10.060495376586914 s
Analytic probability: 0.17582667305655825
Lower prediction bound: 0.17578353772354324
Upper prediction bound: 0.17584602263711024
Contains true value? True
-------------------------
The probability of measuring a multiple of 8 falls between 17.578353772354323 % and 17.584602263711023 %

Conclusion

QAE was put forward as a generalized adaptation of quantum search methods, which hinted at the viability of quantum computers as a notable improvement on existing technologies early on. IQAE pushed forward with this task, showing not only that quantum hardware has the potential to provide advantages in estimation algorithms but that classical feedback loops can be employed to reduce implementation limitations. The implementation shown here does not achieve full quadratic speedup compared to classical estimation methods, but advancements have been subsequently made to meet this metric [3].

As alluded to, IQAE has the potential to make notable impacts in various applications, such as in the calculation of expectation values in applications like quantum chemistry. Proposals have also been put forward to employ IQAE in areas such as risk analysis, financial portfolio modelling and optimization, and power grid analysis. The true feasibility of using IQAE (and, really, any other quantum algorithm) in these applications will become clearer as quantum hardware advances in capability and usability.

There are, of course, limitations to IQAE’s feasibility, particularly that the circuit depth grows with the number of iterations. Despite this, IQAE shows an example of how quantum-classical hybrid approaches can be used to achieve speedup compared to classical methods and resource requirement reduction compared to other quantum methods.

References