Add matrix_power to quantum_info.linalg_operations

doc/source/api-reference/qibo.rst

@@ -1952,6 +1952,12 @@ Matrix exponentiation
 .. autofunction:: qibo.quantum_info.matrix_exponentiation
 
 
+Matrix power
+""""""""""""
+
+.. autofunction:: qibo.quantum_info.matrix_power
+
+
 Quantum Networks
 ^^^^^^^^^^^^^^^^
 

pyproject.toml

@@ -57,7 +57,7 @@ ipython = "^8.10.0"
 qulacs = { version = "^0.6.4", markers = "(sys_platform == 'darwin' and python_version > '3.9') or sys_platform != 'darwin'" }
 seaborn = "^0.13.2"
 ipykernel = "^6.29.4"
-qibojit = { git = "https://github.com/qiboteam/qibojit.git" }
+qibojit = { git = "https://github.com/qiboteam/qibojit.git" , branch = "matrix_power" }
 
 [tool.poetry.group.tests]
 optional = true
@@ -72,7 +72,7 @@ pylint = "3.1.0"
 matplotlib = "^3.7.0"
 tensorflow = { version = "^2.16.1", markers = "sys_platform == 'linux'" }
 torch = "^2.1.1,<2.4"
-qibojit = { git = "https://github.com/qiboteam/qibojit.git" }
+qibojit = { git = "https://github.com/qiboteam/qibojit.git" , branch = "matrix_power" }
 qibotn = { git = "https://github.com/qiboteam/qibotn.git" }
 stim = "^1.12.0"
 qulacs = { version = "^0.6.4", markers = "(sys_platform == 'darwin' and python_version > '3.9') or sys_platform != 'darwin'" }
@@ -92,7 +92,7 @@ optional = true
 [tool.poetry.group.cuda11.dependencies]
 cupy-cuda11x = "^12.0.0"
 cuquantum-python-cu11 = "^23.3.0"
-qibojit = { git = "https://github.com/qiboteam/qibojit.git" }
+qibojit = { git = "https://github.com/qiboteam/qibojit.git" , branch = "matrix_power" }
 qibotn = { git = "https://github.com/qiboteam/qibotn.git" }
 
 
@@ -102,7 +102,7 @@ optional = true
 [tool.poetry.group.cuda12.dependencies]
 cupy-cuda12x = "^12.0.0"
 cuquantum-python-cu12 = "^23.3.0"
-qibojit = { git = "https://github.com/qiboteam/qibojit.git" }
+qibojit = { git = "https://github.com/qiboteam/qibojit.git" , branch = "matrix_power" }
 qibotn = { git = "https://github.com/qiboteam/qibotn.git" }
 
 [tool.poetry.extras]

src/qibo/backends/abstract.py

@@ -1,4 +1,5 @@
 import abc
+from typing import Union
 
 from qibo.config import raise_error
 
