Merge pull request #1684 from qiboteam/metrics

Add `Backend.calculate_matrix_log` and `Backend.calculate_matrix_sqrt` + refactor `Backend.calculate_matrix_exp` + add `quantum_info.matrix_logarithm` and `quantum_info.matrix_sqrt`

Author: renatomello

doc/source/api-reference/qibo.rst

@@ -1814,27 +1814,12 @@ von Neumann entropy
 
 .. autofunction:: qibo.quantum_info.von_neumann_entropy
 
-.. note::
-    ``check_hermitian`` flag allows the user to choose if the function will check if input
-    ``state`` is Hermitian or not. Default option is ``check_hermitian=False``, i.e. the
-    assumption of Hermiticity. This is faster and, more importantly,
-    this function are intended to be used on Hermitian inputs. When ``check_hermitian=True``
-    and ``state`` is non-Hermitian, an error will be raised when using `cupy` backend.
-
 
 Relative von Neumann entropy
 """"""""""""""""""""""""""""
 
 .. autofunction:: qibo.quantum_info.relative_von_neumann_entropy
 
-.. note::
-    ``check_hermitian`` flag allows the user to choose if the function will check if input
-    ``state`` is Hermitian or not. Default option is ``check_hermitian=False``, i.e. the
-    assumption of Hermiticity. This is faster and, more importantly,
-    this function are intended to be used on Hermitian inputs. When ``check_hermitian=True``
-    and either ``state`` or ``target`` is non-Hermitian,
-    an error will be raised when using `cupy` backend.
-
 
 Mutual information
 """"""""""""""""""
@@ -1871,15 +1856,6 @@ Entanglement entropy
 
 .. autofunction:: qibo.quantum_info.entanglement_entropy
 
-.. note::
-    ``check_hermitian`` flag allows the user to choose if the function will check if
-    the reduced density matrix resulting from tracing out ``bipartition`` from input
-    ``state`` is Hermitian or not. Default option is ``check_hermitian=False``, i.e. the
-    assumption of Hermiticity. This is faster and, more importantly,
-    this function are intended to be used on Hermitian inputs. When ``check_hermitian=True``
-    and the reduced density matrix is non-Hermitian, an error will be raised
-    when using `cupy` backend.
-
 
 Metrics
 ^^^^^^^
@@ -1904,14 +1880,6 @@ Trace distance
 
 .. autofunction:: qibo.quantum_info.trace_distance
 
-.. note::
-    ``check_hermitian`` flag allows the user to choose if the function will check if difference
-    between inputs, ``state - target``, is Hermitian or not. Default option is
-    ``check_hermitian=False``, i.e. the assumption of Hermiticity, because it is faster and,
-    more importantly, the functions are intended to be used on Hermitian inputs.
-    When ``check_hermitian=True`` and ``state - target`` is non-Hermitian, an error will be
-    raised when using `cupy` backend.
-
 
 Hilbert-Schmidt inner product
 """""""""""""""""""""""""""""
@@ -2033,12 +2001,24 @@ Matrix exponentiation
 .. autofunction:: qibo.quantum_info.matrix_exponentiation
 
 
+Matrix logarithm
+""""""""""""""""
+
+.. autofunction:: qibo.quantum_info.matrix_logarithm
+
+
 Matrix power
 """"""""""""
 
 .. autofunction:: qibo.quantum_info.matrix_power
 
 
+Matrix square root
+""""""""""""""""""
+
+.. autofunction:: qibo.quantum_info.matrix_sqrt
+
+
 Singular value decomposition
 """"""""""""""""""""""""""""
 

src/qibo/backends/abstract.py

