monte_carlo#

Functionality#

The Quantum Monte Carlo Sampling algorithm aims to estimate the expected value of some function of a random variable \(\mathbb{E}_{x \sim p}[f(x)]\) given quantum oracle access to the probability distribution of the variable \(p(x)\) and the function of interest \(f(\cdot)\). For a distribution which is potentially not analytically known or otherwise difficult to calculate directly, it might not be feasible to get an exact result for the expectation value of the function applied to the distribution. Both the classical and quantum versions of the Monte Carlo Sampling algorithm utilize methods to approximate this expectation value by utilizing calls to the distribution.

`monte_carlo.py`#

Callee: QuantumMonteCarlo

With the quantum algorithm, the probability distribution and function are instead coded into calls to a quantum oracle. We can then use the ideas of Quantum Amplitude Estimation to return an estimate of the exectation value which has variance decreasing quadratically faster than ordinary classical sampling.

The Monte Carlo algorithm is benchmarked in two different ways:

method=1: General Monte Carlo benchmark with user-specified distribution and objective function. Example distributions for benchmarking may include Gaussians, but any discrete probability distribution can be used, provided it can be encoded into the quantum state.
method=2: We run only a single circuit with the specific distribution: \(p(x) = \frac{1}{2^n}\) for all \(x\) in \([0, 2^n - 1]\), i.e., a uniform distribution. The function \(f(x)\) is set to be \(f(x) = x \text{ mod } 2\), aiming to indicate whether \(x\) is even or odd.

The following table shows the qubit layout for the Monte Carlo circuit:

Index	\(0 \sim (n-1)\)	\(n\)	\((n+1) \sim (n+m-1)\)
Register	State qubits	Object qubit	Counting qubits

To encode information about our distribution and function we define the following operators which will be called as oracles; \(\mathcal{R}\), \(\mathcal{F}\) defined on \(n+1\) qubits such that (the following mapping only includes the state and objective qubits):

\[ (\mathbb{I}_{2\times 2}\otimes\mathcal{R}) \ket{0}_1 \ket{0}_n = \sum_{x} \sqrt{p(x)} \ket{0}_1 \ket{x}_n \]

\[ \mathcal{F} \ket{0}_1 \ket{x}_n = \left( \sqrt{1 - f(x)} \ket{0}_1 + \sqrt{f(x)} \ket{1}_1 \right) \ket{x}_n \]

In this context, prepares the zero state into the all real-valued superposition according to the desired measurement probabilities, while \(\mathcal{F}\) encodes the value of the function \(f(x)\) of the state of the \(n\) qubit register \(\ket{x}\) into the objective qubit at the end. The product of these two is:

\[ \mathcal{F}(\mathbb{I}_{2\times 2}\otimes\mathcal{R}) \ket{0}_1 \ket{0}_n = \sum_{x} \sqrt{p(x)} \left( \sqrt{1 - f(x)} \ket{0}_1 + \sqrt{f(x)} \ket{1}_1 \right) \ket{x}_n \]

If we remember that the expection value we are looking for is defined as

\[ \mathbb{E}_{x \sim p}[f(x)] \equiv \sum_{x} p(x) f(x), \]

we can see that if there is an efficient way of estimating the amplitude on the objective qubit in the state, there is an efficient way to estimate the expected value of the function on the random variable \(a = \mathbb{E}_{x \sim p}[f(x)]\). Then, defining the operator

\[ \mathcal{F}(\mathbb{I}_{2\times 2}\otimes\mathcal{R}) \ket{0}_{n+1} = \sqrt{1 - a} \ket{0}_1 \ket{\psi_0}_{n} + \sqrt{a} \ket{1}_1 \ket{\psi_1}_{n}, \]

Amplitude Estimation will allow us to find the value \(a\) by introducing \(m\) counting qubits and measuring them appropriately. Let \(\ket{\psi_1}\) is the good state, while \(\ket{\psi_0}\) is the bad state, where the amplitude on the good state is what we want to estimate. Amplitude Estimation defines the Grover operator as:

\[ \mathcal{Q} = -\mathcal{A} \mathcal{S}_0 \mathcal{A}^\dagger \mathcal{S}_\chi, \]

where

