Merge pull request #1470 from qiboteam/tsallis_quantum

Add `relative_tsallis_entropy` to `quantum_info.entropies` + fix to `relative_von_neumann_entropy`

doc/source/api-reference/qibo.rst

@@ -1826,6 +1826,12 @@ Tsallis entropy
 .. autofunction:: qibo.quantum_info.tsallis_entropy
 
 
+Relative Tsallis entropy
+""""""""""""""""""""""""
+
+.. autofunction:: qibo.quantum_info.relative_tsallis_entropy
+
+
 Entanglement entropy
 """"""""""""""""""""
 

src/qibo/quantum_info/entropies.py

@@ -561,23 +561,33 @@ def von_neumann_entropy(
 
 
 def relative_von_neumann_entropy(
-    state, target, base: float = 2, check_hermitian: bool = False, backend=None
+    state,
+    target,
+    base: float = 2,
+    check_hermitian: bool = False,
+    precision_tol: float = 1e-14,
+    backend=None,
 ):
-    """Calculates the relative entropy :math:`S(\\rho \\, \\| \\, \\sigma)` between ``state`` :math:`\\rho` and ``target`` :math:`\\sigma`.
+    """Calculates the relative von Neumann entropy  between two quantum states.
 
-    It is given by
+    Also known as *quantum relative entropy*, :math:`S(\\rho \\, \\| \\, \\sigma)` is given by
 
     .. math::
         S(\\rho \\, \\| \\, \\sigma) = \\text{tr}\\left[\\rho \\, \\log(\\rho)\\right]
             - \\text{tr}\\left[\\rho \\, \\log(\\sigma)\\right]
 
+    where ``state`` :math:`\\rho` and ``target`` :math:`\\sigma` are two quantum states.
+
     Args:
         state (ndarray): statevector or density matrix :math:`\\rho`.
         target (ndarray): statevector or density matrix :math:`\\sigma`.
         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``.
+        precision_tol (float, optional): Used when entropy is calculated via engenvalue
+            decomposition. Eigenvalues that are smaller than ``precision_tol`` in absolute value
+            are set to :math:`0`. Defaults to :math:`10^{-14}`.
         backend (:class:`qibo.backends.abstract.Backend`, optional): backend to be used
             in the execution. If ``None``, it uses
             the current backend. Defaults to ``None``.
@@ -627,26 +637,31 @@ def relative_von_neumann_entropy(
     if len(target.shape) == 1:
         target = backend.np.outer(target, backend.np.conj(target))
 
-    eigenvalues_state = backend.calculate_eigenvalues(
+    eigenvalues_state, eigenvectors_state = backend.calculate_eigenvectors(
         state,
         hermitian=(not check_hermitian or _check_hermitian(state, backend=backend)),
     )
-    eigenvalues_target = backend.calculate_eigenvalues(
+    eigenvalues_target, eigenvectors_target = backend.calculate_eigenvectors(
         target,
         hermitian=(not check_hermitian or _check_hermitian(target, backend=backend)),
     )
 
+    overlaps = backend.np.conj(eigenvectors_state.T) @ eigenvectors_target
+    overlaps = backend.np.abs(overlaps) ** 2
+
     log_state = backend.np.where(
-        backend.np.real(eigenvalues_state) > 0,
+        backend.np.real(eigenvalues_state) > precision_tol,
         backend.np.log2(eigenvalues_state) / np.log2(base),
         0.0,
     )
     log_target = backend.np.where(
-        backend.np.real(eigenvalues_target) > 0,
+        backend.np.real(eigenvalues_target) > precision_tol,
         backend.np.log2(eigenvalues_target) / np.log2(base),
-        -np.inf,
+        0.0,
     )
 
+    log_target = overlaps @ log_target
+
     log_target = backend.np.where(eigenvalues_state != 0.0, log_target, 0.0)
 
     entropy_state = backend.np.sum(eigenvalues_state * log_state)
@@ -791,7 +806,7 @@ def relative_renyi_entropy(
             \\sigma^{1 - \\alpha} \\right) \\right) \\, .
 
     A special case is the limit :math:`\\alpha \\to 1`, in which the Rényi entropy
-    coincides with the :func:`qibo.quantum_info.entropies.relative_entropy`.
+    coincides with the :func:`qibo.quantum_info.entropies.relative_von_neumann_entropy`.
 
     In the limit :math:`\\alpha \\to \\infty`, the function reduces to
     :math:`-2 \\, \\log(\\|\\sqrt{\\rho} \\, \\sqrt{\\sigma}\\|_{1})`,
@@ -943,6 +958,87 @@ def tsallis_entropy(state, alpha: float, base: float = 2, backend=None):
     )
 
 
+def relative_tsallis_entropy(
+    state,
+    target,
+    alpha: Union[float, int],
+    base: float = 2,
+    check_hermitian: bool = False,
+    backend=None,
+):
+    """Calculate the relative Tsallis entropy between two quantum states.
+
+    For :math:`\\alpha \\in [0, \\, 2]` and quantum states :math:`\\rho` and
+    :math:`\\sigma`, the relative Tsallis entropy is defined as
+
+    .. math::
+        \\Delta_{\\alpha}^{\\text{ts}}(\\rho, \\, \\sigma) = \\frac{1 -
+            \\text{tr}\\left(\\rho^{\\alpha} \\, \\sigma^{1 - \\alpha}\\right)}{1 - \\alpha} \\, .
+
+    A special case is the limit :math:`\\alpha \\to 1`, in which the Tsallis entropy
+    coincides with the :func:`qibo.quantum_info.entropies.relative_von_neumann_entropy`.
+
+    Args:
+        state (ndarray): statevector or density matrix :math:`\\rho`.
+        target (ndarray): statevector or density matrix :math:`\\sigma`.
+        alpha (float or int): entropic index :math:`\\alpha \\in [0, \\, 2]`.
+        base (float, optional): the base of the log used when :math:`\\alpha = 1`.
+            Defaults to :math:`2`.
+        check_hermitian (bool, optional): Used when :math:`\\alpha = 1`.
+            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: Relative Tsallis entropy :math:`\\Delta_{\\alpha}^{\\text{ts}}`.
+
+    References:
+        1. S. Abe, *Nonadditive generalization of the quantum Kullback-Leibler
+        divergence for measuring the degree of purification*,
+        `Phys. Rev. A 68, 032302 <https://doi.org/10.1103/PhysRevA.68.032302>`_.
+
+        2. S. Furuichi, K. Yanagi, and K. Kuriyama,
+        *Fundamental properties of Tsallis relative entropy*,
+        `J. Math. Phys., Vol. 45, Issue 12, pp. 4868-4877 (2004)
+        <https://doi.org/10.1063/1.1805729>`_ .
+    """
+    if alpha == 1.0:
+        return relative_von_neumann_entropy(
+            state, target, base=base, check_hermitian=check_hermitian, backend=backend
+        )
+
+    if not isinstance(alpha, (float, int)):
+        raise_error(
+            TypeError,
+            f"``alpha`` must be type float or int, but it is type {type(alpha)}.",
+        )
+
+    if alpha < 0.0 or alpha > 2.0:
+        raise_error(
+            ValueError, f"``alpha`` must be in the interval [0, 2], but it is {alpha}."
+        )
+
+    if alpha < 1.0:
+        alpha = 2 - alpha
+
+    factor = 1 - alpha
+
+    if len(state.shape) == 1:
+        state = backend.np.outer(state, backend.np.conj(state.T))
+
+    if len(target.shape) == 1:
+        target = backend.np.outer(target, backend.np.conj(target.T))
+
+    trace = matrix_power(state, alpha, backend=backend)
+    trace = trace @ matrix_power(target, factor, backend=backend)
+    trace = backend.np.trace(trace)
+
+    return (1 - trace) / factor
+
+
 def entanglement_entropy(
     state,
     bipartition,

tests/test_quantum_info_entropies.py

@@ -12,6 +12,7 @@
     entanglement_entropy,
     mutual_information,
     relative_renyi_entropy,
+    relative_tsallis_entropy,
     relative_von_neumann_entropy,
     renyi_entropy,
     shannon_entropy,
@@ -476,9 +477,10 @@ def test_von_neumann_entropy(backend, base, check_hermitian):
     )
 
 
+@pytest.mark.parametrize("statevector", [False, True])
 @pytest.mark.parametrize("check_hermitian", [False, True])
 @pytest.mark.parametrize("base", [2, 10, np.e, 5])
-def test_relative_entropy(backend, base, check_hermitian):
+def test_relative_von_neumann_entropy(backend, base, check_hermitian, statevector):
     with pytest.raises(TypeError):
         state = np.random.rand(2, 3)
         state = backend.cast(state, dtype=state.dtype)
@@ -513,43 +515,36 @@ def test_relative_entropy(backend, base, check_hermitian):
     nqubits = 2
     dims = 2**nqubits
 
-    state = random_density_matrix(dims, backend=backend)
+    state = (
+        random_statevector(dims, backend=backend)
+        if statevector
+        else random_density_matrix(dims, backend=backend)
+    )
     target = backend.identity_density_matrix(nqubits, normalize=True)
 
+    entropy = relative_von_neumann_entropy(
+        state, target, base=base, check_hermitian=check_hermitian, backend=backend
+    )
+
+    if statevector:
+        state = backend.np.outer(state, backend.np.conj(state.T))
+
+    entropy_target = von_neumann_entropy(
+        state, base=base, check_hermitian=check_hermitian, backend=backend
+    )
+
     backend.assert_allclose(
-        relative_von_neumann_entropy(state, target, base, check_hermitian, backend),
-        np.log(dims) / np.log(base)
-        - von_neumann_entropy(
-            state, base=base, check_hermitian=check_hermitian, backend=backend
-        ),
-        atol=1e-5,
+        entropy, np.log(dims) / np.log(base) - entropy_target, atol=1e-5
     )
 
-    state = backend.cast([1.0, 0.0], dtype=np.float64)
+    state = random_density_matrix(2, seed=8, pure=False, backend=backend)
     target = backend.cast([0.0, 1.0], dtype=np.float64)
 
-    assert relative_von_neumann_entropy(state, target, backend=backend) == 0.0
-
-    # for coverage when GPUs are present
-    if backend.__class__.__name__ in ["CupyBackend", "CuQuantumBackend"]:
-        with pytest.raises(NotImplementedError):
-            state = random_unitary(4, backend=backend)
-            target = random_density_matrix(4, backend=backend)
-            test = relative_von_neumann_entropy(
-                state, target, base=base, check_hermitian=True, backend=backend
-            )
-        with pytest.raises(NotImplementedError):
-            target = random_unitary(4, backend=backend)
-            state = random_density_matrix(4, backend=backend)
-            test = relative_von_neumann_entropy(
-                state, target, base=base, check_hermitian=True, backend=backend
-            )
-    else:
-        state = random_unitary(4, backend=backend)
-        target = random_unitary(4, backend=backend)
-        test = relative_von_neumann_entropy(
-            state, target, base=base, check_hermitian=True, backend=backend
-        )
+    backend.assert_allclose(
+        relative_von_neumann_entropy(state, target, backend=backend),
+        -0.21227801,
+        atol=1e-8,
+    )
 
 
 @pytest.mark.parametrize("check_hermitian", [False, True])
@@ -674,7 +669,9 @@ def test_relative_renyi_entropy(backend, alpha, base, state_flag, target_flag):
                 relative_renyi_entropy(state, target, alpha, base, backend)
         else:
             if alpha == 1.0:
-                log = relative_von_neumann_entropy(state, target, base, backend=backend)
+                log = relative_von_neumann_entropy(
+                    state, target, base=base, backend=backend
+                )
             elif alpha == np.inf:
                 state_outer = (
                     backend.np.outer(state, backend.np.conj(state.T))
@@ -766,6 +763,63 @@ def test_tsallis_entropy(backend, alpha, base):
     )
 
 
+@pytest.mark.parametrize(
+    ["state_flag", "target_flag"], [[True, True], [False, True], [True, False]]
+)
+@pytest.mark.parametrize("check_hermitian", [False, True])
+@pytest.mark.parametrize("base", [2, 10, np.e, 5])
+@pytest.mark.parametrize("alpha", [0, 0.5, 1, 1.9])
+def test_relative_tsallis_entropy(
+    backend, alpha, base, check_hermitian, state_flag, target_flag
+):
+    state = random_statevector(4, backend=backend)
+    target = random_statevector(4, backend=backend)
+
+    with pytest.raises(TypeError):
+        test = relative_tsallis_entropy(state, target, alpha=1j, backend=backend)
+
+    with pytest.raises(ValueError):
+        test = relative_tsallis_entropy(state, target, alpha=3, backend=backend)
+
+    with pytest.raises(ValueError):
+        test = relative_tsallis_entropy(state, target, alpha=-1.0, backend=backend)
+
+    state = (
+        random_statevector(4, seed=10, backend=backend)
+        if state_flag
+        else random_density_matrix(4, seed=10, backend=backend)
+    )
+    target = (
+        random_statevector(4, seed=11, backend=backend)
+        if target_flag
+        else random_density_matrix(4, seed=11, backend=backend)
+    )
+
+    value = relative_tsallis_entropy(
+        state, target, alpha, base, check_hermitian, backend
+    )
+
+    if alpha == 1.0:
+        target_value = relative_von_neumann_entropy(
+            state, target, base=base, check_hermitian=check_hermitian, backend=backend
+        )
+    else:
+        if alpha < 1.0:
+            alpha = 2 - alpha
+
+        if state_flag:
+            state = backend.np.outer(state, backend.np.conj(state.T))
+
+        if target_flag:
+            target = backend.np.outer(target, backend.np.conj(target.T))
+
+        target_value = matrix_power(state, alpha, backend=backend)
+        target_value = target_value @ matrix_power(target, 1 - alpha, backend=backend)
+        target_value = (1 - backend.np.trace(target_value)) / (1 - alpha)
+
+    backend.assert_allclose(value, target_value, atol=1e-10)
+
+
 @pytest.mark.parametrize("check_hermitian", [False, True])
 @pytest.mark.parametrize("base", [2, 10, np.e, 5])
 @pytest.mark.parametrize("bipartition", [[0], [1]])