@@ -418,9 +418,27 @@ def calculate_expectation_density_matrix(
 
     @abc.abstractmethod
     def calculate_matrix_exp(
-        self, a, matrix, eigenvectors=None, eigenvalues=None
+        self,
+        matrix,
+        phase: Union[float, int, complex] = 1,
+        eigenvectors=None,
+        eigenvalues=None,
+    ):  # pragma: no cover
+        """Calculate the exponential :math:`e^{\\theta \\, A}` of a matrix :math:`A`
+        and ``phase`` :math:`\\theta`.
+
+        If the eigenvectors and eigenvalues are given the matrix diagonalization is
+        used for exponentiation.
+        """
+        raise_error(NotImplementedError)
+
+    @abc.abstractmethod
+    def calculate_matrix_log(
+        self, matrix, base: Union[float, int] = 2, eigenvectors=None, eigenvalues=None
     ):  # pragma: no cover
-        """Calculate matrix exponential of a matrix.
+        """Calculate the logarithm :math:`\\log_{b}(A)` with a ``base`` :math:`b`
+        of a matrix :math:`A`.
+
         If the eigenvectors and eigenvalues are given the matrix diagonalization is
         used for exponentiation.
         """
@@ -441,6 +459,19 @@ def calculate_matrix_power(
         """
         raise_error(NotImplementedError)
 
+    @abc.abstractmethod
+    def calculate_matrix_sqrt(
+        self, matrix, precision_singularity: float = 1e-14
+    ):  # pragma: no cover
+        """Calculate the square root of ``matrix`` :math:`A`, i.e. :math:`A^{1/2}`.
+
+        .. note::
+            For the ``pytorch`` backend, this method relies on a copy of the original tensor.
+            This may break the gradient flow. For the GPU backends (i.e. ``cupy`` and
+            ``cuquantum``), this method falls back to CPU.
+        """
+        raise_error(NotImplementedError)
+
     @abc.abstractmethod
     def calculate_singular_value_decomposition(self, matrix):  # pragma: no cover
         """Calculate the Singular Value Decomposition of ``matrix``."""

src/qibo/backends/numpy.py

@@ -6,7 +6,7 @@
 
 import numpy as np
 from scipy import sparse
-from scipy.linalg import block_diag, fractional_matrix_power
+from scipy.linalg import block_diag, fractional_matrix_power, logm
 
 from qibo import __version__
 from qibo.backends import einsum_utils
@@ -793,16 +793,39 @@ def calculate_expectation_density_matrix(self, hamiltonian, state, normalize):
             ev /= norm
         return ev
 
-    def calculate_matrix_exp(self, a, matrix, eigenvectors=None, eigenvalues=None):
+    def calculate_matrix_exp(
+        self,
+        matrix,
+        phase: Union[float, int, complex] = 1,
+        eigenvectors=None,
+        eigenvalues=None,
+    ):
         if eigenvectors is None or self.is_sparse(matrix):
             if self.is_sparse(matrix):
                 from scipy.sparse.linalg import expm
             else:
                 from scipy.linalg import expm
-            return expm(-1j * a * matrix)
-        expd = self.np.diag(self.np.exp(-1j * a * eigenvalues))
+
+            return expm(phase * matrix)
+
+        expd = self.np.exp(phase * eigenvalues)
+        ud = self.np.transpose(np.conj(eigenvectors))
+
+        return (eigenvectors * expd) @ ud
+
+    def calculate_matrix_log(self, matrix, base=2, eigenvectors=None, eigenvalues=None):
+        if eigenvectors is None:
+            # to_numpy and cast needed for GPUs
+            matrix_log = logm(self.to_numpy(matrix)) / float(np.log(base))
+            matrix_log = self.cast(matrix_log, dtype=matrix_log.dtype)
+
+            return matrix_log
+
+        # log = self.np.diag(self.np.log(eigenvalues) / float(np.log(base)))
+        log = self.np.log(eigenvalues) / float(np.log(base))
         ud = self.np.transpose(np.conj(eigenvectors))
-        return self.np.matmul(eigenvectors, self.np.matmul(expd, ud))
+
+        return (eigenvectors * log) @ ud
 
     def calculate_matrix_power(
         self,
@@ -826,6 +849,9 @@ def calculate_matrix_power(
 
         return fractional_matrix_power(matrix, power)
 
+    def calculate_matrix_sqrt(self, matrix):
+        return self.calculate_matrix_power(matrix, power=0.5)
+
     def calculate_singular_value_decomposition(self, matrix):
         return self.np.linalg.svd(matrix)
 

src/qibo/callbacks.py

@@ -101,13 +101,11 @@ def __init__(
         partition: Optional[List[int]] = None,
         compute_spectrum: bool = False,
         base: float = 2,
-        check_hermitian: bool = False,
     ):
         super().__init__()
         self.partition = partition
         self.compute_spectrum = compute_spectrum
         self.base = base
-        self.check_hermitian = check_hermitian
         self.spectrum = list()
 
     @Callback.nqubits.setter
@@ -132,7 +130,6 @@ def apply(self, backend, state):
             state,
             bipartition=self.partition,
             base=self.base,
-            check_hermitian=self.check_hermitian,
             return_spectrum=True,
             backend=backend,
         )

src/qibo/hamiltonians/hamiltonians.py

@@ -92,7 +92,11 @@ def exp(self, a):
         if self._exp.get("a") != a:
             self._exp["a"] = a
             self._exp["result"] = matrix_exponentiation(
-                a, self.matrix, self._eigenvectors, self._eigenvalues, self.backend
+                self.matrix,
+                -1j * a,
+                self._eigenvectors,
+                self._eigenvalues,
+                self.backend,
             )
         return self._exp.get("result")
 

src/qibo/hamiltonians/terms.py

@@ -62,7 +62,7 @@ def gate(self):
 
     def exp(self, x):
         """Matrix exponentiation of the term."""
-        return self.backend.calculate_matrix_exp(x, self.matrix)
+        return self.backend.calculate_matrix_exp(self.matrix, phase=-1j * x)
 
     def expgate(self, x):
         """:class:`qibo.gates.gates.Unitary` gate implementing the action of exp(term) on states."""

src/qibo/quantum_info/entanglement.py

@@ -153,9 +153,7 @@ def negativity(state, bipartition, backend=None):
     return backend.np.real((norm - 1) / 2)
 
 
-def entanglement_fidelity(
-    channel, nqubits: int, state=None, check_hermitian: bool = False, backend=None
-):
+def entanglement_fidelity(channel, nqubits: int, state=None, backend=None):
     """Entanglement fidelity :math:`F_{\\mathcal{E}}` of a ``channel`` :math:`\\mathcal{E}`
     on ``state`` :math:`\\rho` is given by
 
@@ -175,9 +173,6 @@ def entanglement_fidelity(
             by ``channel``. If ``None``, defaults to the maximally entangled state
             :math:`\\frac{1}{2^{n}} \\, \\sum_{k} \\, \\ket{k}\\ket{k}`, where
             :math:`n` is ``nqubits``. Defaults to ``None``.
-        check_hermitian (bool, optional): if ``True``, checks if the final state
-            :math:`\\rho_{f}` is Hermitian. If ``False``, it assumes it is Hermitian.
-            Defaults to ``False``.
         backend (:class:`qibo.backends.abstract.Backend`, optional): backend to be used
             in the execution. If ``None``, it uses the current backend.
             Defaults to ``None``.
@@ -205,12 +200,6 @@ def entanglement_fidelity(
             f"state must have dims either (k,) or (k,k), but have dims {state.shape}.",
         )
 
-    if not isinstance(check_hermitian, bool):
-        raise_error(
-            TypeError,
-            f"check_hermitian must be type bool, but it is type {type(check_hermitian)}.",
-        )
-
     backend = _check_backend(backend)
 
     if state is None:
@@ -223,9 +212,7 @@ def entanglement_fidelity(
 
     state_final = backend.apply_channel_density_matrix(channel, state, nqubits)
 
-    entang_fidelity = fidelity(
-        state_final, state, check_hermitian=check_hermitian, backend=backend
-    )
+    entang_fidelity = fidelity(state_final, state, backend=backend)
 
     return entang_fidelity