@@ -351,6 +352,14 @@ def calculate_matrix_exp(
         """
         raise_error(NotImplementedError)
 
+    @abc.abstractmethod
+    def calculate_matrix_power(
+        self, matrix, power: Union[float, int]
+    ):  # pragma: no cover
+        """Calculates the (fractional) ``power`` :math:`\\alpha` of ``matrix`` :math:`A`,
+        i.e. :math:`A^{\\alpha}`."""
+        raise_error(NotImplementedError)
+
     @abc.abstractmethod
     def calculate_hamiltonian_matrix_product(
         self, matrix1, matrix2

src/qibo/backends/numpy.py

@@ -1,9 +1,10 @@
 import collections
 import math
+from typing import Union
 
 import numpy as np
 from scipy import sparse
-from scipy.linalg import block_diag
+from scipy.linalg import block_diag, fractional_matrix_power
 
 from qibo import __version__
 from qibo.backends import einsum_utils
@@ -768,6 +769,14 @@ 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]):
+        if not isinstance(power, (float, int)):
+            raise_error(
+                TypeError,
+                f"``power`` must be either float or int, but it is type {type(power)}.",
+            )
+        return fractional_matrix_power(matrix, power)
+
     # TODO: remove this method
     def calculate_hamiltonian_matrix_product(self, matrix1, matrix2):
         return matrix1 @ matrix2

src/qibo/backends/pytorch.py

@@ -195,6 +195,11 @@ 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)
+        return self.cast(copied, dtype=copied.dtype)
+
     def _test_regressions(self, name):
         if name == "test_measurementresult_apply_bitflips":
             return [

src/qibo/quantum_info/entropies.py

@@ -3,12 +3,10 @@
 from typing import Union
 
 import numpy as np
-from scipy.linalg import fractional_matrix_power
 
 from qibo.backends import _check_backend
-from qibo.backends.pytorch import PyTorchBackend
 from qibo.config import PRECISION_TOL, raise_error
-from qibo.quantum_info.linalg_operations import partial_trace
+from qibo.quantum_info.linalg_operations import matrix_power, partial_trace
 from qibo.quantum_info.metrics import _check_hermitian_or_not_gpu, purity
 
 
@@ -652,7 +650,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)))
 
     return (1 / (1 - alpha)) * log / np.log2(base)
 
@@ -751,17 +749,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)
+        new_target = matrix_power(target, 0.5, 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)
+    log = log @ matrix_power(target, 1 - alpha, backend)
     log = backend.np.log2(backend.np.trace(log))
 
     return (1 / (alpha - 1)) * log / np.log2(base)
@@ -819,7 +817,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)) - 1
     )
 
 
@@ -881,19 +879,3 @@ def entanglement_entropy(
     )
 
     return entropy_entanglement
-
-
-def _matrix_power(matrix, alpha, backend):
-    """Calculates ``matrix ** alpha`` according to backend."""
-    if backend.__class__.__name__ in [
-        "CupyBackend",
-        "CuQuantumBackend",
-    ]:  # pragma: no cover
-        new_matrix = backend.to_numpy(matrix)
-    else:
-        new_matrix = backend.np.copy(matrix)
-
-    if len(new_matrix.shape) == 1:
-        new_matrix = backend.np.outer(new_matrix, backend.np.conj(new_matrix))
-
-    return backend.cast(fractional_matrix_power(backend.to_numpy(new_matrix), alpha))

src/qibo/quantum_info/linalg_operations.py

@@ -196,3 +196,28 @@ def matrix_exponentiation(
     backend = _check_backend(backend)
 
     return backend.calculate_matrix_exp(phase, matrix, eigenvectors, eigenvalues)
+
+
+def matrix_power(matrix, power: Union[float, int], backend=None):
+    """Given a ``matrix`` :math:`A` and power :math:`\\alpha`, calculates :math:`A^{\\alpha}`.
+
+    Args:
+        matrix (ndarray): matrix whose power to calculate.
+        power (float or int): power to raise ``matrix`` to.
+        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:
+        ndarray: matrix power :math:`A^{\\alpha}`.
+    """
+
+    if not isinstance(power, (float, int)):
+        raise_error(
+            TypeError,
+            f"``power`` must be either float or int, but it is type {type(power)}.",
+        )
+
+    backend = _check_backend(backend)
+
+    return backend.calculate_matrix_power(matrix, power)

tests/test_quantum_info_entropies.py

@@ -1,10 +1,8 @@
 import numpy as np
 import pytest
-from scipy.linalg import sqrtm
 
 from qibo.config import PRECISION_TOL
 from qibo.quantum_info.entropies import (
-    _matrix_power,
     classical_relative_entropy,
     classical_relative_renyi_entropy,
     classical_renyi_entropy,
@@ -17,6 +15,7 @@
     tsallis_entropy,
     von_neumann_entropy,
 )
+from qibo.quantum_info.linalg_operations import matrix_power
 from qibo.quantum_info.random_ensembles import (
     random_density_matrix,
     random_statevector,
@@ -604,8 +603,18 @@ def test_relative_renyi_entropy(backend, alpha, base, state_flag, target_flag):
             if alpha == 1.0:
                 log = relative_von_neumann_entropy(state, target, base, backend=backend)
             elif alpha == np.inf:
-                new_state = _matrix_power(state, 0.5, backend)
-                new_target = _matrix_power(target, 0.5, backend)
+                state_outer = (
+                    backend.np.outer(state, backend.np.conj(state.T))
+                    if state_flag
+                    else state
+                )
+                target_outer = (
+                    backend.np.outer(target, backend.np.conj(target.T))
+                    if target_flag
+                    else target
+                )
+                new_state = matrix_power(state_outer, 0.5, backend)
+                new_target = matrix_power(target_outer, 0.5, backend)
 
                 log = backend.np.log2(
                     backend.calculate_norm_density_matrix(
@@ -621,8 +630,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)
+                log = log @ matrix_power(target, 1 - alpha, backend)
                 log = backend.np.log2(backend.np.trace(log))
 
                 log = (1 / (alpha - 1)) * log / np.log2(base)
@@ -668,7 +677,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)) - 1
         )
 
     backend.assert_allclose(

tests/test_quantum_info_operations.py

@@ -5,8 +5,10 @@
 from qibo.quantum_info.linalg_operations import (
     anticommutator,
     commutator,
+    matrix_power,
     partial_trace,
 )
+from qibo.quantum_info.metrics import purity
 from qibo.quantum_info.random_ensembles import random_density_matrix, random_statevector
 
 
@@ -109,3 +111,16 @@ def test_partial_trace(backend, density_matrix):
     Id = backend.identity_density_matrix(1, normalize=True)
 
     backend.assert_allclose(traced, Id)
+
+
+def test_matrix_power(backend):
+    nqubits, power = 2, 2
+    dims = 2**nqubits
+
+    state = random_density_matrix(dims, backend=backend)
+
+    power = matrix_power(state, power, backend)
+
+    backend.assert_allclose(
+        float(backend.np.real(backend.np.trace(power))), purity(state, backend=backend)
+    )