Merge branch 'master' into plot_circuit_resolution

doc/source/api-reference/qibo.rst

@@ -401,6 +401,11 @@ with the last qubit as the control qubit and the first qubit as a target qubit.
 .. autofunction:: qibo.models.encodings.entangling_layer
 
 
+Greenberger-Horne-Zeilinger (GHZ) state
+"""""""""""""""""""""""""""""""""""""""
+
+.. autofunction:: qibo.models.encodings.ghz_state
+
 .. _error-mitigation:
 
 Error Mitigation
@@ -1826,6 +1831,12 @@ Tsallis entropy
 .. autofunction:: qibo.quantum_info.tsallis_entropy
 
 
+Relative Tsallis entropy
+""""""""""""""""""""""""
+
+.. autofunction:: qibo.quantum_info.relative_tsallis_entropy
+
+
 Entanglement entropy
 """"""""""""""""""""
 

src/qibo/models/encodings.py

@@ -376,6 +376,35 @@ def entangling_layer(
     return circuit
 
 
+def ghz_state(nqubits: int, **kwargs):
+    """Generates an :math:`n`-qubit Greenberger-Horne-Zeilinger (GHZ) state that takes the form
+
+    .. math::
+        \\ket{\\text{GHZ}} = \\frac{\\ket{0}^{\\otimes n} + \\ket{1}^{\\otimes n}}{\\sqrt{2}}
+
+    where :math:`n` is the number of qubits.
+
+    Args:
+        nqubits (int): number of qubits :math:`n >= 2`.
+        kwargs (dict, optional): additional arguments used to initialize a Circuit object.
+            For details, see the documentation of :class:`qibo.models.circuit.Circuit`.
+
+    Returns:
+        :class:`qibo.models.circuit.Circuit`: Circuit that prepares the GHZ state.
+    """
+    if nqubits < 2:
+        raise_error(
+            ValueError,
+            f"nqubits given as {nqubits}. nqubits needs to be >= 2.",
+        )
+
+    circuit = Circuit(nqubits, **kwargs)
+    circuit.add(gates.H(0))
+    circuit.add(gates.CNOT(qubit, qubit + 1) for qubit in range(nqubits - 1))
+
+    return circuit
+
+
 def _generate_rbs_pairs(nqubits: int, architecture: str, **kwargs):
     """Generating list of indexes representing the RBS connections
 

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_models_encodings.py

@@ -11,6 +11,7 @@
 from qibo.models.encodings import (
     comp_basis_encoder,
     entangling_layer,
+    ghz_state,
     phase_encoder,
     unary_encoder,
     unary_encoder_random_gaussian,
@@ -279,3 +280,24 @@ def test_circuit_kwargs(density_matrix):
 
     test = unary_encoder_random_gaussian(4, density_matrix=density_matrix)
     assert test.density_matrix is density_matrix
+
+
+@pytest.mark.parametrize("density_matrix", [False, True])
+@pytest.mark.parametrize("nqubits", [1, 2, 3, 4])
+def test_ghz_circuit(backend, nqubits, density_matrix):
+    if nqubits < 2:
+        with pytest.raises(ValueError):
+            GHZ_circ = ghz_state(nqubits, density_matrix=density_matrix)
+    else:
+        target = np.zeros(2**nqubits, dtype=complex)
+        target[0] = 1 / np.sqrt(2)
+        target[2**nqubits - 1] = 1 / np.sqrt(2)
+        target = backend.cast(target, dtype=target.dtype)
+
+        GHZ_circ = ghz_state(nqubits, density_matrix=density_matrix)
+        state = backend.execute_circuit(GHZ_circ).state()
+
+        if density_matrix:
+            target = backend.np.outer(target, backend.np.conj(target.T))
+
+        backend.assert_allclose(state, target)