\(\mathcal{A} = \mathcal{F}(\mathbb{I}_{2\times 2}\otimes\mathcal{R}) \in \mathrm{U}(2^{n+1})\);
\(\mathcal{S}_0\) is the reflection about the zero state, i.e., \(\mathcal{S}_0 = \mathbb{I}_{2^{n+1} \times 2^{n+1}} - 2 \ket{0}_{n+1} \langle 0|_{n+1}\);
\(\mathcal{S}_\chi\) is the reflection about the good state, i.e., \(\mathcal{S}_\chi = \mathbb{I}_{2^{n+1} \times 2^{n+1}} - 2 P_{\text{good}}\), where \(P_{\text{good}}\) is the projector onto the good state, i.e., \(P_{\text{good}} = \ket{1}_1 \bra{1}_1 \otimes \mathbb{I}_{2^n \times 2^n}\).

In the two-dimensional subspace spanned by \(\{ \ket{\psi_0}, \ket{\psi_1} \}\), the operator \(\mathcal{Q}\) is denoted as

\[\begin{split} \mathcal{Q} = \begin{bmatrix} \cos(2\theta) & -\sin(2\theta) \\ \sin(2\theta) & \cos(2\theta) \end{bmatrix}, \end{split}\]

where \(\sin^2(\theta) = a\).

Initialize the counting qubits in the \(\ket{0}^{\otimes m}\) state uising Hadamard gates. Then, apply the controlled-\(\mathcal{Q}^{2^j}\) operations for \(j = 0, 1, \ldots, m-1\). The joint state of the system is:

\[ \ket{\Psi} = \frac{1}{2^{m/2}} \sum_{k=0}^{2^m - 1} \ket{k}_m \otimes \mathcal{Q}^k \mathcal{A} \ket{0}_{n+1} = \frac{1}{2^{m/2}} \sum_{k=0}^{2^m - 1} \ket{k}_m \otimes \left( \cos((2k + 1)\theta) \ket{\psi_0} + \sin((2k + 1)\theta) \ket{\psi_1} \right) \]

Next, perform an inverse Quantum Fourier Transform, and the final state approximately turns out to be:

\[ \ket{\Phi} \approx \frac{1}{\sqrt{2}} \left( \ket{y_+}_m \otimes \ket{\psi_+}_{n+1} - \ket{y_-}_m \otimes \ket{\psi_-}_{n+1} \right). \]

Like most cases using quantum phase estimation, measuring the counting qubits will approximately yield either \(y_{+}\) or \(y_{-}\), whereas this does not guarantee exact two values. Here, we let \(\delta\) as the error term, and

\[ y_{\pm} = \left\lfloor 2^m \left( \frac{\theta}{\pi} \pm \delta \right) \right\rfloor \mod 2^m. \]

Approximately,

\[ y_{+} \approx 2^m \frac{\theta}{\pi}, \quad y_{-} \approx 2^m \left( 1 - \frac{\theta}{\pi} \right). \]

Then, we can derive an estimation \(\hat{a}\) of \(a\) as follows:

\[ \hat{a} = \sin^2 \left( \frac{\pi y_{+}}{2^m} \right) \quad \text{or} \quad \hat{a} = \sin^2 \left( \frac{\pi (2^m - y_{-})}{2^m} \right). \]

Note:

Considering the validity of \(\sqrt{1-f(x)}\) and \(\sqrt{f(x)}\), the function \(f(x)\) must satisfy \(0 \leq f(x) \leq 1\) for all \(x\) in the domain.
According to the literature [23], the error bound of the estimation \(\hat{a}\) can be proven as

\[ |\hat{a} - a| \leq \frac{2\pi \sqrt{a(1-a)}}{M} + \frac{\pi^2}{M^2} \leq \frac{\pi}{M} + \frac{\pi^2}{M^2}=O(M^{-1}), \]

with probability at least \(8/\pi^2\), where \(M\) is the number of samples, i.e., \(M=2^m\).

Doc references: [24]

Technical Debt for Testing#

Regarding the provided test cases in test_monte_carlo.py, it could be found that several tests do not strictly pass for method=1, although five tests can pass. This does not always imply a bug in the implementation of the QuantumMonteCarlo class itself without more convincing evidence. The test cases and the test oracles may not be perfectly aligned with each other, owing to the substantial difficulty in capturing the mathematically approximate operation involved in the code. For example, the root cause may lie in the fact that the function _f_on_objective employs polynomial fitting to approximate amplitudes, and the resulting error is influenced by the parameters degree and epsilon. Thus, further investigation and refinement of the test cases and test oracles may be necessary to ensure their validity and reliability.

The following code appends several test inputs that cannot strictly pass in the current test suite for method=1.

