diff --git a/doc/source/api-reference/qibo.rst b/doc/source/api-reference/qibo.rst index 5f80e763c1..0018faa93a 100644 --- a/doc/source/api-reference/qibo.rst +++ b/doc/source/api-reference/qibo.rst @@ -1271,7 +1271,7 @@ Callbacks Callbacks provide a way to calculate quantities on the state vector as it propagates through the circuit. Example of such quantity is the entanglement entropy, which is currently the only callback implemented in -:class:`qibo.callbacks.EntanglementEntropy`. +:class:`qibo.callbacks.Entanglemententropy`. The user can create custom callbacks by inheriting the :class:`qibo.callbacks.Callback` class. The point each callback is calculated inside the circuit is defined by adding a :class:`qibo.gates.CallbackGate`. @@ -1284,7 +1284,7 @@ This can be added similarly to a standard gate and does not affect the state vec Entanglement entropy ^^^^^^^^^^^^^^^^^^^^ -.. autoclass:: qibo.callbacks.EntanglementEntropy +.. autoclass:: qibo.callbacks.Entanglemententropy :members: :member-order: bysource @@ -1537,18 +1537,24 @@ Shannon entropy .. autofunction:: qibo.quantum_info.shannon_entropy -Classical Relative Entropy +Classical relative entropy """""""""""""""""""""""""" .. autofunction:: qibo.quantum_info.classical_relative_entropy -Classical Rényi Entropy +Classical Rényi entropy """"""""""""""""""""""" .. autofunction:: qibo.quantum_info.classical_renyi_entropy +Classical Rényi relative entropy +"""""""""""""""""""""""""""""""" + +.. autofunction:: qibo.quantum_info.classical_relative_renyi_entropy + + Entropy """"""" @@ -1562,7 +1568,7 @@ Entropy and ``state`` is non-Hermitian, an error will be raised when using `cupy` backend. -Relative Entropy +Relative entropy """""""""""""""" .. autofunction:: qibo.quantum_info.relative_entropy @@ -1576,13 +1582,19 @@ Relative Entropy an error will be raised when using `cupy` backend. -Rényi Entropy +Rényi entropy """"""""""""" .. autofunction:: qibo.quantum_info.renyi_entropy -Entanglement Entropy +Rényi relative entropy +"""""""""""""""""""""" + +.. autofunction:: qibo.quantum_info.renyi_relative_entropy + + +Entanglement entropy """""""""""""""""""" .. autofunction:: qibo.quantum_info.entanglement_entropy diff --git a/src/qibo/quantum_info/entropies.py b/src/qibo/quantum_info/entropies.py index 8c674b37a8..606be7d221 100644 --- a/src/qibo/quantum_info/entropies.py +++ b/src/qibo/quantum_info/entropies.py @@ -170,7 +170,7 @@ def classical_renyi_entropy( Another special case is the limit :math:`\\alpha \\to 0`, where the function is reduced to :math:`\\log\\left(|\\mathbf{p}|\\right)`, with :math:`|\\mathbf{p}|` being the support of :math:`\\mathbf{p}`. - This is knows as the `Hartley entropy ` + This is knows as the `Hartley entropy `_ (also known as *Hartley function* or *max-entropy*). In the limit :math:`\\alpha \\to \\infty`, the function reduces to @@ -238,6 +238,106 @@ def classical_renyi_entropy( return renyi_ent +def classical_relative_renyi_entropy( + prob_dist_p, prob_dist_q, alpha: Union[float, int], base: float = 2, backend=None +): + """Calculates the classical relative Rényi entropy between two discrete probability distributions. + + This function is also known as + `Rényi divergence `_. + + For :math:`\\alpha \\in (0, \\, 1) \\cup (1, \\, \\infty)` and probability distributions + :math:`\\mathbf{p}` and :math:`\\mathbf{q}`, the classical relative Rényi entropy is defined as + + .. math:: + H_{\\alpha}(\\mathbf{p} \\, \\| \\, \\mathbf{q}) = \\frac{1}{\\alpha - 1} \\, + \\log\\left( \\sum_{x} \\, \\frac{\\mathbf{p}^{\\alpha}(x)} + {\\mathbf{q}^{\\alpha - 1}(x)} \\right) \\, . + + A special case is the limit :math:`\\alpha \\to 1`, in which the classical Rényi divergence + coincides with the :func:`qibo.quantum_info.entropies.classical_relative_entropy`. + + Another special case is the limit :math:`\\alpha \\to 1/2`, where the function is + reduced to :math:`-2 \\log\\left(\\sum_{x} \\, \\sqrt{\\mathbf{p}(x) \\, \\mathbf{q}(x)} \\right)`. + The sum inside the :math:`\\log` is known as the + `Bhattacharyya coefficient `_. + + In the limit :math:`\\alpha \\to \\infty`, the function reduces to + :math:`\\log(\\max_{x}(\\mathbf{p}(x) \\, \\mathbf{q}(x))`. + + Args: + prob_dist_p (ndarray or list): discrete probability distribution :math:`p`. + prob_dist_q (ndarray or list): discrete probability distribution :math:`q`. + alpha (float or int): order of the Rényi entropy. + 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: Classical relative Rényi entropy :math:`H_{\\alpha}(\\mathbf{p} \\, \\| \\, \\mathbf{q})`. + """ + if backend is None: # pragma: no cover + backend = GlobalBackend() + + if isinstance(prob_dist_p, list): + # np.float64 is necessary instead of native float because of tensorflow + prob_dist_p = backend.cast(prob_dist_p, dtype=np.float64) + if isinstance(prob_dist_q, list): + # np.float64 is necessary instead of native float because of tensorflow + prob_dist_q = backend.cast(prob_dist_q, dtype=np.float64) + + if (len(prob_dist_p.shape) != 1) or (len(prob_dist_q.shape) != 1): + raise_error( + TypeError, + "Probability arrays must have dims (k,) but have " + + f"dims {prob_dist_p.shape} and {prob_dist_q.shape}.", + ) + + if (len(prob_dist_p) == 0) or (len(prob_dist_q) == 0): + raise_error(TypeError, "At least one of the arrays is empty.") + + if not isinstance(alpha, (float, int)): + raise_error( + TypeError, f"alpha must be type float, but it is type {type(alpha)}." + ) + + if alpha < 0.0: + raise_error(ValueError, "alpha must a non-negative float.") + + if base <= 0: + raise_error(ValueError, "log base must be non-negative.") + + if (any(prob_dist_p < 0) or any(prob_dist_p > 1.0)) or ( + any(prob_dist_q < 0) or any(prob_dist_q > 1.0) + ): + raise_error( + ValueError, + "All elements of the probability array must be between 0. and 1..", + ) + if np.abs(np.sum(prob_dist_p) - 1.0) > PRECISION_TOL: + raise_error(ValueError, "First probability array must sum to 1.") + + if np.abs(np.sum(prob_dist_q) - 1.0) > PRECISION_TOL: + raise_error(ValueError, "Second probability array must sum to 1.") + + if alpha == 0.5: + return -2 * np.log2(np.sum(np.sqrt(prob_dist_p * prob_dist_q))) / np.log2(base) + + if alpha == 1.0: + return classical_relative_entropy( + prob_dist_p, prob_dist_q, base=base, backend=backend + ) + + if alpha == np.inf: + return np.log2(max(prob_dist_p / prob_dist_q)) / np.log2(base) + + prob_p = prob_dist_p**alpha + prob_q = prob_dist_q ** (1 - alpha) + + return (1 / (alpha - 1)) * np.log2(np.sum(prob_p * prob_q)) / np.log2(base) + + def entropy( state, base: float = 2, @@ -503,6 +603,101 @@ def renyi_entropy(state, alpha: Union[float, int], base: float = 2, backend=None return (1 / (1 - alpha)) * log / np.log2(base) +def relative_renyi_entropy( + state, target, alpha: Union[float, int], base: float = 2, backend=None +): + if backend is None: # pragma: no cover + backend = GlobalBackend() + + if ( + (len(state.shape) >= 3) + or (len(state) == 0) + or (len(state.shape) == 2 and state.shape[0] != state.shape[1]) + ): + raise_error( + TypeError, + f"state must have dims either (k,) or (k,k), but have dims {state.shape}.", + ) + + if ( + (len(target.shape) >= 3) + or (len(target) == 0) + or (len(target.shape) == 2 and target.shape[0] != target.shape[1]) + ): + raise_error( + TypeError, + f"target must have dims either (k,) or (k,k), but have dims {target.shape}.", + ) + + if not isinstance(alpha, (float, int)): + raise_error( + TypeError, f"alpha must be type float, but it is type {type(alpha)}." + ) + + if alpha < 0.0: + raise_error(ValueError, "alpha must a non-negative float.") + + if base <= 0.0: + raise_error(ValueError, "log base must be non-negative.") + + if ( + abs(purity(state) - 1.0) < PRECISION_TOL + and abs(purity(target) - 1) < PRECISION_TOL + ): + return 0.0 + + if alpha == 1.0: + return relative_entropy(state, target, base, backend=backend) + + if alpha == np.inf: + if backend.__class__.__name__ in ["CupyBackend", "CuQuantumBackend"]: + eigenvalues_state, eigenvectors_state = np.linalg.eigh(state) + new_state = np.zeros_like(state, dtype=complex) + new_state = backend.cast(new_state, dtype=new_state.dtype) + for eigenvalue, eigenstate in zip( + eigenvalues_state, np.transpose(eigenvectors_state) + ): + new_state += np.sqrt(eigenvalue) * np.outer( + eigenstate, np.conj(eigenstate) + ) + + eigenvalues_target, eigenvectors_target = np.linalg.eigh(target) + new_target = np.zeros_like(target, dtype=complex) + new_target = backend.cast(new_target, dtype=new_target.dtype) + for eigenvalue, eigenstate in zip( + eigenvalues_target, np.transpose(eigenvectors_target) + ): + new_target += np.sqrt(eigenstate) * np.outer( + eigenstate, np.conj(eigenstate) + ) + else: + from scipy.linalg import sqrtm # pylint: disable=C0415 + + # astype method needed because of tensorflow + new_state = sqrtm(state).astype("complex128") + new_target = sqrtm(target).astype("complex128") + new_state = backend.cast(new_state, dtype=new_state.dtype) + new_target = backend.cast(new_target, dtype=new_target.dtype) + + log = np.log2( + backend.calculate_norm_density_matrix(new_state @ new_target, order=1) + ) + + return -2 * log / np.log2(base) + + if len(state.shape) == 1: + state = np.outer(state, np.conj(state)) + + if len(target.shape) == 1: + target = np.outer(target, np.conj(target)) + + log = np.linalg.matrix_power(state, alpha) + log = log @ np.linalg.matrix_power(target, 1 - alpha) + log = np.log2(np.trace(log)) + + return (1 / (alpha - 1)) * log / np.log2(base) + + def entanglement_entropy( state, bipartition, diff --git a/tests/test_quantum_info_entropies.py b/tests/test_quantum_info_entropies.py index 702ce8db2e..c7e7ebdf9b 100644 --- a/tests/test_quantum_info_entropies.py +++ b/tests/test_quantum_info_entropies.py @@ -1,13 +1,16 @@ import numpy as np import pytest +from scipy.linalg import sqrtm from qibo.config import PRECISION_TOL from qibo.quantum_info.entropies import ( classical_relative_entropy, + classical_relative_renyi_entropy, classical_renyi_entropy, entanglement_entropy, entropy, relative_entropy, + relative_renyi_entropy, renyi_entropy, shannon_entropy, ) @@ -115,6 +118,103 @@ def test_classical_relative_entropy(backend, base, kind): backend.assert_allclose(divergence, target, atol=1e-5) +@pytest.mark.parametrize("kind", [None, list]) +@pytest.mark.parametrize("base", [2, 10, np.e, 5]) +@pytest.mark.parametrize("alpha", [0, 1 / 2, 1, 2, 3, np.inf]) +def test_classical_relative_renyi_entropy(backend, alpha, base, kind): + with pytest.raises(TypeError): + prob = np.random.rand(1, 2) + prob_q = np.random.rand(1, 5) + prob = backend.cast(prob, dtype=prob.dtype) + prob_q = backend.cast(prob_q, dtype=prob_q.dtype) + test = classical_relative_renyi_entropy( + prob, prob_q, alpha, base, backend=backend + ) + with pytest.raises(TypeError): + prob = np.random.rand(1, 2)[0] + prob_q = np.array([]) + prob = backend.cast(prob, dtype=prob.dtype) + prob_q = backend.cast(prob_q, dtype=prob_q.dtype) + test = classical_relative_renyi_entropy( + prob, prob_q, alpha, base, backend=backend + ) + with pytest.raises(ValueError): + prob = np.array([-1, 2.0]) + prob_q = np.random.rand(1, 5)[0] + prob = backend.cast(prob, dtype=prob.dtype) + prob_q = backend.cast(prob_q, dtype=prob_q.dtype) + test = classical_relative_renyi_entropy( + prob, prob_q, alpha, base, backend=backend + ) + with pytest.raises(ValueError): + prob = np.random.rand(1, 2)[0] + prob_q = np.array([1.0, 0.0]) + prob = backend.cast(prob, dtype=prob.dtype) + prob_q = backend.cast(prob_q, dtype=prob_q.dtype) + test = classical_relative_renyi_entropy( + prob, prob_q, alpha, base, backend=backend + ) + with pytest.raises(ValueError): + prob = np.array([1.0, 0.0]) + prob_q = np.random.rand(1, 2)[0] + prob = backend.cast(prob, dtype=prob.dtype) + prob_q = backend.cast(prob_q, dtype=prob_q.dtype) + test = classical_relative_renyi_entropy( + prob, prob_q, alpha, base, backend=backend + ) + with pytest.raises(ValueError): + prob = np.array([1.0, 0.0]) + prob_q = np.array([0.0, 1.0]) + prob = backend.cast(prob, dtype=prob.dtype) + prob_q = backend.cast(prob_q, dtype=prob_q.dtype) + test = classical_relative_renyi_entropy( + prob, prob_q, alpha, base=-2, backend=backend + ) + with pytest.raises(TypeError): + prob = np.array([1.0, 0.0]) + prob_q = np.array([0.0, 1.0]) + prob = backend.cast(prob, dtype=prob.dtype) + prob_q = backend.cast(prob_q, dtype=prob_q.dtype) + test = classical_relative_renyi_entropy( + prob, prob_q, alpha="1", base=base, backend=backend + ) + with pytest.raises(ValueError): + prob = np.array([1.0, 0.0]) + prob_q = np.array([0.0, 1.0]) + prob = backend.cast(prob, dtype=prob.dtype) + prob_q = backend.cast(prob_q, dtype=prob_q.dtype) + test = classical_relative_renyi_entropy( + prob, prob_q, alpha=-2, base=base, backend=backend + ) + + prob_p = np.random.rand(10) + prob_q = np.random.rand(10) + prob_p /= np.sum(prob_p) + prob_q /= np.sum(prob_q) + + if alpha == 0.5: + target = -2 * np.log2(np.sum(np.sqrt(prob_p * prob_q))) / np.log2(base) + elif alpha == 1.0: + target = classical_relative_entropy(prob_p, prob_q, base=base, backend=backend) + elif alpha == np.inf: + target = np.log2(max(prob_p / prob_q)) / np.log2(base) + else: + target = ( + (1 / (alpha - 1)) + * np.log2(np.sum(prob_p**alpha * prob_q ** (1 - alpha))) + / np.log2(base) + ) + + if kind is not None: + prob_p, prob_q = kind(prob_p), kind(prob_q) + + divergence = classical_relative_renyi_entropy( + prob_p, prob_q, alpha=alpha, base=base, backend=backend + ) + + backend.assert_allclose(divergence, target, atol=1e-5) + + @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]) @@ -343,6 +443,104 @@ def test_renyi_entropy(backend, alpha, base): ) +@pytest.mark.parametrize("base", [2, 10, np.e, 5]) +@pytest.mark.parametrize("alpha", [0, 1, 2, 3, np.inf]) +def test_relative_renyi_entropy(backend, alpha, base): + with pytest.raises(TypeError): + state = np.random.rand(2, 3) + state = backend.cast(state, dtype=state.dtype) + target = random_density_matrix(4, backend=backend) + test = relative_renyi_entropy( + state, target, alpha=alpha, base=base, backend=backend + ) + with pytest.raises(TypeError): + target = np.random.rand(2, 3) + target = backend.cast(target, dtype=target.dtype) + state = random_density_matrix(4, backend=backend) + test = relative_renyi_entropy( + state, target, alpha=alpha, base=base, backend=backend + ) + with pytest.raises(TypeError): + state = random_statevector(4, backend=backend) + target = random_statevector(4, backend=backend) + test = relative_renyi_entropy( + state, target, alpha="2", base=base, backend=backend + ) + with pytest.raises(ValueError): + state = random_statevector(4, backend=backend) + target = random_statevector(4, backend=backend) + test = relative_renyi_entropy( + state, target, alpha=-1, base=base, backend=backend + ) + with pytest.raises(ValueError): + state = random_statevector(4, backend=backend) + target = random_statevector(4, backend=backend) + test = relative_renyi_entropy( + state, target, alpha=alpha, base=0, backend=backend + ) + + state = random_density_matrix(4, backend=backend) + target = random_density_matrix(4, backend=backend) + + if alpha == 1.0: + log = relative_entropy(state, target, base, backend=backend) + elif alpha == np.inf: + if backend.__class__.__name__ in ["CupyBackend", "CuQuantumBackend"]: + eigenvalues_state, eigenvectors_state = np.linalg.eigh(state) + new_state = np.zeros_like(state, dtype=complex) + new_state = backend.cast(new_state, dtype=new_state.dtype) + for eigenvalue, eigenstate in zip( + eigenvalues_state, np.transpose(eigenvectors_state) + ): + new_state += np.sqrt(eigenvalue) * np.outer( + eigenstate, np.conj(eigenstate) + ) + + eigenvalues_target, eigenvectors_target = np.linalg.eigh(target) + new_target = np.zeros_like(target, dtype=complex) + new_target = backend.cast(new_target, dtype=new_target.dtype) + for eigenvalue, eigenstate in zip( + eigenvalues_target, np.transpose(eigenvectors_target) + ): + new_target += np.sqrt(eigenstate) * np.outer( + eigenstate, np.conj(eigenstate) + ) + else: + new_state, new_target = sqrtm(state).astype("complex128"), sqrtm( + target + ).astype("complex128") + new_state = backend.cast(new_state, dtype=new_state.dtype) + new_target = backend.cast(new_target, dtype=new_target.dtype) + + log = np.log2( + backend.calculate_norm_density_matrix(new_state @ new_target, order=1) + ) + + log = -2 * log / np.log2(base) + + else: + log = np.linalg.matrix_power(state, alpha) + log = log @ np.linalg.matrix_power(target, 1 - alpha) + log = np.log2(np.trace(log)) + + log = (1 / (alpha - 1)) * log / np.log2(base) + + backend.assert_allclose( + relative_renyi_entropy(state, target, alpha=alpha, base=base, backend=backend), + log, + atol=1e-5, + ) + + # test pure states + state = random_density_matrix(4, pure=True, backend=backend) + target = random_density_matrix(4, pure=True, backend=backend) + backend.assert_allclose( + relative_renyi_entropy(state, target, alpha=alpha, base=base, backend=backend), + 0.0, + atol=1e-8, + ) + + @pytest.mark.parametrize("check_hermitian", [False, True]) @pytest.mark.parametrize("base", [2, 10, np.e, 5]) @pytest.mark.parametrize("bipartition", [[0], [1]])