Merge branch 'master' into clifford_unpack_fix

doc/source/api-reference/qibo.rst

@@ -1740,6 +1740,12 @@ Classical relative entropy
 .. autofunction:: qibo.quantum_info.classical_relative_entropy
 
 
+Classical mutual information
+""""""""""""""""""""""""""""
+
+.. autofunction:: qibo.quantum_info.classical_mutual_information
+
+
 Classical Rényi entropy
 """""""""""""""""""""""
 
@@ -1785,6 +1791,12 @@ Relative von Neumann entropy
     an error will be raised when using `cupy` backend.
 
 
+Mutual information
+""""""""""""""""""
+
+.. autofunction:: qibo.quantum_info.mutual_information
+
+
 Rényi entropy
 """""""""""""
 

src/qibo/backends/abstract.py

@@ -318,12 +318,16 @@ def calculate_overlap_density_matrix(self, state1, state2):  # pragma: no cover
         raise_error(NotImplementedError)
 
     @abc.abstractmethod
-    def calculate_eigenvalues(self, matrix, k=6):  # pragma: no cover
+    def calculate_eigenvalues(
+        self, matrix, k: int = 6, hermitian: bool = True
+    ):  # pragma: no cover
         """Calculate eigenvalues of a matrix."""
         raise_error(NotImplementedError)
 
     @abc.abstractmethod
-    def calculate_eigenvectors(self, matrix, k=6):  # pragma: no cover
+    def calculate_eigenvectors(
+        self, matrix, k: int = 6, hermitian: bool = True
+    ):  # pragma: no cover
         """Calculate eigenvectors of a matrix."""
         raise_error(NotImplementedError)
 

src/qibo/backends/numpy.py

@@ -719,7 +719,7 @@ def calculate_overlap_density_matrix(self, state1, state2):
             self.np.matmul(self.np.conj(self.cast(state1)).T, self.cast(state2))
         )
 
-    def calculate_eigenvalues(self, matrix, k=6, hermitian=True):
+    def calculate_eigenvalues(self, matrix, k: int = 6, hermitian: bool = True):
         if self.is_sparse(matrix):
             log.warning(
                 "Calculating sparse matrix eigenvectors because "
@@ -730,7 +730,7 @@ def calculate_eigenvalues(self, matrix, k=6, hermitian=True):
             return np.linalg.eigvalsh(matrix)
         return np.linalg.eigvals(matrix)
 
-    def calculate_eigenvectors(self, matrix, k=6, hermitian=True):
+    def calculate_eigenvectors(self, matrix, k: int = 6, hermitian: bool = True):
         if self.is_sparse(matrix):
             if k < matrix.shape[0]:
                 from scipy.sparse.linalg import eigsh

src/qibo/backends/pytorch.py

@@ -180,12 +180,12 @@ def sample_shots(self, probabilities, nshots):
             self.cast(probabilities, dtype="float"), nshots, replacement=True
         )
 
-    def calculate_eigenvalues(self, matrix, k=6, hermitian=True):
+    def calculate_eigenvalues(self, matrix, k: int = 6, hermitian: bool = True):
         if hermitian:
             return self.np.linalg.eigvalsh(matrix)  # pylint: disable=not-callable
         return self.np.linalg.eigvals(matrix)  # pylint: disable=not-callable
 
-    def calculate_eigenvectors(self, matrix, k=6, hermitian=True):
+    def calculate_eigenvectors(self, matrix, k: int = 6, hermitian: int = True):
         if hermitian:
             return self.np.linalg.eigh(matrix)  # pylint: disable=not-callable
         return self.np.linalg.eig(matrix)  # pylint: disable=not-callable

src/qibo/backends/tensorflow.py

@@ -177,12 +177,12 @@ def calculate_norm_density_matrix(self, state, order="nuc"):
             return self.np.trace(state)
         return self.tf.norm(state, ord=order)
 
-    def calculate_eigenvalues(self, matrix, k=6, hermitian=True):
+    def calculate_eigenvalues(self, matrix, k: int = 6, hermitian: bool = True):
         if hermitian:
             return self.tf.linalg.eigvalsh(matrix)
         return self.tf.linalg.eigvals(matrix)
 
-    def calculate_eigenvectors(self, matrix, k=6, hermitian=True):
+    def calculate_eigenvectors(self, matrix, k: int = 6, hermitian: bool = True):
         if hermitian:
             return self.tf.linalg.eigh(matrix)
         return self.tf.linalg.eig(matrix)

src/qibo/quantum_info/entropies.py

@@ -9,7 +9,7 @@
 from qibo.backends.pytorch import PyTorchBackend
 from qibo.config import PRECISION_TOL, raise_error
 from qibo.quantum_info.linalg_operations import partial_trace
-from qibo.quantum_info.metrics import _check_hermitian_or_not_gpu, purity
+from qibo.quantum_info.metrics import _check_hermitian, purity
 
 
 def shannon_entropy(prob_dist, base: float = 2, backend=None):
