Skip to main content

View on GitHub

Open this notebook in GitHub to run it yourself
The quantum counting algorithm [1] efficiently estimates the number of valid solutions to a search problem, based on the amplitude estimation algorithm. It demonstrates a quadratic improvement with regard to a classical algorithm with black box oracle access to the function ff. More precisely, given a Boolean function f:{0,1}n→{0,1}f :\{0, 1\}^n\rightarrow\{0,1\}, the counting problem estimates the number of inputs xx to ff such that f(x)=1f(x)=1. This tutorial demonstrates how to estimate the counting problem using a specific variant of the amplitude estimation algorithm: the Iterative Quantum Amplitude Estimation (IQAE) [2]. The IQAE does not rely on the Quantum Phase Estimation algorithm [3], but purely on applications of the grover operator: Q≡−AS0A†Sψ1,Q\equiv - A S_0 A^{\dagger} S_{\psi_1}, thereby reducing the required number of qubits and gates of the circuit, at the expense of additional multiplicative factor polylogarithmic in the error ϵ\epsilon.

​
Setting Up the Problem

We choose this equation: (a+b)<=2(a + b) <= 2 where aa, bb are 2-bit unsigned integers. This equation has six solutions. The goal is to estimate the number of valid solutions out of the 16 possible inputs, with precision 0.50.5.

​
Amplitude Estimation Using Phase Estimation

We first show how to use quantum phase estimation algorithm for quantum counting [3], then solve it using the IQAE method. Given a state ∣ψ⟩|\psi\rangle such that ∣ψ⟩=a∣ψ1⟩+1−a∣ψ0⟩|\psi\rangle=\sqrt{a}|\psi_1\rangle+\sqrt{1-a}|\psi_0\rangle we can measure aa up to arbitrary precision, given the following building blocks:
  1. State preparation: A unitary AA such that: A∣0⟩=∣ψ⟩=a∣ψ1⟩+1−a∣ψ0⟩A|0\rangle = |\psi\rangle = \sqrt{a}|\psi_1\rangle+\sqrt{1-a}|\psi_0\rangle.
  2. Oracle: A unitary Sψ1S_{\psi_1} such that Sψ1=I−2∣ψ1⟩⟨ψ1∣S_{\psi_1}=I-2|\psi_1\rangle\langle\psi_1|, which adds a (−1)(-1) phase to ∣ψ1∣ψ⟩⟩|\psi_1|\psi\rangle\rangle and does nothing to any orthognal states to ∣ψ1⟩|\psi_1\rangle.
This is effectively a reflection around the “good” state ∣ψ1⟩|\psi_1\rangle. Given these two functions, we can construct the Grover operator: Q≡−AS0A†Sψ1,Q\equiv - A S_0 A^{\dagger} S_{\psi_1} , which is exactly the same operator as for the Grover’s search algorithm. In the subspace spanned by ∣ψ1⟩|\psi_1\rangle and ∣ψ0⟩|\psi_0\rangle, QQ has two eigenvalues: λ±=exp⁡(±i2πθ),sin⁡2(πθ)≡a.\lambda_{\pm}=\exp\left(\pm i2\pi \theta \right), \qquad \sin^2 \left(\pi \theta\right)\equiv a. Therefore, if we apply a QPE on A∣0⟩A|0\rangle we have these two eigenvalues encoded in the QPE register; however, both give the value of aa, so there is no ambiguity.

​

​
Arithmetic Oracle

We define the Sψ1S_{\psi_1} oracle: Sψ1∣a⟩∣b⟩=(−1)f(a,b)∣a⟩∣b⟩.S_{\psi_1}|a\rangle|b\rangle= (-1)^{f(a,b)}|a\rangle|b\rangle. 
from classiq import *

A_SIZE = 2
B_SIZE = 2
DOMAIN_SIZE = A_SIZE + B_SIZE


class OracleVars(QStruct):
    a: QNum[A_SIZE]
    b: QNum[B_SIZE]


@qperm
def arith_equation(state: Const[OracleVars], res: QBit):
    res ^= state.a + state.b <= 2


# use phase kickback for turning the arith_equation to an oracle
@qfunc
def arith_oracle(state: OracleVars):
    phase_oracle(arith_equation, state)

​

​
State Preparation Oracle

The state preparation function AA reflects knowledge about the solution space and can be used to eliminate invalid assignments. Here we assume no knowledge of the solution space; hence, we use the uniform superposition state preparation. 
sp_oracle = hadamard_transform

​

​
Wrapping All to the Phase Estimation

We will achieve the desired precision only in the IQAE phase. Here, we compute the worst-case precision for five phase qubits: 
import numpy as np

NUM_PHASE_QUBITS = 5

x = np.linspace(0, 1, 100)
(2**DOMAIN_SIZE) * max(
    np.abs(
        np.sin(np.pi * x) ** 2 - np.sin(np.pi * (x - 1 / (2**NUM_PHASE_QUBITS))) ** 2
    )
)
Output:
1.5681439279637486
We use the grover_operator library function, where we plug in the defined oracles. We wrap it by the qpe function, that accepts the grover operator as a unitary, and additional phase register. 
@qfunc
def main(
    phase_reg: Output[QNum[NUM_PHASE_QUBITS, UNSIGNED, NUM_PHASE_QUBITS]],
) -> None:
    state_reg = OracleVars()
    allocate(state_reg)
    allocate(phase_reg)

    sp_oracle(state_reg)
    qpe(
        unitary=lambda: grover_operator(
            arith_oracle,
            sp_oracle,
            state_reg,
        ),
        phase=phase_reg,
    )
    drop(state_reg)

