Merge pull request #1052 from qiboteam/haar_integral

Calculation of `quantum_info.haar_integral` exactly

src/qibo/models/error_mitigation.py

@@ -1,5 +1,7 @@
 """Error Mitigation Methods."""
 
+from math import factorial
+
 import numpy as np
 from scipy.optimize import curve_fit
 
@@ -34,11 +36,11 @@ def get_gammas(c, solve: bool = True):
         gammas = np.array(
             [
                 1
-                / (2 ** (2 * cmax) * np.math.factorial(i))
+                / (2 ** (2 * cmax) * factorial(i))
                 * (-1) ** i
                 / (1 + 2 * i)
-                * np.math.factorial(1 + 2 * cmax)
-                / (np.math.factorial(cmax) * np.math.factorial(cmax - i))
+                * factorial(1 + 2 * cmax)
+                / (factorial(cmax) * factorial(cmax - i))
                 for i in c
             ]
         )

src/qibo/quantum_info/metrics.py

@@ -1,5 +1,7 @@
 """Submodule with distances, metrics, and measures for quantum states and channels."""
 
+from typing import Optional, Union
+
 import numpy as np
 from scipy import sparse
 
@@ -1086,8 +1088,14 @@ def entangling_capability(circuit, samples: int, backend=None):
     return capability
 
 
-def expressibility(circuit, t: int, samples: int, backend=None):
-    """Returns the expressibility :math:`\\|A\\|_{HS}` of a parametrized
+def expressibility(
+    circuit,
+    power_t: int,
+    samples: int,
+    order: Optional[Union[int, float, str]] = 2,
+    backend=None,
+):
+    """Returns the expressibility :math:`\\|A\\|` of a parametrized
     circuit, where
 
     .. math::
@@ -1098,8 +1106,11 @@ def expressibility(circuit, t: int, samples: int, backend=None):
 
     Args:
         circuit (:class:`qibo.models.Circuit`): Parametrized circuit.
-        t (int): power that defines the :math:`t`-design.
+        power_t (int): power that defines the :math:`t`-design.
         samples (int): number of samples to estimate the integrals.
+        order (int or float or str, optional): order of the norm :math:`\\|A\\|`.
+            For specifications, see :meth:`qibo.backends.abstract.calculate_norm`.
+            Defaults to :math:`2`.
         backend (:class:`qibo.backends.abstract.Backend`, optional): backend to be used
             in the execution. If ``None``, it uses :class:`qibo.backends.GlobalBackend`.
             Defaults to ``None``.
@@ -1108,8 +1119,10 @@ def expressibility(circuit, t: int, samples: int, backend=None):
         float: Entangling capability.
     """
 
-    if isinstance(t, int) is False:
-        raise_error(TypeError, f"t must be type int, but it is type {type(t)}.")
+    if isinstance(power_t, int) is False:
+        raise_error(
+            TypeError, f"power_t must be type int, but it is type {type(power_t)}."
+        )
 
     if isinstance(samples, int) is False:
         raise_error(
@@ -1124,11 +1137,11 @@ def expressibility(circuit, t: int, samples: int, backend=None):
     if backend is None:  # pragma: no cover
         backend = GlobalBackend()
 
-    expr = haar_integral(circuit.nqubits, t, samples, backend=backend) - pqc_integral(
-        circuit, t, samples, backend=backend
-    )
+    deviation = haar_integral(
+        circuit.nqubits, power_t, samples=None, backend=backend
+    ) - pqc_integral(circuit, power_t, samples, backend=backend)
 
-    fid = np.trace(expr @ expr)
+    fid = float(backend.calculate_norm(deviation, order=order))
 
     return fid
 

src/qibo/quantum_info/utils.py

@@ -1,7 +1,10 @@
 """Utility functions for the Quantum Information module."""
 
 from functools import reduce
+from itertools import permutations
+from math import factorial
 from re import finditer
+from typing import Optional
 
 import numpy as np
 
@@ -276,7 +279,12 @@ def hellinger_fidelity(prob_dist_p, prob_dist_q, validate: bool = False, backend
     return (1 - distance**2) ** 2
 
 
-def haar_integral(nqubits: int, power_t: int, samples: int, backend=None):
+def haar_integral(
+    nqubits: int,
+    power_t: int,
+    samples: Optional[int] = None,
+    backend=None,
+):
     """Returns the integral over pure states over the Haar measure.
 
     .. math::
@@ -286,13 +294,19 @@ def haar_integral(nqubits: int, power_t: int, samples: int, backend=None):
     Args:
         nqubits (int): Number of qubits.
         power_t (int): power that defines the :math:`t`-design.
-        samples (int): number of samples to estimate the integral.
+        samples (int, optional): If ``None``, estimated the integral exactly.
+            Otherwise, number of samples to estimate the integral via sampling.
+            Defaults to ``None``.
         backend (:class:`qibo.backends.abstract.Backend`, optional): backend to be
             used in the execution. If ``None``, it uses
             :class:`qibo.backends.GlobalBackend`. Defaults to ``None``.
 
     Returns:
         array: Estimation of the Haar integral.
