# Source code for pennylane.templates.subroutines.qmc

# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at

# Unless required by applicable law or agreed to in writing, software
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
"""
Contains the QuantumMonteCarlo template and utility functions.
"""
# pylint: disable=too-many-arguments
import numpy as np
import pennylane as qml
from pennylane.operation import AnyWires, Operation
from pennylane.ops import QubitUnitary

def probs_to_unitary(probs):
r"""Calculates the unitary matrix corresponding to an input probability distribution.

For a given distribution :math:p(i), this function returns the unitary :math:\mathcal{A}
that transforms the :math:|0\rangle state as

.. math::

\mathcal{A} |0\rangle = \sum_{i} \sqrt{p(i)} |i\rangle,

so that measuring the resulting state in the computational basis will give the state
:math:|i\rangle with probability :math:p(i). Note that the returned unitary matrix is
real and hence an orthogonal matrix.

Args:
probs (array): input probability distribution as a flat array

Returns:
array: unitary

Raises:
ValueError: if the input array is not flat or does not correspond to a probability
distribution

**Example:**

>>> p = np.ones(4) / 4
>>> probs_to_unitary(p)
array([[ 0.5       ,  0.5       ,  0.5       ,  0.5       ],
[ 0.5       , -0.83333333,  0.16666667,  0.16666667],
[ 0.5       ,  0.16666667, -0.83333333,  0.16666667],
[ 0.5       ,  0.16666667,  0.16666667, -0.83333333]])
"""

if not qml.math.is_abstract(
sum(probs)
):  # skip check and error if jitting to avoid JAX tracer errors
if not qml.math.allclose(sum(probs), 1) or min(probs) < 0:
raise ValueError(
"A valid probability distribution of non-negative numbers that sum to one "
"must be input"
)

# Using the approach discussed here:
# https://quantumcomputing.stackexchange.com/questions/10239/how-can-i-fill-a-unitary-knowing-only-its-first-column
psi = qml.math.sqrt(probs)
overlap = psi[0]
denominator = qml.math.sqrt(2 + 2 * overlap)
psi = qml.math.set_index(psi, 0, psi[0] + 1)  # psi[0] += 1, but JAX-JIT compatible
psi /= denominator

dim = len(probs)
return 2 * qml.math.outer(psi, psi) - np.eye(dim)

def func_to_unitary(func, M):
r"""Calculates the unitary that encodes a function onto an ancilla qubit register.

Consider a function defined on the set of integers :math:X = \{0, 1, \ldots, M - 1\} whose
output is bounded in the interval :math:[0, 1], i.e., :math:f: X \rightarrow [0, 1].

The func_to_unitary function returns a unitary :math:\mathcal{R} that performs the
transformation:

.. math::

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

In other words, for a given input state :math:|i\rangle \otimes |0\rangle, this unitary
encodes the amplitude :math:\sqrt{f(i)} onto the :math:|1\rangle state of the ancilla qubit.
Hence, measuring the ancilla qubit will result in the :math:|1\rangle state with probability
:math:f(i).

Args:
func (callable): a function defined on the set of integers
:math:X = \{0, 1, \ldots, M - 1\} with output value inside :math:[0, 1]
M (int): the number of integers that the function is defined on

Returns:
array: the :math:\mathcal{R} unitary

Raises:
ValueError: if func is not bounded with :math:[0, 1] for all :math:X

**Example:**

>>> func = lambda i: np.sin(i) ** 2
>>> M = 16
>>> func_to_unitary(func, M)
array([[ 1.        ,  0.        ,  0.        , ...,  0.        ,
0.        ,  0.        ],
[ 0.        , -1.        ,  0.        , ...,  0.        ,
0.        ,  0.        ],
[ 0.        ,  0.        ,  0.54030231, ...,  0.        ,
0.        ,  0.        ],
...,
[ 0.        ,  0.        ,  0.        , ..., -0.13673722,
0.        ,  0.        ],
[ 0.        ,  0.        ,  0.        , ...,  0.        ,
0.75968791,  0.65028784],
[ 0.        ,  0.        ,  0.        , ...,  0.        ,
0.65028784, -0.75968791]])
"""
unitary = np.zeros((2 * M, 2 * M))

fs = [func(i) for i in range(M)]
if not qml.math.is_abstract(
fs[0]
):  # skip check and error if jitting to avoid JAX tracer errors
if min(fs) < 0 or max(fs) > 1:
raise ValueError(
"func must be bounded within the interval [0, 1] for the range of input values"
)