​

​
Synthesizing the Model to a Quantum Program

qprog_qpe = synthesize(
    main,
    constraints=Constraints(max_width=14),
    preferences=Preferences(optimization_level=1),
)
show(qprog_qpe)
Output:

Quantum program link: https://platform.classiq.io/circuit/32pSFY9VF2TayfnH2DMe7YFeGcl

​

​
Executing the Quantum Program

result = execute(qprog_qpe).result_value()
Upon plotting the resulting histogram, we see two phase values with high probability (however, both correspond to the same amplitude). Note that phase_reg is already coded as fixed QNum in the range [0,1]. 
import matplotlib.pyplot as plt

phases_counts = dict(
    (sampled_state.state["phase_reg"], sampled_state.shots)
    for sampled_state in result.parsed_counts
)
plt.bar(phases_counts.keys(), phases_counts.values(), width=0.1)
plt.xticks(rotation=90)
print("phase with max probability: ", max(phases_counts, key=phases_counts.get))
Output:
phase with max probability:  0.78125
output From the phase, we can extract the number of solutions: 
expected_num_solutions = 6
solutions_ratio_qpe = np.sin(np.pi * max(phases_counts, key=phases_counts.get)) ** 2
print(
    "Number of solutions: ",
    (2**DOMAIN_SIZE) * solutions_ratio_qpe,
)
assert np.isclose(
    (2**DOMAIN_SIZE) * solutions_ratio_qpe, expected_num_solutions, atol=1.5
)
Output:

Number of solutions:  6.439277423870974

​
Amplitude Estimation Using Iterative Quantum Amplitude Estimation

Now we are ready for the iterative method. Instead of QPE, the algorithm applies the unitary (Q)mA(Q)^mA where mm, the number of repetitions, changes between iterations of the algorithm. There is one subtlety that changes the way we work with the Grover operator. The classical algorithm expects an additional indicator qubit that marks the “good” states, i.e.: ∣a⟩∣b⟩∣f(a,b)⟩|a\rangle|b\rangle|f(a,b)\rangle So now, most of our logic goes into the state preparation oracle (AA). It combines the loading of the solution space with setting the indicator qubit. 
@qfunc
def iqae_state_preparation(vars: OracleVars, ind: QBit):
    hadamard_transform(vars)
    arith_equation(vars, ind)

​

​
Wrapping All to the Iterative Quantum Amplitude Estimation Algorithm

We use the built-in IQAE class that the quantum code for the algorithm as well as the classical execution code. The circuit starts with the state A∣0⟩A|0\rangle, then applies iterations of the Grover operator. Note that the algorithm applies a varied number of Grover iterations on each execution. The number of iterations is chosen dynamically based on previous execution results, using statistical inference methods. It expects a state preparation function that creates the following state: ∣Ψ⟩=a∣Ψ1⟩∣1⟩ind+1−a2∣Ψ0⟩∣0⟩ind|\Psi\rangle = a|\Psi_1\rangle|1\rangle_{ind} + \sqrt{1-a^2}|\Psi_0\rangle|0\rangle_{ind} Where the indicator qubit is marking the wanted state. 
from classiq.applications.iqae.iqae import IQAE

iqae = IQAE(
    state_prep_op=iqae_state_preparation,
    problem_vars_size=DOMAIN_SIZE,
    constraints=Constraints(optimization_parameter="width"),
    preferences=Preferences(optimization_level=1),
)

​

​
Synthesizing the Model to a Quantum Program

qprog_iqae = iqae.get_qprog()
show(qprog_iqae)
Output:

Quantum program link: https://platform.classiq.io/circuit/32pSIs3apHDldfDrbzHHHCWW7UY

​

​
Executing the Quantum Program

iqae_result = iqae.run(
    epsilon=1 / (2**DOMAIN_SIZE * 2),  # desired error
    alpha=0.01,  # desired probability for error
)
We set ϵ=1/24⋅0.5=1/32\epsilon = 1/{2^4} \cdot 0.5 = 1/32. alpha is the tail probability of estimating the result with accuracy ϵ\epsilon. 
print(
    f"Number of solutions: {(2**DOMAIN_SIZE) * iqae_result.estimation}, accuracy: "
    f"{(2**DOMAIN_SIZE)*(iqae_result.confidence_interval[1]-iqae_result.confidence_interval[0])}"
)
Output:

Number of solutions: 6.011647242198923, accuracy: 0.09716976875854222
  


assert np.isclose(
    (2**DOMAIN_SIZE) * iqae_result.estimation, 6, atol=1 / (2**DOMAIN_SIZE - 1)
)
We can also see the statistics of the IQAE execution: 
for i, iteration in enumerate(iqae_result.iterations_data):
    print(
        f"iteration_id: {i}, num grover iterations: {iteration.grover_iterations}, counts: {iteration.sample_results.counts}"
    )
Output:
iteration_id: 0, num grover iterations: 0, counts: {'0': 1268, '1': 780}
  iteration_id: 1, num grover iterations: 6, counts: {'0': 898, '1': 1150}

​
References

[1]: Quantum Counting Algorithm, Wikipedia. [2]: Grinko, D., Gacon, J., Zoufal, C. et al. Iterative quantum amplitude estimation. npj Quantum Inf 7, 52 (2021). [3]: Brassard, G., Hoyer, P., Mosca, M., & Tapp, A. (2002). Quantum Amplitude Amplification and Estimation. Contemporary Mathematics, 305, 53-74.