diff --git a/doc/source/api-reference/qibo.rst b/doc/source/api-reference/qibo.rst index ee62633812..a322f5c4f9 100644 --- a/doc/source/api-reference/qibo.rst +++ b/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 """""""""""""""""""" diff --git a/src/qibo/quantum_info/entropies.py b/src/qibo/quantum_info/entropies.py index e424a09089..48767896d0 100644 --- a/src/qibo/quantum_info/entropies.py +++ b/src/qibo/quantum_info/entropies.py @@ -561,16 +561,23 @@ 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`. @@ -578,6 +585,9 @@ def relative_von_neumann_entropy( 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 `_. + + 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) + `_ . + """ + 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, diff --git a/tests/test_quantum_info_entropies.py b/tests/test_quantum_info_entropies.py index b955608c39..7c289009f9 100644 --- a/tests/test_quantum_info_entropies.py +++ b/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]])