for i, f in enumerate(fs):
# array = set_index(array, idx, val) is a JAX-JIT compatible version of array[idx] = val
unitary = qml.math.set_index(unitary, (2 * i, 2 * i), qml.math.sqrt(1 - f))
unitary = qml.math.set_index(unitary, (2 * i + 1, 2 * i), qml.math.sqrt(f))
unitary = qml.math.set_index(unitary, (2 * i, 2 * i + 1), qml.math.sqrt(f))
unitary = qml.math.set_index(unitary, (2 * i + 1, 2 * i + 1), -qml.math.sqrt(1 - f))

return unitary

def _make_V(dim):
r"""Calculates the :math:\mathcal{V} unitary which performs a reflection along the
:math:|1\rangle state of the end ancilla qubit.

Args:
dim (int): dimension of :math:\mathcal{V}

Returns:
array: the :math:\mathcal{V} unitary
"""
assert dim % 2 == 0, "dimension for _make_V() must be even"

one = np.array([[0, 0], [0, 1]])
dim_without_qubit = int(dim / 2)

return 2 * np.kron(np.eye(dim_without_qubit), one) - np.eye(dim)

def _make_Z(dim):
r"""Calculates the :math:\mathcal{Z} unitary which performs a reflection along the all
:math:|0\rangle state.

Args:
dim (int): dimension of :math:\mathcal{Z}

Returns:
array: the :math:\mathcal{Z} unitary
"""
Z = -np.eye(dim)
Z[0, 0] = 1
return Z

def make_Q(A, R):
r"""Calculates the :math:\mathcal{Q} matrix that encodes the expectation value according to
the probability unitary :math:\mathcal{A} and the function-encoding unitary
:math:\mathcal{R}.

Following this <https://journals.aps.org/pra/abstract/10.1103/PhysRevA.98.022321>__ paper,
the expectation value is encoded as the phase of an eigenvalue of :math:\mathcal{Q}. This
phase can be estimated using quantum phase estimation using the
:func:~.QuantumPhaseEstimation template. See :func:~.QuantumMonteCarlo for more details,
which loads make_Q() internally and applies phase estimation.

Args:
A (array): The unitary matrix of :math:\mathcal{A} which encodes the probability
distribution
R (array): The unitary matrix of :math:\mathcal{R} which encodes the function

Returns:
array: the :math:\mathcal{Q} unitary
"""
A_big = qml.math.kron(A, np.eye(2))
F = R @ A_big
F_dagger = F.conj().T

dim = len(R)
V = _make_V(dim)
Z = _make_Z(dim)
UV = F @ Z @ F_dagger @ V

return UV @ UV