@pytest.mark.parametrize(
    "num_state_qubits, num_counting_qubits, method, target_dist, f, expected_a", [
        # Delta distribution with indicator function (expectation = 1)
        (3, 5, 1, dist_delta(3, 5), f_first_quarter, 1.0),
        # Uniform distribution + parity function (expectation = 1/2)
        (3, 4, 1, dist_uniform(3), f_parity, 0.5),
        # Uniform distribution + first quarter indicator (expectation = 1/4)
        (4, 5, 1, dist_uniform(4), f_first_quarter, 0.25),
    ],
    ids=["tf0", "tf1", "tf2"],
)
def test_monte_carlo(num_state_qubits, num_counting_qubits, method, target_dist, f, expected_a):
    ...

API#

qolumbina.programs.monte_carlo.monte_carlo#

Monte Carlo Sampling Benchmark Program via Amplitude Estimation- Qiskit

class QuantumMonteCarlo(num_state_qubits, num_counting_qubits, method, target_dist=None, f=None, epsilon=0.05, degree=2, if_barrier=False, name=None)[source]#

Bases: QuantumCircuit

Monte Carlo Sampling via Amplitude Estimation Circuit

Parameters:

num_state_qubits (int) – Number of state qubits, denoted as \(n\).
num_counting_qubits (int) – Number of counting qubits, denoted as \(m\).
method (Literal[1, 2]) –
Benchmark method (1 or 2)

method=1:
General Monte Carlo benchmark with user-specified distribution and objective function.

method=2:
Fixed, analytically tractable benchmark instance with a predefined distribution and objective function.
target_dist (Callable | None) –
Function representing the target probability distribution \(p(x)\).
1. Required if method=1 and must be explicitly provided by the caller. For convenience, a Gaussian distribution function generator is provided in the module [monte_carlo._common.mc_utils](./_common/mc_utils.py), which can directly generate Gaussian distributions.
2. Must be None if method=2, in which case a uniform distribution over all computational basis states is used internally.
f (Callable | None) –
Function representing the objective function \(f(x)\).
1. Required if method=1 and must be explicitly provided by the caller. For convenience, a polynomial approximation function generator is provided in the module [monte_carlo._common.mc_utils](./_common/mc_utils.py).
2. Must be None if method=2, in which case a fixed binary-valued function \(f(x) = x \bmod 2\) is used internally.
epsilon (float | None) –
Approximation error tolerance for the polynomial approximation of \(f(x)\).
- Used only if method=1.
- Must be None if method=2.
- Default value is 0.05 when method=1.
degree (int | None) –
Degree of the polynomial used to approximate the objective function \(f(x)\).
- Used only if method=1.
- Must be None if method=2.
- Default value is 2 when method=1.
if_barrier (bool) – Whether to insert barriers.
name (str | None) – Optional name of the quantum circuit.

Raises:

ValueError – If method is not 1 or 2.
ValueError – If required arguments for the specified method are missing or if arguments that should be None for the specified method are not None.

AE_Subroutine(num_state_qubits, num_counting_qubits, A_circuit, if_barrier, method)[source]#

Parameters:

num_state_qubits (int)
num_counting_qubits (int)
A_circuit (QuantumCircuit)
if_barrier (bool)
method (int)

Return type:

QuantumCircuit

Ctrl_Q(num_state_qubits, A_circ)[source]#

Construct the Grover-like operator Q for amplitude estimation

Parameters:

num_state_qubits (int)
A_circ (QuantumCircuit)

Return type:

tuple[QuantumCircuit, QuantumCircuit]

f_on_objective(qc, qr, f, epsilon, degree)[source]#

Assume last qubit is the objective. Function f is evaluated on first n-1 qubits

Parameters:

qc (QuantumCircuit)
qr (QuantumRegister)
f (Callable)
epsilon (float)
degree (int)

Return type:

None

square_on_objective(qc)[source]#

Assume last qubit is the objective. Shifted square wave function: if \(x\) is even, \(f(x) = 0\); if \(x\) is odd, \(f(x) = 1\)

Parameters:: qc (QuantumCircuit)
Return type:: None

state_prep(qc, qr, target_dist, num_state_qubits)[source]#

Use controlled Ry gates to construct the superposition \(\sum_i \sqrt{p_i} \ket{i}\) where \(p_i\) is the probability of the \(i^{th}\) state defined by the target distribution.

Parameters:

qc (QuantumCircuit)
qr (QuantumRegister)
target_dist (Callable)
num_state_qubits (int)

Return type:

None

uniform_prep(qc, num_state_qubits)[source]#

Generates a uniform distribution over all states

Parameters:

qc (QuantumCircuit)
num_state_qubits (int)

Return type:

None

monte_carlo

Contents

monte_carlo#

Functionality#

monte_carlo.py#

Technical Debt for Testing#

API#

qolumbina.programs.monte_carlo.monte_carlo#

`monte_carlo.py`#