quantum_monte_carlo(tape, wires, target_wire, estimation_wires)[source]

Applies the transform quantum Monte Carlo estimation algorithm.

The input tape` should be the quantum circuit corresponding to the \(\mathcal{F}\) unitary in the paper above. This unitary encodes the probability distribution and random variable onto wires so that measurement of the target_wire provides the expectation value to be estimated. The quantum Monte Carlo algorithm then estimates the expectation value using quantum phase estimation (check out QuantumPhaseEstimation for more details), using the estimation_wires.


A complementary approach for quantum Monte Carlo is available with the QuantumMonteCarlo template.

The quantum_monte_carlo transform is intended for use when you already have the circuit for performing \(\mathcal{F}\) set up, and is compatible with resource estimation and potential hardware implementation. The QuantumMonteCarlo template is only compatible with simulators, but may perform faster and is suited to quick prototyping.

  • tape (QNode or QuantumTape or Callable) – the quantum circuit that applies quantum operations according to the \(\mathcal{F}\) unitary used as part of quantum Monte Carlo estimation

  • wires (Union[Wires or Sequence[int]]) – the wires acted upon by the fn circuit

  • target_wire (Union[Wires, int]) – The wire in which the expectation value is encoded. Must be contained within wires.

  • estimation_wires (Union[Wires, Sequence[int], or int]) – the wires used for phase estimation


The transformed circuit as described in qml.transform. Executing this circuit will perform the quantum Monte Carlo estimation.

Return type

qnode (QNode) or quantum function (Callable) or tuple[List[QuantumTape], function]


ValueError – if wires and estimation_wires share a common wire

Consider an input quantum circuit fn that performs the unitary

\[\mathcal{F} = \mathcal{R} \mathcal{A}.\]

Here, the unitary \(\mathcal{A}\) prepares a probability distribution \(p(i)\) of dimension \(M = 2^{m}\) over \(m \geq 1\) qubits:

\[\mathcal{A}|0\rangle^{\otimes m} = \sum_{i \in X} p(i) |i\rangle,\]

where \(X = \{0, 1, \ldots, M - 1\}\) and \(|i\rangle\) is the basis state corresponding to \(i\). The \(\mathcal{R}\) unitary imprints the result of a function \(f: X \rightarrow [0, 1]\) onto an ancilla qubit:

\[\mathcal{R}|i\rangle |0\rangle = |i\rangle \left(\sqrt{1 - f(i)} |0\rangle + \sqrt{f(i)}|1\rangle\right).\]

Following this paper, the probability of measuring the state \(|1\rangle\) in the final qubit is

\[\mu = \sum_{i \in X} p(i) f(i).\]

However, it is possible to measure \(\mu\) more efficiently using quantum Monte Carlo estimation. This function transforms an input quantum circuit fn that performs the unitary \(\mathcal{F}\) to a larger circuit for measuring \(\mu\) using the quantum Monte Carlo algorithm.


The algorithm proceeds as follows:

  1. The probability distribution \(p(i)\) is encoded using a unitary \(\mathcal{A}\) applied to the first \(m\) qubits specified by wires.

  2. The function \(f(i)\) is encoded onto the target_wire using a unitary \(\mathcal{R}\).

  3. The unitary \(\mathcal{Q}\) is defined with eigenvalues \(e^{\pm 2 \pi i \theta}\) such that the phase \(\theta\) encodes the expectation value through the equation \(\mu = (1 + \cos (\pi \theta)) / 2\). The circuit in steps 1 and 2 prepares an equal superposition over the two states corresponding to the eigenvalues \(e^{\pm 2 \pi i \theta}\).

  4. The circuit returned by this function is applied so that \(\pm\theta\) can be estimated by finding the probabilities of the \(n\) estimation wires. This in turn allows for the estimation of \(\mu\).

Visit Rebentrost et al. (2018) for further details. In this algorithm, the number of applications \(N\) of the \(\mathcal{Q}\) unitary scales as \(2^{n}\). However, due to the use of quantum phase estimation, the error \(\epsilon\) scales as \(\mathcal{O}(2^{-n})\). Hence,

\[N = \mathcal{O}\left(\frac{1}{\epsilon}\right).\]

This scaling can be compared to standard Monte Carlo estimation, where \(N\) samples are generated from the probability distribution and the average over \(f\) is taken. In that case,

\[N = \mathcal{O}\left(\frac{1}{\epsilon^{2}}\right).\]

Hence, the quantum Monte Carlo algorithm has a quadratically improved time complexity with \(N\).


Consider a standard normal distribution \(p(x)\) and a function \(f(x) = \sin ^{2} (x)\). The expectation value of \(f(x)\) is \(\int_{-\infty}^{\infty}f(x)p(x)dx \approx 0.432332\). This number can be approximated by discretizing the problem and using the quantum Monte Carlo algorithm.

First, the problem is discretized:

from scipy.stats import norm

m = 5
M = 2 ** m

xmax = np.pi  # bound to region [-pi, pi]
xs = np.linspace(-xmax, xmax, M)

probs = np.array([norm().pdf(x) for x in xs])
probs /= np.sum(probs)

func = lambda i: np.sin(xs[i]) ** 2
r_rotations = np.array([2 * np.arcsin(np.sqrt(func(i))) for i in range(M)])

The quantum_monte_carlo transform can then be used:

from pennylane.templates.state_preparations.mottonen import (
    _apply_uniform_rotation_dagger as r_unitary,

n = 6
N = 2 ** n

a_wires = range(m)
wires = range(m + 1)
target_wire = m
estimation_wires = range(m + 1, n + m + 1)

dev = qml.device("default.qubit", wires=(n + m + 1))

def fn():
    qml.templates.MottonenStatePreparation(np.sqrt(probs), wires=a_wires)
    r_unitary(qml.RY, r_rotations, control_wires=a_wires[::-1], target_wire=target_wire)

def qmc():
    qml.quantum_monte_carlo(fn, wires, target_wire, estimation_wires)()
    return qml.probs(estimation_wires)

phase_estimated = np.argmax(qmc()[:int(N / 2)]) / N

The estimated value can be retrieved using the formula \(\mu = (1-\cos(\pi \theta))/2\)

>>> (1 - np.cos(np.pi * phase_estimated)) / 2

It is also possible to explore the resources required to perform the quantum Monte Carlo algorithm

>>> qml.specs(qmc, expansion_strategy="device")()
{'resources': Resources(
    gate_types=defaultdict(<class 'int'>, {'RY': 7747, 'CNOT': 7874, 'Hadamard': 258, 'CZ': 126, 'Adjoint(CNOT)': 7812, 'Adjoint(RY)': 7686, 'PauliX': 252, 'MultiControlledX': 126, 'Adjoint(QFT)': 1}),
    gate_sizes=defaultdict(<class 'int'>, {1: 15943, 2: 15812, 7: 126, 6: 1}), depth=30610, shots=Shots(total_shots=None, shot_vector=()),
 'num_observables': 1,
 'num_diagonalizing_gates': 0,
 'num_trainable_params': 15433,
 'num_device_wires': 12,
 'device_name': 'default.qubit',
 'expansion_strategy': 'gradient',
 'gradient_options': {},
 'interface': 'auto',
 'diff_method': 'best',
 'gradient_fn': 'backprop'}