Unary encoder for random Gaussian states

doc/source/api-reference/qibo.rst

@@ -236,6 +236,12 @@ Unary Encoder
 .. autofunction:: qibo.models.encodings.unary_encoder
 
 
+Unary Encoder for Random Gaussian States
+""""""""""""""""""""""""""""""""""""""""
+
+.. autofunction:: qibo.models.encodings.unary_encoder_random_gaussian
+
+
 .. _error-mitigation:
 
 Error Mitigation

src/qibo/models/encodings.py

@@ -3,6 +3,7 @@
 import math
 
 import numpy as np
+from scipy.stats import rv_continuous
 
 from qibo import gates
 from qibo.config import raise_error
@@ -73,6 +74,86 @@ def unary_encoder(data):
     return circuit
 
 
+def unary_encoder_random_gaussian(nqubits: int, seed=None):
+    """Creates a circuit that performs the unary encoding of a random Gaussian state.
+
+    Given :math:`d` qubits, encodes the quantum state
+    :math:`\\ket{\\psi} \\in \\mathcal{H}` such that
+
+
+    .. math::
+        \\ket{\\psi} = \\frac{1}{\\|\\mathbf{x}\\|_{\\textup{HS}}} \\,
+            \\sum_{k=1}^{d} \\, x_{k} \\, \\ket{k}
+
+    where :math:`x_{k}` are independent Gaussian random variables,
+    :math:`\\mathcal{H} \\cong \\mathbb{C}^{d}` is a :math:`d`-qubit Hilbert space,
+    and :math:`\\|\\cdot\\|_{\\textup{HS}}` being the Hilbert-Schmidt norm.
+    Here, :math:`\\ket{k}` is a unary representation of the number :math:`1` through
+    :math:`d`.
+
+    At depth :math:`h`, the angles :math:`\\theta_{k} \\in [0, 2\\pi]` of the the
+    gates :math:`RBS(\\theta_{k})` are sampled from the following probability density function:
+
+    .. math::
+        p_{h}(\\theta) = \\frac{1}{2} \\, \\frac{\\Gamma(2^{h-1})}{\\Gamma^{2}(2^{h-2})}
+            \\abs{\\sin(\\theta) \\, \\cos(\\theta)}^{2^{h-1} - 1} \\, ,
+
+    where :math:`\\Gamma(\\cdot)` is the
+    `Gamma function <https://en.wikipedia.org/wiki/Gamma_function>`_.
+
+    Args:
+        nqubits (int): number of qubits.
+        seed (int or :class:`numpy.random.Generator`, optional): Either a generator of
+            random numbers or a fixed seed to initialize a generator. If ``None``,
+            initializes a generator with a random seed. Defaults to ``None``.
+
+    Returns:
+        :class:`qibo.models.circuit.Circuit`: circuit that loads a random Gaussian array in unary representation.
+
+    References:
+        1. A. Bouland, A. Dandapani, and A. Prakash, *A quantum spectral method for simulating
+        stochastic processes, with applications to Monte Carlo*.
+        `arXiv:2303.06719v1 [quant-ph] <https://arxiv.org/abs/2303.06719>`_
+    """
+    if not isinstance(nqubits, int):
+        raise_error(
+            TypeError, f"nqubits must be type int, but it is type {type(nqubits)}."
+        )
+    elif nqubits <= 0.0:
+        raise_error(
+            ValueError, f"nqubits must be a positive integer, but it is {nqubits}."
+        )
+    elif not math.log2(nqubits).is_integer():
+        raise_error(ValueError, f"nqubits must be a power of 2, but it is {nqubits}.")
+
+    if (
+        seed is not None
+        and not isinstance(seed, int)
+        and not isinstance(seed, np.random.Generator)
+    ):
+        raise_error(
+            TypeError, "seed must be either type int or numpy.random.Generator."
+        )
+
+    local_state = (
+        np.random.default_rng(seed) if seed is None or isinstance(seed, int) else seed
+    )
+
+    sampler = _probability_distribution_gaussian_loader(
+        a=0, b=2 * math.pi, seed=local_state
+    )
+
+    circuit, pairs_rbs = _generate_rbs_pairs(nqubits)
+
+    phases = []
+    for depth, row in enumerate(pairs_rbs, 1):
+        phases += list(sampler.rvs(depth=depth, size=len(row)))
+
+    circuit.set_parameters(phases)
+
+    return circuit
+
+
 def _generate_rbs_pairs(nqubits):
     """Generating list of indexes representing the RBS connections
     and creating circuit with all RBS initialised with 0.0 phase."""
@@ -92,3 +173,15 @@ def _generate_rbs_pairs(nqubits):
             circuit.add(gates.RBS(*pair, 0.0, trainable=True))
 
     return circuit, pairs_rbs
+
+
+class _probability_distribution_gaussian_loader(rv_continuous):
+    """Probability density function for sampling phases of
+    the RBS gates as a function of circuit depth."""
+
+    def _pdf(self, theta: float, depth: int):
+        amplitude = 2 * math.gamma(2 ** (depth - 1)) / math.gamma(2 ** (depth - 2)) ** 2
+
+        probability = abs(math.sin(theta) * math.cos(theta)) ** (2 ** (depth - 1) - 1)
+
+        return amplitude * probability / 4

tests/test_models_encodings.py

@@ -1,8 +1,16 @@
 """Tests for qibo.models.encodings"""
+import math
+
 import numpy as np
 import pytest
+from scipy.optimize import curve_fit
+
+from qibo.models.encodings import unary_encoder, unary_encoder_random_gaussian
+
 
-from qibo.models.encodings import unary_encoder
+def gaussian(x, a, b, c):
+    """Gaussian used in the `unary_encoder_random_gaussian test"""
+    return np.exp(a * x**2 + b * x + c)
 
 
 @pytest.mark.parametrize("nqubits", [2, 4, 8])
@@ -31,3 +39,47 @@ def test_unary_encoder(backend, nqubits):
     backend.assert_allclose(
         state, np.sort(data) / backend.calculate_norm(data, order=2)
     )
+
+
+@pytest.mark.parametrize("seed", [None, 10, np.random.default_rng(10)])
+@pytest.mark.parametrize("nqubits", [8, 16])
+def test_unary_encoder_random_gaussian(backend, nqubits, seed):
+    """Tests if encoded vector are random variables sampled from
+    Gaussian distribution with 0.0 mean and variance close to the norm
+    of the random Gaussian vector that was encoded."""
+    with pytest.raises(TypeError):
+        unary_encoder_random_gaussian("1", seed=seed)
+    with pytest.raises(ValueError):
+        unary_encoder_random_gaussian(-1, seed=seed)
+    with pytest.raises(ValueError):
+        unary_encoder_random_gaussian(3, seed=seed)
+    with pytest.raises(TypeError):
+        unary_encoder_random_gaussian(nqubits, seed="seed")
+
+    samples = int(1e2)
+
+    local_state = np.random.default_rng(seed) if seed in [None, 10] else seed
+
+    amplitudes = []
+    for _ in range(samples):
+        circuit = unary_encoder_random_gaussian(nqubits, seed=local_state)
+        state = backend.execute_circuit(circuit).state()
+        indexes = np.flatnonzero(state)
+        state = np.real(state[indexes])
+        amplitudes += [float(elem) for elem in list(state)]
+
+    y, x = np.histogram(amplitudes, bins=50, density=True)
+    x = (x[:-1] + x[1:]) / 2
+
+    params, _ = curve_fit(gaussian, x, y)
+
+    stddev = np.sqrt(-1 / (2 * params[0]))
+    mean = stddev**2 * params[1]
+
+    theoretical_norm = (
+        math.sqrt(2) * math.gamma((nqubits + 1) / 2) / math.gamma(nqubits / 2)
+    )
+    theoretical_norm = 1.0 / theoretical_norm
+
+    backend.assert_allclose(0.0, mean, atol=1e-1)
+    backend.assert_allclose(stddev, theoretical_norm, atol=1e-1)