Merge branch 'master' into fix_autodiff

src/qibo/backends/abstract.py

@@ -358,7 +358,7 @@ def calculate_matrix_exp(
 
     @abc.abstractmethod
     def calculate_matrix_power(
-        self, matrix, power: Union[float, int]
+        self, matrix, power: Union[float, int], precision_singularity: float = 1e-14
     ):  # pragma: no cover
         """Calculate the (fractional) ``power`` :math:`\\alpha` of ``matrix`` :math:`A`,
         i.e. :math:`A^{\\alpha}`.

src/qibo/backends/numpy.py

@@ -415,8 +415,12 @@ def execute_circuit(self, circuit, initial_state=None, nshots=1000):
             if circuit.measurements or circuit.has_collapse:
                 return self.execute_circuit_repeated(circuit, nshots, initial_state)
             else:
-                raise RuntimeError(
-                    "Attempting to perform noisy simulation with `density_matrix=False` and no Measurement gate in the Circuit. If you wish to retrieve the statistics of the outcomes please include measurements in the circuit, otherwise set `density_matrix=True` to recover the final state."
+                raise_error(
+                    RuntimeError,
+                    "Attempting to perform noisy simulation with `density_matrix=False` "
+                    + "and no Measurement gate in the Circuit. If you wish to retrieve the "
+                    + "statistics of the outcomes please include measurements in the circuit, "
+                    + "otherwise set `density_matrix=True` to recover the final state.",
                 )
 
         if circuit.accelerators:  # pragma: no cover
@@ -495,8 +499,11 @@ def execute_circuit_repeated(self, circuit, nshots, initial_state=None):
             and not circuit.measurements
             and not circuit.density_matrix
         ):
-            raise RuntimeError(
-                "The circuit contains only collapsing measurements (`collapse=True`) but `density_matrix=False`. Please set `density_matrix=True` to retrieve the final state after execution."
+            raise_error(
+                RuntimeError,
+                "The circuit contains only collapsing measurements (`collapse=True`) but "
+                + "`density_matrix=False`. Please set `density_matrix=True` to retrieve "
+                + "the final state after execution.",
             )
 
         results, final_states = [], []
@@ -767,12 +774,26 @@ def calculate_matrix_exp(self, a, matrix, eigenvectors=None, eigenvalues=None):
         ud = self.np.transpose(np.conj(eigenvectors))
         return self.np.matmul(eigenvectors, self.np.matmul(expd, ud))
 
-    def calculate_matrix_power(self, matrix, power: Union[float, int]):
+    def calculate_matrix_power(
+        self,
+        matrix,
+        power: Union[float, int],
+        precision_singularity: float = 1e-14,
+    ):
         if not isinstance(power, (float, int)):
             raise_error(
                 TypeError,
                 f"``power`` must be either float or int, but it is type {type(power)}.",
             )
+
+        if power < 0.0:
+            # negative powers of singular matrices via SVD
+            determinant = self.np.linalg.det(matrix)
+            if abs(determinant) < precision_singularity:
+                return _calculate_negative_power_singular_matrix(
+                    matrix, power, precision_singularity, self.np, self
+                )
+
         return fractional_matrix_power(matrix, power)
 
     def calculate_singular_value_decomposition(self, matrix):
@@ -818,3 +839,15 @@ def _test_regressions(self, name):
                 {5: 18, 4: 5, 7: 4, 1: 2, 6: 1},
                 {4: 8, 2: 6, 5: 5, 1: 3, 3: 3, 6: 2, 7: 2, 0: 1},
             ]
+
+
+def _calculate_negative_power_singular_matrix(
+    matrix, power: Union[float, int], precision_singularity: float, engine, backend
+):
+    """Calculate negative power of singular matrix."""
+    U, S, Vh = backend.calculate_singular_value_decomposition(matrix)
+    # cast needed because of different dtypes in `torch`
+    S = backend.cast(S)
+    S_inv = engine.where(engine.abs(S) < precision_singularity, 0.0, S**power)
+
+    return engine.linalg.inv(Vh) @ backend.np.diag(S_inv) @ engine.linalg.inv(U)

src/qibo/backends/pytorch.py

@@ -1,5 +1,7 @@
 """PyTorch backend."""
 
+from typing import Union
+
 import numpy as np
 
 from qibo import __version__
@@ -226,9 +228,15 @@ def calculate_matrix_exp(self, a, matrix, eigenvectors=None, eigenvalues=None):
         ud = self.np.conj(eigenvectors).T
         return self.np.matmul(eigenvectors, self.np.matmul(expd, ud))
 
-    def calculate_matrix_power(self, matrix, power):
-        copied = self.to_numpy(self.np.copy(matrix))
-        copied = super().calculate_matrix_power(copied, power)
+    def calculate_matrix_power(
+        self,
+        matrix,
+        power: Union[float, int],
+        precision_singularity: float = 1e-14,
+    ):
+        copied = self.cast(matrix, copy=True)
+        copied = self.to_numpy(copied) if power >= 0.0 else copied.detach()
+        copied = super().calculate_matrix_power(copied, power, precision_singularity)
         return self.cast(copied, dtype=copied.dtype)
 
     def _test_regressions(self, name):

src/qibo/backends/tensorflow.py

@@ -1,11 +1,12 @@
 import collections
 import os
+from typing import Union
 
 import numpy as np
 
 from qibo import __version__
 from qibo.backends.npmatrices import NumpyMatrices
-from qibo.backends.numpy import NumpyBackend
+from qibo.backends.numpy import NumpyBackend, _calculate_negative_power_singular_matrix
 from qibo.config import TF_LOG_LEVEL, log, raise_error
 
 
@@ -192,6 +193,28 @@ def calculate_matrix_exp(self, a, matrix, eigenvectors=None, eigenvalues=None):
             return self.tf.linalg.expm(-1j * a * matrix)
         return super().calculate_matrix_exp(a, matrix, eigenvectors, eigenvalues)
 
+    def calculate_matrix_power(
+        self,
+        matrix,
+        power: Union[float, int],
+        precision_singularity: float = 1e-14,
+    ):
+        if not isinstance(power, (float, int)):
+            raise_error(
+                TypeError,
+                f"``power`` must be either float or int, but it is type {type(power)}.",
+            )
+
+        if power < 0.0:
+            # negative powers of singular matrices via SVD
+            determinant = self.tf.linalg.det(matrix)
+            if abs(determinant) < precision_singularity:
+                return _calculate_negative_power_singular_matrix(
+                    matrix, power, precision_singularity, self.tf, self
+                )
+
+        return super().calculate_matrix_power(matrix, power, precision_singularity)
+
     def calculate_singular_value_decomposition(self, matrix):
         # needed to unify order of return
         S, U, V = self.tf.linalg.svd(matrix)

src/qibo/quantum_info/entanglement.py

@@ -148,9 +148,9 @@ def negativity(state, bipartition, backend=None):
 
     reduced = partial_transpose(state, bipartition, backend)
     reduced = backend.np.conj(reduced.T) @ reduced
-    norm = backend.np.trace(matrix_power(reduced, 1 / 2, backend))
+    norm = backend.np.trace(matrix_power(reduced, 1 / 2, backend=backend))
 
-    return float(backend.np.real((norm - 1) / 2))
+    return backend.np.real((norm - 1) / 2)
 
 
 def entanglement_fidelity(

src/qibo/quantum_info/entropies.py

@@ -772,7 +772,7 @@ def renyi_entropy(state, alpha: Union[float, int], base: float = 2, backend=None
             / np.log2(base)
         )
 
-    log = backend.np.log2(backend.np.trace(matrix_power(state, alpha, backend)))
+    log = backend.np.log2(backend.np.trace(matrix_power(state, alpha, backend=backend)))
 
     return (1 / (1 - alpha)) * log / np.log2(base)
 
@@ -871,17 +871,17 @@ def relative_renyi_entropy(
         return relative_von_neumann_entropy(state, target, base, backend=backend)
 
     if alpha == np.inf:
-        new_state = matrix_power(state, 0.5, backend)
-        new_target = matrix_power(target, 0.5, backend)
+        new_state = matrix_power(state, 0.5, backend=backend)
+        new_target = matrix_power(target, 0.5, backend=backend)
 
         log = backend.np.log2(
             backend.calculate_norm_density_matrix(new_state @ new_target, order=1)
         )
 
         return -2 * log / np.log2(base)
 
-    log = matrix_power(state, alpha, backend)
-    log = log @ matrix_power(target, 1 - alpha, backend)
+    log = matrix_power(state, alpha, backend=backend)
+    log = log @ matrix_power(target, 1 - alpha, backend=backend)
     log = backend.np.log2(backend.np.trace(log))
 
     return (1 / (alpha - 1)) * log / np.log2(base)
@@ -939,7 +939,7 @@ def tsallis_entropy(state, alpha: float, base: float = 2, backend=None):
         return von_neumann_entropy(state, base=base, backend=backend)
 
     return (1 / (1 - alpha)) * (
-        backend.np.trace(matrix_power(state, alpha, backend)) - 1
+        backend.np.trace(matrix_power(state, alpha, backend=backend)) - 1
     )
 
 

src/qibo/quantum_info/linalg_operations.py

@@ -277,12 +277,17 @@ def matrix_exponentiation(
     return backend.calculate_matrix_exp(phase, matrix, eigenvectors, eigenvalues)
 
 
-def matrix_power(matrix, power: Union[float, int], backend=None):
+def matrix_power(
+    matrix, power: Union[float, int], precision_singularity: float = 1e-14, backend=None
+):
     """Given a ``matrix`` :math:`A` and power :math:`\\alpha`, calculate :math:`A^{\\alpha}`.
 
     Args:
         matrix (ndarray): matrix whose power to calculate.
         power (float or int): power to raise ``matrix`` to.
+        precision_singularity (float, optional): If determinant of ``matrix`` is smaller than
+            ``precision_singularity``, then matrix is considered to be singular.
+            Used when ``power`` is negative.
         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``.
@@ -292,7 +297,7 @@ def matrix_power(matrix, power: Union[float, int], backend=None):
     """
     backend = _check_backend(backend)
 
-    return backend.calculate_matrix_power(matrix, power)
+    return backend.calculate_matrix_power(matrix, power, precision_singularity)
 
 
 def singular_value_decomposition(matrix, backend=None):
@@ -314,7 +319,8 @@ def singular_value_decomposition(matrix, backend=None):
             :class:`qibo.backends.GlobalBackend`. Defaults to ``None``.
 
     Returns:
-        ndarray, ndarray, ndarray: Singular value decomposition of :math:`A`.
+        ndarray, ndarray, ndarray: Singular value decomposition of :math:`A`, i.e.
+        :math:`U`, :math:`S`, and :math:`V^{\\dagger}`, in that order.
     """
     backend = _check_backend(backend)
 

tests/test_quantum_info_entropies.py

@@ -686,8 +686,8 @@ def test_relative_renyi_entropy(backend, alpha, base, state_flag, target_flag):
                     if target_flag
                     else target
                 )