@@ -30,7 +30,7 @@ def shannon_entropy(prob_dist, base: float = 2, backend=None):
             Defaults to ``None``.
 
     Returns:
-        (float): Shannon entropy :math:`H(\\mathcal{p})`.
+        float: Shannon entropy :math:`H(\\mathcal{p})`.
     """
     backend = _check_backend(backend)
 
@@ -143,6 +143,40 @@ def classical_relative_entropy(prob_dist_p, prob_dist_q, base: float = 2, backen
     return entropy_p - relative
 
 
+def classical_mutual_information(
+    prob_dist_joint, prob_dist_p, prob_dist_q, base: float = 2, backend=None
+):
+    """Calculates the classical mutual information of two random variables.
+
+    Given two random variables :math:`(X, \\, Y)`, their mutual information is given by
+
+    .. math::
+        I(X, \\, Y) \\equiv H(p(x)) + H(q(y)) - H(p(x, \\, y)) \\, ,
+
+    where :math:`p(x, \\, y)` is the joint probability distribution of :math:`(X, Y)`,
+    :math:`p(x)` is the marginal probability distribution of :math:`X`,
+    :math:`q(y)` is the marginal probability distribution of :math:`Y`,
+    and :math:`H(\\cdot)` is the :func:`qibo.quantum_info.entropies.shannon_entropy`.
+
+    Args:
+        prob_dist_joint (ndarray): joint probability distribution :math:`p(x, \\, y)`.
+        prob_dist_p (ndarray): marginal probability distribution :math:`p(x)`.
+        prob_dist_q (ndarray): marginal probability distribution :math:`q(y)`.
+        base (float): the base of the log. 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``.
+
+    Returns:
+        float: Mutual information :math:`I(X, \\, Y)`.
+    """
+    return (
+        shannon_entropy(prob_dist_p, base, backend)
+        + shannon_entropy(prob_dist_q, base, backend)
+        - shannon_entropy(prob_dist_joint, base, backend)
+    )
+
+
 def classical_renyi_entropy(
     prob_dist, alpha: Union[float, int], base: float = 2, backend=None
 ):
@@ -458,9 +492,7 @@ def von_neumann_entropy(
 
     eigenvalues = backend.calculate_eigenvalues(
         state,
-        hermitian=(
-            not check_hermitian or _check_hermitian_or_not_gpu(state, backend=backend)
-        ),
+        hermitian=(not check_hermitian or _check_hermitian(state, backend=backend)),
     )
 
     log_prob = backend.np.where(
@@ -549,15 +581,11 @@ def relative_von_neumann_entropy(
 
     eigenvalues_state = backend.calculate_eigenvalues(
         state,
-        hermitian=(
-            not check_hermitian or _check_hermitian_or_not_gpu(state, backend=backend)
-        ),
+        hermitian=(not check_hermitian or _check_hermitian(state, backend=backend)),
     )
     eigenvalues_target = backend.calculate_eigenvalues(
         target,
-        hermitian=(
-            not check_hermitian or _check_hermitian_or_not_gpu(target, backend=backend)
-        ),
+        hermitian=(not check_hermitian or _check_hermitian(target, backend=backend)),
     )
 
     log_state = backend.np.where(
@@ -580,6 +608,50 @@ def relative_von_neumann_entropy(
     return float(backend.np.real(entropy_state - relative))
 
 
+def mutual_information(
+    state, partition, base: float = 2, check_hermitian: bool = False, backend=None
+):
+    """Calculates the mutual information of a bipartite state.
+
+    Given a qubit ``partition`` :math:`A`, the mutual information
+    of state :math:`\\rho` is given by
+
+    .. math::
+        I(\\rho}) \\equiv S(\\rho_{A}) + S(\\rho_{B}) - S(\\rho) \\, ,
+
+    where :math:`B` is the remaining qubits that are not in partition :math:`A`,
+    and :math:`S(\\cdot)` is the :func:`qibo.quantum_info.von_neumann_entropy`.
+
+    Args:
+        state (ndarray): statevector or density matrix.
+        partition (Union[List[int], Tuple[int]]): indices of qubits in partition :math:`A`.
+        base (float, optional): the base of the log. Defaults to :math:`2`.
+        check_hermitian (bool, optional): if ``True``, checks if ``state`` is Hermitian.
+            If ``False``, it assumes ``state`` is Hermitian . Defaults to ``False``.
+        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:
+        float: Mutual information :math:`I(\\rho)` of ``state`` :math:`\\rho`.
+    """
+    nqubits = np.log2(len(state))
+
+    if not nqubits.is_integer():
+        raise_error(ValueError, f"dimensions of ``state`` must be a power of 2.")
+
+    partition_b = set(list(range(int(nqubits)))) ^ set(list(partition))
+
+    state_a = partial_trace(state, partition_b, backend)
+    state_b = partial_trace(state, partition, backend)
+
+    return (
+        von_neumann_entropy(state_a, base, check_hermitian, False, backend)
+        + von_neumann_entropy(state_b, base, check_hermitian, False, backend)
+        - von_neumann_entropy(state, base, check_hermitian, False, backend)
+    )
+
+
 def renyi_entropy(state, alpha: Union[float, int], base: float = 2, backend=None):
     """Calculates the Rényi entropy :math:`H_{\\alpha}` of a quantum state :math:`\\rho`.
 

