Skip to main content

Documentation Index

Fetch the complete documentation index at: https://prod-mint.classiq.io/llms.txt

Use this file to discover all available pages before exploring further.

View on GitHub

Open this notebook in GitHub to run it yourself

​
Introduction

In this tutorial, we implement a Quantum Monte Carlo Integration (QMCI) circuit to estimate π\pi. The approach samples lattice points on a discretized grid and uses Quantum Amplitude Estimation (QAE) to estimate the fraction of points that lie inside a quarter circle. This fraction directly determines the value of π\pi. Following Ref.~[1], lattice-based uniform sampling is appropriate for the present π\pi-estimation problem because the task reduces to estimating the fraction of marked grid points. In this setting, QAE estimates the corresponding amplitude with a quadratic speedup over classical Monte Carlo methods. This should be distinguished from classical grid-based quadrature, whose cost generally scales exponentially with the dimension. In contrast, the QAE-based QMCI approach can avoid this curse of dimensionality at the level of amplitude estimation. In more general QMCI applications, such as finance, the practical advantage further depends on whether the required probability distribution and oracle can be implemented efficiently.

​
Method

The goal of this tutorial is to estimate π\pi from the area ratio between a square and an inscribed quarter circle. We sample lattice points (x,y)(x,y) on a discretized square and determine whether each point lies inside the quarter circle. If the fraction of points inside the quarter circle is denoted by α\alpha, then α=π4,\alpha = \frac{\pi}{4}, and therefore π=4α.\pi = 4\alpha. Two nn-qubit registers are used to represent the xx- and yy-coordinates, respectively: 0≤x≤2n−1,0≤y≤2n−1.0 \le x \le 2^n-1, \qquad 0 \le y \le 2^n-1. Thus, the computational basis corresponds to 22n2^{2n} lattice points in a square of side length 2n2^n. We define the indicator function f(x,y)={1if x2+y2<22n,0if x2+y2≥22n.f(x,y)= \begin{cases} 1 & \text{if } x^2+y^2 < 2^{2n},\\ 0 & \text{if } x^2+y^2 \ge 2^{2n}. \end{cases} With uniform sampling, p(x,y)=122n,p(x,y)=\frac{1}{2^{2n}}, the target quantity is the expectation value α=∑x=02n−1∑y=02n−1f(x,y)p(x,y).\alpha = \sum_{x=0}^{2^n-1} \sum_{y=0}^{2^n-1} f(x,y)p(x,y). This is the probability that a uniformly sampled lattice point lies inside the quarter circle. To estimate this quantity on a quantum computer, we prepare the uniform superposition over all lattice points using P^=H⊗2n,\hat{P}=H^{\otimes 2n}, so that (P^⊗I)∣0⟩2n∣0⟩=122n∑x=02n−1∑y=02n−1∣x,y⟩2n∣0⟩.(\hat{P}\otimes I)|0\rangle_{2n}|0\rangle = \frac{1}{\sqrt{2^{2n}}} \sum_{x=0}^{2^n-1} \sum_{y=0}^{2^n-1} |x,y\rangle_{2n}|0\rangle . Then an oracle R^\hat{R} marks whether each point is inside the quarter circle by writing the value of f(x,y)f(x,y) into an ancilla qubit: R^∣x,y⟩2n∣0⟩=∣x,y⟩2n∣f(x,y)⟩.\hat{R}|x,y\rangle_{2n}|0\rangle= |x,y\rangle_{2n}|f(x,y)\rangle . Thus, R^(P^⊗I)∣0⟩2n∣0⟩=122n∑x=02n−1∑y=02n−1∣x,y⟩2n∣f(x,y)⟩.\hat{R}(\hat{P}\otimes I)|0\rangle_{2n}|0\rangle = \frac{1}{\sqrt{2^{2n}}} \sum_{x=0}^{2^n-1} \sum_{y=0}^{2^n-1} |x,y\rangle_{2n}|f(x,y)\rangle . Therefore, the probability of measuring the ancilla qubit in the state ∣1⟩|1\rangle is exactly α\alpha. By applying Quantum Amplitude Estimation (QAE), we estimate this amplitude more efficiently than by classical Monte Carlo integration, and finally compute π=4α.\pi = 4\alpha. Compared with classical grid-based quadrature, whose cost grows exponentially with the dimension, this QAE-based approach avoids the curse of dimensionality. Compared with classical Monte Carlo integration, QAE provides a quadratic speedup in the estimation error. In this π\pi-estimation task, lattice-based uniform sampling is appropriate because the problem reduces to estimating the fraction of marked grid points. In more general QMCI applications, such as finance, the advantage additionally depends on whether the required probability distribution and oracle can be implemented efficiently.

​
Dataset

We first generate the grid data. nn is the parameter for number of data. 
import matplotlib.pyplot as plt
import numpy as np


def plot_pi_sampling(GRID_POINTS_PER_AXIS):
    limit = 2**GRID_POINTS_PER_AXIS
    x_range = np.arange(limit)
    y_range = np.arange(limit)

    X, Y = np.meshgrid(x_range, y_range)

    # constraint: f(x, y) = 1 (True) if x^2 + y^2 < limit^2
    inside = X**2 + Y**2 < limit**2

    plt.figure(figsize=(3, 3))

    # data plot
    plt.scatter(X[inside], Y[inside], color="blue", s=10, label="f(x,y)=1 (Inside)")
    plt.scatter(X[~inside], Y[~inside], color="red", s=10, label="f(x,y)=0 (Outside)")

    theta = np.linspace(0, np.pi / 2, 100)
    plt.plot(limit * np.cos(theta), limit * np.sin(theta), color="black", linewidth=2)

    plt.xlim(-0.5, limit)
    plt.ylim(-0.5, limit)
    plt.gca().set_aspect("equal", adjustable="box")
    plt.title(f"Lattice Points (n={GRID_POINTS_PER_AXIS})")
    plt.xlabel("x-axis")
    plt.ylabel("y-axis")
    plt.legend()
    plt.grid(True, linestyle="--", alpha=0.6)
    plt.show()


plot_pi_sampling(GRID_POINTS_PER_AXIS=4)
output

​
Simulation: Classical Simulation

Below, as a point of comparison, we first present an example of a classical simulation based on uniformly spaced sampling points. 
def estimate_pi_grid(GRID_POINTS_PER_AXIS):
    limit = GRID_POINTS_PER_AXIS
    inside_circle = 0
    total_samples = limit * limit

    for y in range(limit):
        for x in range(limit):
            if x**2 + y**2 < limit**2:
                inside_circle += 1

    pi_estimated = (inside_circle / total_samples) * 4
    return pi_estimated


LIST_GRID_POINTS_PER_AXIS = [3, 5, 6]
for mi in LIST_GRID_POINTS_PER_AXIS:
    print(f"(Total Sampling: {mi}²={mi**2}), Estimated Pi: {estimate_pi_grid(mi)}")
Output:
(Total Sampling: 3²=9), Estimated Pi: 4.0
  (Total Sampling: 5²=25), Estimated Pi: 3.52
  (Total Sampling: 6²=36), Estimated Pi: 3.6666666666666665

​
Quantum Simulation

​

​
Option 1: QAE based QMCI

import numpy as np

from classiq import *

N_QUBITS = 3


class position(QStruct):
    x: QNum[N_QUBITS]
    y: QNum[N_QUBITS]


@qfunc
def encode_prob(p: position):
    hadamard_transform(p)


@qperm
def oracle_func(p: Const[position], res: QBit):
    fn = 2 ** (2 * N_QUBITS)
    res ^= p.x**2 + p.y**2 < fn


@qfunc
def oracle_op(state: position):
    phase_oracle(oracle_func, state)


@qfunc
def my_grover_operator(state: position):
    grover_operator(
        oracle_op,
        hadamard_transform,
        state,
    )


QPE_SIZE = 3


@qfunc
def main(phase_reg: Output[QNum[QPE_SIZE, UNSIGNED, QPE_SIZE]]):
    state_reg = position()
    allocate(state_reg.size, state_reg)

    allocate(phase_reg)
    encode_prob(state_reg)
    qpe(
        unitary=lambda: my_grover_operator(state_reg),
        phase=phase_reg,
    )
    drop(state_reg)


qprog = synthesize(main)
show(qprog)
Output:

Quantum program link: https://platform.classiq.io/circuit/3CvTYWFwNzHPWUxN8jrhrtd6MAx
Output:
https://platform.classiq.io/circuit/3CvTYWFwNzHPWUxN8jrhrtd6MAx?login=True&version=18
  


result = execute(qprog).result_value()
phases_counts = dict(
    (sampled_state.state["phase_reg"], sampled_state.shots)
    for sampled_state in result.parsed_counts
)


expected_alpha = np.sin(np.pi * max(phases_counts, key=phases_counts.get))
print("QMCI pi = ", expected_alpha * 4)
print("exact pi = ", np.pi)
Output:

QMCI pi =  3.695518130045147
  exact pi =  3.141592653589793
This result indicates that QMCI successfully approaches the value of π\pi, although the current estimate still has a significant error compared with the exact value. In the next cell, we investigate how accurately π\pi can be estimated by applying IQAE and examining its achievable precision.

​

​
Option2: Iterative QAE based QMCI

from classiq.applications.iqae.iqae import IQAE


def iqae_circuit(N_QUBITS):
    class position(QStruct):
        x: QNum[N_QUBITS]
        y: QNum[N_QUBITS]

    @qfunc
    def encode_prob(p: position):
        hadamard_transform(p)

    @qperm
    def oracle_func(p: Const[position], res: QBit):
        fn = 2 ** (2 * N_QUBITS)
        res ^= p.x**2 + p.y**2 < fn

    @qfunc
    def iqae_state_preparation(p: position, ind: QBit):
        encode_prob(p)
        oracle_func(p, ind)

    iqae = IQAE(
        state_prep_op=iqae_state_preparation,
        problem_vars_size=2 * N_QUBITS,
        constraints=Constraints(max_width=30),
        preferences=Preferences(machine_precision=N_QUBITS),
    )
    return iqae


N_QUBITS = 3
qprog = iqae_circuit(N_QUBITS).get_qprog()
show(qprog)
Output:

Quantum program link: https://platform.classiq.io/circuit/3CvUMKYNGscQI6wWJm3Bjqf3v9q
Output:
https://platform.classiq.io/circuit/3CvUMKYNGscQI6wWJm3Bjqf3v9q?login=True&version=18

​

​
IAQE Parameters

For the numerical test, we use the standard Classiq IAQE routine. The two main parameters are the target accuracy ϵ\epsilon and the failure probability α\alpha. We set ϵ=0.03,α=0.01.\epsilon = 0.03, \qquad \alpha = 0.01. Here, ϵ\epsilon specifies the target additive accuracy of the amplitude estimation, while α\alpha specifies the allowed failure probability. Thus, α=0.01\alpha=0.01 corresponds to a 9999% confidence level. This setting is suitable for an initial validation because it keeps the estimation cost moderate while still requiring a reliable confidence level. Once the implementation is verified, ϵ\epsilon can be reduced for a higher-precision estimation. 
iqae_res = iqae_circuit(N_QUBITS).run(epsilon=0.03, alpha=0.01)


iqae_res.confidence_interval
Output:
[0.8735326790922084, 0.8783898719147034]
  


iqae_res.estimation * 4
Output:
3.5038451020138237
The factor of 44 is required because IQAE estimates π/4\pi/4, corresponding to the area ratio of the quarter circle inside the square.

​

​
Experiment

In this experiment, we examine whether the estimated value approaches the exact solution as the number of qubits increases. Because a larger number of qubits allows exponentially more sampling points to be represented, the discretization becomes finer, and the estimate is expected to become more accurate. By tracking this trend, we can directly evaluate how increasing the qubit number improves the precision of the π\pi estimation. 
list_expected_value = []
list_conf_err = []
for N_QUBITS in range(1, 5):
    iqae_res = iqae_circuit(N_QUBITS).run(epsilon=0.03, alpha=0.01)
    print("expected pi =", iqae_res.estimation * 4)
    list_expected_value.append(iqae_res.estimation * 4)
    list_conf_err.append(iqae_res.confidence_interval)
Output:
expected pi = 3.998328879248948
  expected pi = 3.744334365027651
  expected pi = 3.499632929603465
  expected pi = 3.34139467918454
When you run the cell above, you should obtain results similar to those shown below: 
# list_expected_value = [3.998328879248948, 3.7572795311895337, 3.5043386707175443, 3.3476167564481774]
# list_conf_err = [[0.999164439624474, 1.0], [0.9328714567271077, 0.9457683088676591], [0.8720696958863934, 0.8800996394723789], [0.8298663183601778, 0.8439420598639109]]


import matplotlib.pyplot as plt
import numpy as np

list_N_QUBITS = range(1, 5)
total_sampling = [2 ** (2 * N_QUBITS) for N_QUBITS in list_N_QUBITS]

pi_ests = list_expected_value

yerr_lower = [pi_ests[i] - (list_conf_err[i][0] * 4) for i in range(len(pi_ests))]
yerr_upper = [(list_conf_err[i][1] * 4) - pi_ests[i] for i in range(len(pi_ests))]

plt.figure(figsize=(4, 4))
plt.tick_params(direction="in", which="both", top=True, right=True)

plt.errorbar(
    total_sampling,
    pi_ests,
    yerr=[yerr_lower, yerr_upper],
    fmt="o",
    color="red",
    capsize=5,
    markersize=6,
    label="IQAE based QMCI",
)

LIST_GRID_POINTS_PER_AXIS = np.array(
    [2, 3, 4, 5, 6, 7, 8, 9, 10, 20, 30, 40, 50, 60, 70, 80]
)
classical_pi_vals = [estimate_pi_grid(mi) for mi in LIST_GRID_POINTS_PER_AXIS]
plt.plot(
    LIST_GRID_POINTS_PER_AXIS**2,
    classical_pi_vals,
    "--",
    color="black",
    label="Classical Calculation",
)

plt.axhline(y=np.pi, color="green", linestyle="-.", alpha=0.6, label="π (True Value)")

plt.xlabel("Total sampling(= Grid Size)")
plt.ylabel("Estimated π")
plt.ylim(2.8, 4.2)
plt.xscale("log")
plt.legend(loc="upper right")
plt.tight_layout()
plt.show()
output We observe convergence toward the exact value in both the classical and IQAE-based quantum estimates of π\pi as the number of sampling points increases.

​
Discussion

How much computational advantage can be achieved? According to Ref. [2], the query complexity of QMCI is given by Nquery≈C⋅ln⁡!(2α)ϵ.N_{\mathrm{query}} \approx \frac{C \cdot \ln!\left(\frac{2}{\alpha}\right)}{\epsilon}. By contrast, in the classical grid-based calculation considered here, all lattice points are evaluated explicitly. Hence, the computational cost is tclassical=22n=O(1/ϵ2).t_{\mathrm{classical}} = 2^{2n} = \mathcal{O}(1/\epsilon^2). The table below shows a theoretical comparison of the required number of evaluations and queries.
ItemClassical Grid-Based IntegrationQuantum Monte Carlo Integration (IQAE)
n=2n=224=162^4 = 1612,71612,716
n=4n=428=2562^8 = 25612,71612,716
n=10n=10220=1,048,5762^{20} = 1,048,57612,71612,716
This comparison indicates that the quantum method does not provide an advantage at small scales, but becomes more efficient than the classical approach once the problem size is sufficiently large.

​
Reference

[1] Quantum circuit to estimate pi using quantum amplitude estimation [2] Iterative Quantum Amplitude Estimation, Granko et al., 2019