[docs]class QuantumMonteCarlo(Operation): r"""Performs the quantum Monte Carlo estimation <https://arxiv.org/abs/1805.00109>__ algorithm. Given a probability distribution :math:p(i) of dimension :math:M = 2^{m} for some :math:m \geq 1 and a function :math:f: X \rightarrow [0, 1] defined on the set of integers :math:X = \{0, 1, \ldots, M - 1\}, this function implements the algorithm that allows the following expectation value to be estimated: .. math:: \mu = \sum_{i \in X} p(i) f(i). .. figure:: ../../_static/templates/subroutines/qmc.svg :align: center :width: 60% :target: javascript:void(0); Args: probs (array): input probability distribution as a flat array func (callable): input function :math:f defined on the set of integers :math:X = \{0, 1, \ldots, M - 1\} such that :math:f(i)\in [0, 1] for :math:i \in X target_wires (Union[Wires, Sequence[int], or int]): the target wires estimation_wires (Union[Wires, Sequence[int], or int]): the estimation wires Raises: ValueError: if probs is not flat or has a length that is not compatible with target_wires .. note:: This template is only compatible with simulators because the algorithm is performed using unitary matrices. Additionally, this operation is not differentiable. To implement the quantum Monte Carlo algorithm on hardware requires breaking down the unitary matrices into hardware-compatible gates, check out the :func:~.quantum_monte_carlo transformation for more details. .. details:: :title: Usage Details The algorithm proceeds as follows: #. The probability distribution :math:p(i) is encoded using a unitary :math:\mathcal{A} applied to the first :math:m qubits specified by target_wires. #. The function :math:f(i) is encoded onto the last qubit of target_wires using a unitary :math:\mathcal{R}. #. The unitary :math:\mathcal{Q} is defined with eigenvalues :math:e^{\pm 2 \pi i \theta} such that the phase :math:\theta encodes the expectation value through the equation :math:\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 :math:e^{\pm 2 \pi i \theta}. #. The :func:~.QuantumPhaseEstimation circuit is applied so that :math:\pm\theta can be estimated by finding the probabilities of the :math:n estimation wires. This in turn allows for the estimation of :math:\mu. Visit Rebentrost et al. (2018) <https://arxiv.org/abs/1805.00109>__ for further details. In this algorithm, the number of applications :math:N of the :math:\mathcal{Q} unitary scales as :math:2^{n}. However, due to the use of quantum phase estimation, the error :math:\epsilon scales as :math:\mathcal{O}(2^{-n}). Hence, .. math:: N = \mathcal{O}\left(\frac{1}{\epsilon}\right). This scaling can be compared to standard Monte Carlo estimation, where :math:N samples are generated from the probability distribution and the average over :math:f is taken. In that case, .. math:: N = \mathcal{O}\left(\frac{1}{\epsilon^{2}}\right). Hence, the quantum Monte Carlo algorithm has a quadratically improved time complexity with :math:N. An example use case is given below. Consider a standard normal distribution :math:p(x) and a function :math:f(x) = \sin ^{2} (x). The expectation value of :math:f(x) is :math:\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: .. code-block:: python 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 The QuantumMonteCarlo template can then be used: .. code-block:: n = 10 N = 2 ** n target_wires = range(m + 1) estimation_wires = range(m + 1, n + m + 1) dev = qml.device("default.qubit", wires=(n + m + 1)) @qml.qnode(dev) def circuit(): qml.templates.QuantumMonteCarlo( probs, func, target_wires=target_wires, estimation_wires=estimation_wires, ) return qml.probs(estimation_wires) phase_estimated = np.argmax(circuit()[:int(N / 2)]) / N The estimated value can be retrieved using the formula :math:\mu = (1-\cos(\pi \theta))/2 >>> (1 - np.cos(np.pi * phase_estimated)) / 2 0.4327096457464369 """ num_wires = AnyWires grad_method = None @classmethod def _unflatten(cls, data, metadata): new_op = cls.__new__(cls) new_op._hyperparameters = dict(metadata[1]) # pylint: disable=protected-access # call operation.__init__ to initialize private properties like _name, _id, _pauli_rep, etc. Operation.__init__(new_op, *data, wires=metadata[0]) return new_op def __init__(self, probs, func, target_wires, estimation_wires, id=None): if isinstance(probs, np.ndarray) and probs.ndim != 1: raise ValueError("The probability distribution must be specified as a flat array") dim_p = len(probs) num_target_wires_ = np.log2(2 * dim_p) num_target_wires = int(num_target_wires_) if not np.allclose(num_target_wires_, num_target_wires): raise ValueError( "The probability distribution must have a length that is a power of two" ) target_wires = list(target_wires) estimation_wires = list(estimation_wires) wires = target_wires + estimation_wires if num_target_wires != len(target_wires): raise ValueError( f"The probability distribution of dimension {dim_p} requires" f" {num_target_wires} target wires" ) self._hyperparameters = {"estimation_wires": estimation_wires, "target_wires": target_wires} A = probs_to_unitary(probs) R = func_to_unitary(func, dim_p) Q = make_Q(A, R) super().__init__(A, R, Q, wires=wires, id=id) @property def num_params(self): return 3
[docs] @staticmethod def compute_decomposition( A, R, Q, wires, estimation_wires, target_wires ): # pylint: disable=arguments-differ,unused-argument r"""Representation of the operator as a product of other operators. .. math:: O = O_1 O_2 \dots O_n. .. seealso:: :meth:~.QuantumMonteCarlo.decomposition. Args: A (array): unitary matrix corresponding to an input probability distribution R (array): unitary that encodes the function applied to the ancilla qubit register Q (array): matrix that encodes the expectation value according to the probability unitary and the function-encoding unitary wires (Any or Iterable[Any]): full set of wires that the operator acts on target_wires (Iterable[Any]): the target wires estimation_wires (Iterable[Any]): the estimation wires Returns: list[.Operator]: decomposition of the operator """ op_list = [ QubitUnitary(A, wires=target_wires[:-1]), QubitUnitary(R, wires=target_wires), qml.templates.QuantumPhaseEstimation( Q, target_wires=target_wires, estimation_wires=estimation_wires ), ] return op_list

Using PennyLane

Release news

Development

API

Internals