-                new_state = matrix_power(state_outer, 0.5, backend)
-                new_target = matrix_power(target_outer, 0.5, backend)
+                new_state = matrix_power(state_outer, 0.5, backend=backend)
+                new_target = matrix_power(target_outer, 0.5, backend=backend)
 
                 log = backend.np.log2(
                     backend.calculate_norm_density_matrix(
@@ -703,8 +703,8 @@ def test_relative_renyi_entropy(backend, alpha, base, state_flag, target_flag):
                 if len(target.shape) == 1:
                     target = backend.np.outer(target, backend.np.conj(target))
 
-                log = matrix_power(state, alpha, backend)
-                log = log @ matrix_power(target, 1 - alpha, backend)
+                log = matrix_power(state, alpha, backend=backend)
+                log = log @ matrix_power(target, 1 - alpha, backend=backend)
                 log = backend.np.log2(backend.np.trace(log))
 
                 log = (1 / (alpha - 1)) * log / np.log2(base)
@@ -750,7 +750,7 @@ def test_tsallis_entropy(backend, alpha, base):
         target = von_neumann_entropy(state, base=base, backend=backend)
     else:
         target = (1 / (1 - alpha)) * (
-            backend.np.trace(matrix_power(state, alpha, backend)) - 1
+            backend.np.trace(matrix_power(state, alpha, backend=backend)) - 1
         )
 
     backend.assert_allclose(

tests/test_quantum_info_operations.py

@@ -1,5 +1,6 @@
 import numpy as np
 import pytest
+from scipy.linalg import sqrtm
 
 from qibo import Circuit, gates, matrices
 from qibo.quantum_info.linalg_operations import (
@@ -202,18 +203,30 @@ def test_partial_transpose(backend, p, statevector, batch):
             backend.assert_allclose(transposed, target)
 
 
-@pytest.mark.parametrize("power", [2, 2.0, "2"])
-def test_matrix_power(backend, power):
+@pytest.mark.parametrize("singular", [False, True])
+@pytest.mark.parametrize("power", [-0.5, 0.5, 2, 2.0, "2"])
+def test_matrix_power(backend, power, singular):
     nqubits = 2
     dims = 2**nqubits
 
-    state = random_density_matrix(dims, backend=backend)
+    state = random_density_matrix(dims, pure=singular, backend=backend)
 
     if isinstance(power, str):
         with pytest.raises(TypeError):
-            test = matrix_power(state, power, backend)
+            test = matrix_power(state, power, backend=backend)
+    elif power == -0.5 and singular:
+        # When the singular matrix is a state, this power should be itself
+        backend.assert_allclose(matrix_power(state, power, backend=backend), state)
+    elif abs(power) == 0.5 and not singular:
+        # Should be equal to the (inverse) square root
+        sqrt = sqrtm(backend.to_numpy(state)).astype(complex)
+        if power == -0.5:
+            sqrt = np.linalg.inv(sqrt)
+        sqrt = backend.cast(sqrt)
+
+        backend.assert_allclose(matrix_power(state, power, backend=backend), sqrt)
     else:
-        power = matrix_power(state, power, backend)
+        power = matrix_power(state, power, backend=backend)
 
         backend.assert_allclose(
             float(backend.np.real(backend.np.trace(power))),