+
+    .. note::
+        The ``exact=True`` method is implemented using Lemma 34 of
+        `Kliesch and Roth (2020) <https://arxiv.org/abs/2010.05925>`_.
     """
 
     if isinstance(nqubits, int) is False:
@@ -305,32 +319,54 @@ def haar_integral(nqubits: int, power_t: int, samples: int, backend=None):
             TypeError, f"power_t must be type int, but it is type {type(power_t)}."
         )
 
-    if isinstance(samples, int) is False:
+    if samples is not None and isinstance(samples, int) is False:
         raise_error(
             TypeError, f"samples must be type int, but it is type {type(samples)}."
         )
 
-    from qibo.quantum_info.random_ensembles import (  # pylint: disable=C0415
-        random_statevector,
-    )
-
     if backend is None:  # pragma: no cover
         backend = GlobalBackend()
 
     dim = 2**nqubits
 
-    rand_unit_density = np.zeros((dim**power_t, dim**power_t), dtype=complex)
-    rand_unit_density = backend.cast(rand_unit_density, dtype=rand_unit_density.dtype)
-    for _ in range(samples):
-        haar_state = np.reshape(
-            random_statevector(dim, haar=True, backend=backend), (-1, 1)
+    if samples is not None:
+        from qibo.quantum_info.random_ensembles import (  # pylint: disable=C0415
+            random_statevector,
         )
 
-        rho = haar_state @ np.conj(np.transpose(haar_state))
+        rand_unit_density = np.zeros((dim**power_t, dim**power_t), dtype=complex)
+        rand_unit_density = backend.cast(
+            rand_unit_density, dtype=rand_unit_density.dtype
+        )
+        for _ in range(samples):
+            haar_state = np.reshape(
+                random_statevector(dim, haar=True, backend=backend), (-1, 1)
+            )
 
-        rand_unit_density += reduce(np.kron, [rho] * power_t)
+            rho = haar_state @ np.conj(np.transpose(haar_state))
 
-    integral = rand_unit_density / samples
+            rand_unit_density += reduce(np.kron, [rho] * power_t)
+
+        integral = rand_unit_density / samples
+
+        return integral
+
+    normalization = factorial(dim - 1) / factorial(dim - 1 + power_t)
+
+    permutations_list = list(permutations(np.arange(power_t) + power_t))
+    permutations_list = [
+        tuple(np.arange(power_t)) + indices for indices in permutations_list
+    ]
+
+    identity = np.eye(dim**power_t, dtype=float)
+    identity = backend.cast(identity, dtype=identity.dtype)
+    identity = np.reshape(identity, (dim,) * (2 * power_t))
+
+    integral = np.zeros((dim**power_t, dim**power_t), dtype=float)
+    integral = backend.cast(integral, dtype=integral.dtype)
+    for indices in permutations_list:
+        integral += np.reshape(np.transpose(identity, indices), (-1, dim**power_t))
+    integral *= normalization
 
     return integral
 

tests/test_quantum_info_utils.py

@@ -168,25 +168,28 @@ def test_hellinger(backend):
     backend.assert_allclose(hellinger_fidelity(prob, prob_q, backend=backend), 1.0)
 
 
-def test_haar_integral(backend):
+def test_haar_integral_errors(backend):
     with pytest.raises(TypeError):
         nqubits, power_t, samples = 0.5, 2, 10
         test = haar_integral(nqubits, power_t, samples, backend=backend)
     with pytest.raises(TypeError):
         nqubits, power_t, samples = 2, 0.5, 10
         test = haar_integral(nqubits, power_t, samples, backend=backend)
     with pytest.raises(TypeError):
-        nqubits, power_t, samples = 2, 2, 0.5
-        test = haar_integral(nqubits, power_t, samples, backend=backend)
+        nqubits, power_t, samples = 2, 1, 1.2
+        test = haar_integral(nqubits, power_t, samples=samples, backend=backend)
+
 
-    nqubits = 2
-    power_t, samples = 1, 1000
+@pytest.mark.parametrize("power_t", [1, 2, 3])
+@pytest.mark.parametrize("nqubits", [2, 3])
+def test_haar_integral(backend, nqubits, power_t):
+    samples = int(1e4)
 
-    haar_int = haar_integral(nqubits, power_t, samples, backend=backend)
+    haar_int_exact = haar_integral(nqubits, power_t, samples=None, backend=backend)
 
-    fid = fidelity(haar_int, haar_int, backend=backend)
+    haar_int_sampled = haar_integral(nqubits, power_t, samples=samples, backend=backend)
 
-    backend.assert_allclose(fid, 1.0, atol=10 / samples)
+    backend.assert_allclose(haar_int_sampled, haar_int_exact, atol=1e-2)
 
 
 def test_pqc_integral(backend):