This project is archived
Qiskit Finance is no longer maintained and this documentation is no longer updated. It is kept online for reference only, and the code it describes may not work with current versions of Qiskit. See the project repository for more information.
Source code for qiskit_finance.circuit.library.probability_distributions.normal
# This code is part of a Qiskit project.
#
# (C) Copyright IBM 2017, 2024.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.
"""A circuit that encodes a discretized normal probability distribution in qubit amplitudes."""
from typing import Tuple, Union, List, Optional, Any
import numpy as np
from qiskit.circuit import QuantumCircuit
from qiskit.circuit.library import Initialize, Isometry
[docs]class NormalDistribution(QuantumCircuit):
r"""A circuit to encode a discretized normal distribution in qubit amplitudes.
The probability density function of the normal distribution is defined as
.. math::
\mathbb{P}(X = x) = \frac{1}{\sqrt{2\pi\sigma^2}} e^{-\frac{(x - \mu)^2}{\sigma^2}}
.. note::
The parameter ``sigma`` in this class equals the **variance**, :math:`\sigma^2` and not the
standard deviation. This is for consistency with multivariate distributions, where the
uppercase sigma, :math:`\Sigma`, is associated with the covariance.
This circuit considers the discretized version of the normal distribution on
``2 ** num_qubits`` equidistant points, :math:`x_i`, truncated to ``bounds``.
For a one-dimensional random variable, meaning `num_qubits` is a single integer, it applies
the operation
.. math::
\mathcal{P}_X |0\rangle^n = \sum_{i=0}^{2^n - 1} \sqrt{\mathbb{P}(x_i)} |i\rangle
where :math:`n` is `num_qubits`.
.. note::
The circuit loads the **square root** of the probabilities into the qubit amplitudes such
that the sampling probability, which is the square of the amplitude, equals the
probability of the distribution.
In the multi-dimensional case, the distribution is defined as
.. math::
\mathbb{P}(X = x) = \frac{\Sigma^{-1}}{\sqrt{2\pi}} e^{-\frac{(x - \mu)^2}{\Sigma}}
where :math:`\Sigma` is the covariance. To specify a multivariate normal distribution,
``num_qubits`` is a list of integers, each specifying how many
qubits are used to discretize the respective dimension. The arguments ``mu`` and ``sigma``
in this case are a vector and square matrix.
If for instance, ``num_qubits = [2, 3]`` then ``mu`` is a 2d vector and ``sigma`` is the
:math:`2 \times 2` covariance matrix. The first dimension is discretized using 2 qubits, hence
on 4 points, and the second dimension on 3 qubits, hence 8 points. Therefore the random variable
is discretized on :math:`4 \times 8 = 32` points.
Since, in general, it is not yet known how to efficiently prepare the qubit amplitudes to
represent a normal distribution, this class computes the expected amplitudes and then uses
the ``QuantumCircuit.initialize`` method to construct the corresponding circuit.
This circuit is for example used in amplitude estimation applications, such as finance [1, 2],
where customer demand or the return of a portfolio could be modeled using a normal
distribution.
Examples:
>>> from qiskit_finance.circuit.library.probability_distributions import NormalDistribution
>>> circuit = NormalDistribution(3, mu=1, sigma=1, bounds=(0, 2))
>>> circuit.decompose().draw()
»
q_0: ───────────────────────────────────────────────────»
┌────────────────────────┐»
q_1: ─────────────────────────┤0 ├»
┌───────────────────────┐│ multiplex2_reverse_dg │»
q_2: ┤ multiplex1_reverse_dg ├┤1 ├»
└───────────────────────┘└────────────────────────┘»
« ┌────────────────────────┐
«q_0: ┤0 ├
« │ │
«q_1: ┤1 multiplex3_reverse_dg ├
« │ │
«q_2: ┤2 ├
« └────────────────────────┘
>>> mu = [1, 0.9]
>>> sigma = [[1, -0.2], [-0.2, 1]]
>>> circuit = NormalDistribution([2, 3], mu, sigma)
>>> circuit.num_qubits
5
>>> import os
>>> from qiskit import QuantumCircuit
>>> mu = [1, 0.9]
>>> sigma = [[1, -0.2], [-0.2, 1]]
>>> bounds = [(0, 1), (-1, 1)]
>>> p_x = NormalDistribution([2, 3], mu, sigma, bounds)
>>> circuit = QuantumCircuit(6)
>>> _ = circuit.append(p_x, list(range(5)))
>>> for i in range(5):
... _ = circuit.cry(2 ** i, i, 5)
>>> # strip trailing white spaces to match the expected output
>>> print(os.linesep.join([s.rstrip() for s in circuit.draw().lines()]))
┌───────┐
q_0: ┤0 ├────■─────────────────────────────────────────
│ │ │
q_1: ┤1 ├────┼────────■────────────────────────────────
│ │ │ │
q_2: ┤2 P(X) ├────┼────────┼────────■───────────────────────
│ │ │ │ │
q_3: ┤3 ├────┼────────┼────────┼────────■──────────────
│ │ │ │ │ │
q_4: ┤4 ├────┼────────┼────────┼────────┼────────■─────
└───────┘┌───┴───┐┌───┴───┐┌───┴───┐┌───┴───┐┌───┴────┐
q_5: ─────────┤ Ry(1) ├┤ Ry(2) ├┤ Ry(4) ├┤ Ry(8) ├┤ Ry(16) ├
└───────┘└───────┘└───────┘└───────┘└────────┘
References:
[1]: Gacon, J., Zoufal, C., & Woerner, S. (2020).
Quantum-Enhanced Simulation-Based Optimization.
`arXiv:2005.10780 <http://arxiv.org/abs/2005.10780>`_
[2]: Woerner, S., & Egger, D. J. (2018).
Quantum Risk Analysis.
`arXiv:1806.06893 <http://arxiv.org/abs/1806.06893>`_
"""
def __init__(
self,
num_qubits: Union[int, List[int]],
mu: Optional[Union[float, List[float]]] = None,
sigma: Optional[Union[float, List[float]]] = None,
bounds: Optional[Union[Tuple[float, float], List[Tuple[float, float]]]] = None,
upto_diag: bool = False,
name: str = "P(X)",
) -> None:
r"""
Args:
num_qubits: The number of qubits used to discretize the random variable. For a 1d
random variable, ``num_qubits`` is an integer, for multiple dimensions a list
of integers indicating the number of qubits to use in each dimension.
mu: The parameter :math:`\mu`, which is the expected value of the distribution.
Can be either a float for a 1d random variable or a list of floats for a higher
dimensional random variable. Defaults to 0.
sigma: The parameter :math:`\sigma^2` or :math:`\Sigma`, which is the variance or
covariance matrix. Default to the identity matrix of appropriate size.
bounds: The truncation bounds of the distribution as tuples. For multiple dimensions,
``bounds`` is a list of tuples ``[(low0, high0), (low1, high1), ...]``.
If ``None``, the bounds are set to ``(-1, 1)`` for each dimension.
upto_diag: If True, load the square root of the probabilities up to multiplication
with a diagonal for a more efficient circuit.
name: The name of the circuit.
"""
_check_dimensions_match(num_qubits, mu, sigma, bounds)
_check_bounds_valid(bounds)
# set default arguments
dim = 1 if isinstance(num_qubits, int) else len(num_qubits)
if mu is None:
mu = 0 if dim == 1 else [0] * dim
if sigma is None:
sigma = 1 if dim == 1 else np.eye(dim) # type: ignore[assignment]
if bounds is None:
bounds = (-1, 1) if dim == 1 else [(-1, 1)] * dim
if isinstance(num_qubits, int): # univariate case
inner = QuantumCircuit(num_qubits, name=name)
x = np.linspace(bounds[0], bounds[1], num=2**num_qubits) # type: Any
else: # multivariate case
inner = QuantumCircuit(sum(num_qubits), name=name)
# compute the evaluation points using numpy.meshgrid
# indexing 'ij' yields the "column-based" indexing
meshgrid = np.meshgrid(
*[
np.linspace(bound[0], bound[1], num=2 ** num_qubits[i]) # type: ignore
for i, bound in enumerate(bounds)
],
indexing="ij",
)
# flatten into a list of points
x = list(zip(*[grid.flatten() for grid in meshgrid]))
from scipy.stats import multivariate_normal
# compute the normalized, truncated probabilities
probabilities = multivariate_normal.pdf(x, mu, sigma)
normalized_probabilities = probabilities / np.sum(probabilities)
# store the values, probabilities and bounds to make them user accessible
self._values = x
self._probabilities = normalized_probabilities
self._bounds = bounds
super().__init__(*inner.qregs, name=name)
# use default the isometry (or initialize w/o resets) algorithm to construct the circuit
if upto_diag:
inner.append(Isometry(np.sqrt(normalized_probabilities), 0, 0), inner.qubits)
self.append(inner.to_instruction(), inner.qubits) # Isometry is not a Gate
else:
initialize = Initialize(np.sqrt(normalized_probabilities))
circuit = initialize.gates_to_uncompute().inverse()
inner.compose(circuit, inplace=True)
self.append(inner.to_gate(), inner.qubits)
@property
def values(self) -> np.ndarray:
"""Return the discretized points of the random variable."""
return self._values
@property
def probabilities(self) -> np.ndarray:
"""Return the sampling probabilities for the values."""
return self._probabilities
@property
def bounds(self) -> Union[Tuple[float, float], List[Tuple[float, float]]]:
"""Return the bounds of the probability distribution."""
return self._bounds
def _check_dimensions_match(num_qubits, mu, sigma, bounds):
num_qubits = [num_qubits] if not isinstance(num_qubits, (list, np.ndarray)) else num_qubits
dim = len(num_qubits)
if mu is not None:
mu = [mu] if not isinstance(mu, (list, np.ndarray)) else mu
if len(mu) != dim:
raise ValueError(
f"Dimension of mu ({len(mu)}) does not match the dimension of the "
f"random variable specified by the number of qubits ({dim})"
)
if sigma is not None:
sigma = [[sigma]] if not isinstance(sigma, (list, np.ndarray)) else sigma
if len(sigma) != dim or len(sigma[0]) != dim:
raise ValueError(
f"Dimension of sigma ({len(sigma)} x {len(sigma[0])}) does not match the dimension of "
f"the random variable specified by the number of qubits ({dim})"
)
if bounds is not None:
# bit differently to cover the case the users might pass `bounds` as a single list,
# e.g. [0, 1], instead of a tuple
bounds = [bounds] if not isinstance(bounds[0], tuple) else bounds
if len(bounds) != dim:
raise ValueError(
f"Dimension of bounds ({len(bounds)}) does not match the dimension of the "
f"random variable specified by the number of qubits ({dim})"
)
def _check_bounds_valid(bounds):
if bounds is None:
return
bounds = [bounds] if not isinstance(bounds[0], tuple) else bounds
for i, bound in enumerate(bounds):
if not bound[1] - bound[0] > 0:
raise ValueError(
f"Dimension {i} of the bounds are invalid, must be a non-empty "
"interval where the lower bounds is smaller than the upper bound."
)