diff --git a/doc/source/api-reference/qibo.rst b/doc/source/api-reference/qibo.rst index ff8308a5ad..4e67d05199 100644 --- a/doc/source/api-reference/qibo.rst +++ b/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 diff --git a/src/qibo/models/encodings.py b/src/qibo/models/encodings.py index bb2ecf298c..9621d80fcd 100644 --- a/src/qibo/models/encodings.py +++ b/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 `_. + + 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] `_ + """ + 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 = _ProbabilityDistributionGaussianLoader( + 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.extend(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.""" @@ -85,10 +166,26 @@ def _generate_rbs_pairs(nqubits): pairs_rbs += pairs_rbs_per_depth indexes = list(np.array(pairs_rbs_per_depth).flatten()) + pairs_rbs = [ + [(nqubits - 1 - a, nqubits - 1 - b) for a, b in row] for row in pairs_rbs + ] + circuit = Circuit(nqubits) - circuit.add(gates.X(0)) + circuit.add(gates.X(nqubits - 1)) for row in pairs_rbs: for pair in row: circuit.add(gates.RBS(*pair, 0.0, trainable=True)) return circuit, pairs_rbs + + +class _ProbabilityDistributionGaussianLoader(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 diff --git a/tests/test_models_encodings.py b/tests/test_models_encodings.py index 3991c88c39..d0c255ae49 100644 --- a/tests/test_models_encodings.py +++ b/tests/test_models_encodings.py @@ -1,11 +1,19 @@ """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]) +@pytest.mark.parametrize("nqubits", [2, 4, 8, 16]) def test_unary_encoder(backend, nqubits): sampler = np.random.default_rng(1) @@ -26,8 +34,50 @@ def test_unary_encoder(backend, nqubits): circuit = unary_encoder(data) state = backend.execute_circuit(circuit).state() indexes = np.flatnonzero(state) - state = np.sort(state[indexes]) + state = state[indexes] + + backend.assert_allclose(state, data / backend.calculate_norm(data, order=2)) + - 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)