src/qibo/quantum_info/metrics.py

@@ -234,7 +234,7 @@ def fidelity(state, target, check_hermitian: bool = False, backend=None):
             abs(purity_state - 1) > PRECISION_TOL
             and abs(purity_target - 1) > PRECISION_TOL
         ):
-            hermitian = check_hermitian is False or _check_hermitian_or_not_gpu(
+            hermitian = check_hermitian is False or _check_hermitian(
                 state, backend=backend
             )
             # using eigh since rho is supposed to be Hermitian
@@ -812,10 +812,8 @@ def frame_potential(
     return potential / samples**2
 
 
-def _check_hermitian_or_not_gpu(matrix, backend=None):
-    """Checks if a given matrix is Hermitian and whether the backend is neither
-    :class:`qibojit.backends.CupyBackend` nor
-    :class:`qibojit.backends.CuQuantumBackend`.
+def _check_hermitian(matrix, backend=None):
+    """Checks if a given matrix is Hermitian.
 
     Args:
         matrix: input array.
@@ -825,10 +823,6 @@ def _check_hermitian_or_not_gpu(matrix, backend=None):
 
     Returns:
         bool: whether the matrix is Hermitian.
-
-    Raises:
-        NotImplementedError: If `matrix` is not Hermitian and
-        `backend` is not :class:`qibojit.backends.CupyBackend`
     """
     backend = _check_backend(backend)
 
@@ -838,14 +832,4 @@ def _check_hermitian_or_not_gpu(matrix, backend=None):
 
     hermitian = bool(float(norm) <= PRECISION_TOL)
 
-    if hermitian is False and backend.__class__.__name__ in [
-        "CupyBackend",
-        "CuQuantumBackend",
-    ]:  # pragma: no cover
-        raise_error(
-            NotImplementedError,
-            "GPU backends do not support `np.linalg.eig` "
-            + "or `np.linalg.eigvals` for non-Hermitian matrices.",
-        )
-
     return hermitian

tests/test_quantum_info_entropies.py

@@ -5,11 +5,13 @@
 from qibo.config import PRECISION_TOL
 from qibo.quantum_info.entropies import (
     _matrix_power,
+    classical_mutual_information,
     classical_relative_entropy,
     classical_relative_renyi_entropy,
     classical_renyi_entropy,
     classical_tsallis_entropy,
     entanglement_entropy,
+    mutual_information,
     relative_renyi_entropy,
     relative_von_neumann_entropy,
     renyi_entropy,
@@ -125,6 +127,27 @@ def test_classical_relative_entropy(backend, base, kind):
     backend.assert_allclose(divergence, target, atol=1e-5)
 
 
+@pytest.mark.parametrize("base", [2, 10, np.e, 5])
+def test_classical_mutual_information(backend, base):
+    prob_p = np.random.rand(10)
+    prob_q = np.random.rand(10)
+    prob_p /= np.sum(prob_p)
+    prob_q /= np.sum(prob_q)
+
+    joint_dist = np.kron(prob_p, prob_q)
+    joint_dist /= np.sum(joint_dist)
+
+    prob_p = backend.cast(prob_p, dtype=prob_p.dtype)
+    prob_q = backend.cast(prob_q, dtype=prob_q.dtype)
+    joint_dist = backend.cast(joint_dist, dtype=joint_dist.dtype)
+
+    backend.assert_allclose(
+        classical_mutual_information(joint_dist, prob_p, prob_q, base, backend),
+        0.0,
+        atol=1e-10,
+    )
+
+
 @pytest.mark.parametrize("kind", [None, list])
 @pytest.mark.parametrize("base", [2, 10, np.e, 5])
 @pytest.mark.parametrize("alpha", [0, 1, 2, 3, np.inf])
@@ -499,6 +522,25 @@ def test_relative_entropy(backend, base, check_hermitian):
         )
 
 
+@pytest.mark.parametrize("check_hermitian", [False, True])
+@pytest.mark.parametrize("base", [2, 10, np.e, 5])
+def test_mutual_information(backend, base, check_hermitian):
+    with pytest.raises(ValueError):
+        state = np.ones((3, 3))
+        state = backend.cast(state, dtype=state.dtype)
+        test = mutual_information(state, [0], backend)
+
+    state_a = random_density_matrix(4, backend=backend)
+    state_b = random_density_matrix(4, backend=backend)
+    state = backend.np.kron(state_a, state_b)
+
+    backend.assert_allclose(
+        mutual_information(state, [0, 1], base, check_hermitian, backend),
+        0.0,
+        atol=1e-10,
+    )
+
+
 @pytest.mark.parametrize("base", [2, 10, np.e, 5])
 @pytest.mark.parametrize("alpha", [0, 1, 2, 3, np.inf])
 def test_renyi_entropy(backend, alpha, base):