Add calculate_singular_value_decompositon as a backend method

doc/source/api-reference/qibo.rst

@@ -1981,6 +1981,12 @@ Matrix power
 .. autofunction:: qibo.quantum_info.matrix_power
 
 
+Singular value decomposition
+""""""""""""""""""""""""""""
+
+.. autofunction:: qibo.quantum_info.singular_value_decomposition
+
+
 Quantum Networks
 ^^^^^^^^^^^^^^^^
 

src/qibo/backends/abstract.py

@@ -371,6 +371,11 @@ def calculate_matrix_power(
         """
         raise_error(NotImplementedError)
 
+    @abc.abstractmethod
+    def calculate_singular_value_decomposition(self, matrix):  # pragma: no cover
+        """Calculate the Singular Value Decomposition of ``matrix``."""
+        raise_error(NotImplementedError)
+
     @abc.abstractmethod
     def calculate_hamiltonian_matrix_product(
         self, matrix1, matrix2

src/qibo/backends/numpy.py

@@ -777,6 +777,9 @@ def calculate_matrix_power(self, matrix, power: Union[float, int]):
             )
         return fractional_matrix_power(matrix, power)
 
+    def calculate_singular_value_decomposition(self, matrix):
+        return self.np.linalg.svd(matrix)
+
     # TODO: remove this method
     def calculate_hamiltonian_matrix_product(self, matrix1, matrix2):
         return matrix1 @ matrix2

src/qibo/backends/pytorch.py

@@ -16,7 +16,7 @@ class TorchMatrices(NumpyMatrices):
     """
 
     def __init__(self, dtype, requires_grad):
-        import torch  # pylint: disable=import-outside-toplevel
+        import torch  # pylint: disable=import-outside-toplevel  # type: ignore
 
         super().__init__(dtype)
         self.np = torch
@@ -38,7 +38,7 @@ def Unitary(self, u):
 class PyTorchBackend(NumpyBackend):
     def __init__(self):
         super().__init__()
-        import torch  # pylint: disable=import-outside-toplevel
+        import torch  # pylint: disable=import-outside-toplevel  # type: ignore
 
         # Global variable to enable or disable gradient calculation
         self.gradients = True

src/qibo/backends/tensorflow.py

@@ -15,7 +15,7 @@ class TensorflowMatrices(NumpyMatrices):
     def __init__(self, dtype):
         super().__init__(dtype)
         import tensorflow as tf  # pylint: disable=import-error
-        import tensorflow.experimental.numpy as tnp  # pylint: disable=import-error
+        import tensorflow.experimental.numpy as tnp  # pylint: disable=import-error  # type: ignore
 
         self.tf = tf
         self.np = tnp
@@ -35,7 +35,7 @@ def __init__(self):
         os.environ["TF_CPP_MIN_LOG_LEVEL"] = str(TF_LOG_LEVEL)
 
         import tensorflow as tf  # pylint: disable=import-error
-        import tensorflow.experimental.numpy as tnp  # pylint: disable=import-error
+        import tensorflow.experimental.numpy as tnp  # pylint: disable=import-error  # type: ignore
 
         if TF_LOG_LEVEL >= 2:
             tf.get_logger().setLevel("ERROR")
@@ -192,6 +192,11 @@ 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_singular_value_decomposition(self, matrix):
+        # needed to unify order of return
+        S, U, V = self.tf.linalg.svd(matrix)
+        return U, S, self.np.conj(self.np.transpose(V))
+
     def calculate_hamiltonian_matrix_product(self, matrix1, matrix2):
         if self.is_sparse(matrix1) or self.is_sparse(matrix2):
             raise_error(

src/qibo/quantum_info/linalg_operations.py

@@ -293,3 +293,29 @@ def matrix_power(matrix, power: Union[float, int], backend=None):
     backend = _check_backend(backend)
 
     return backend.calculate_matrix_power(matrix, power)
+
+
+def singular_value_decomposition(matrix, backend=None):
+    """Calculate the Singular Value Decomposition (SVD) of ``matrix``.
+
+    Given an :math:`M \\times N` complex matrix :math:`A`, its SVD is given by
+
+    .. math:
+        A = U \\, S \\, V^{\\dagger} \\, ,
+
+    where :math:`U` and :math:`V` are, respectively, an :math:`M \\times M`
+    and an :math:`N \\times N` complex unitary matrices, and :math:`S` is an
+    :math:`M \\times N` diagonal matrix with the singular values of :math:`A`.
+
+    Args:
+        matrix (ndarray): matrix whose SVD to calculate.
+        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, ndarray, ndarray: Singular value decomposition of :math:`A`.
+    """
+    backend = _check_backend(backend)
+
+    return backend.calculate_singular_value_decomposition(matrix)

src/qibo/quantum_info/superoperator_transformations.py

@@ -11,6 +11,7 @@
 from qibo.gates.abstract import Gate
 from qibo.gates.gates import Unitary
 from qibo.gates.special import FusedGate
+from qibo.quantum_info.linalg_operations import singular_value_decomposition
 
 
 def vectorization(state, order: str = "row", backend=None):
@@ -483,10 +484,12 @@ def choi_to_kraus(
         warnings.warn("Input choi_super_op is a non-completely positive map.")
 
         # using singular value decomposition because choi_super_op is non-CP
-        U, coefficients, V = np.linalg.svd(backend.to_numpy(choi_super_op))
-        U = np.transpose(U)
-        coefficients = np.sqrt(coefficients)
-        V = np.conj(V)
+        U, coefficients, V = singular_value_decomposition(
+            choi_super_op, backend=backend
+        )
+        U = U.T
+        coefficients = backend.np.sqrt(coefficients)
+        V = backend.np.conj(V)
 
         kraus_left, kraus_right = [], []
         for coeff, eigenvector_left, eigenvector_right in zip(coefficients, U, V):

tests/test_quantum_info_operations.py

@@ -8,6 +8,7 @@
     matrix_power,
     partial_trace,
     partial_transpose,
+    singular_value_decomposition,
 )
 from qibo.quantum_info.metrics import purity
 from qibo.quantum_info.random_ensembles import random_density_matrix, random_statevector
@@ -218,3 +219,34 @@ def test_matrix_power(backend, power):
             float(backend.np.real(backend.np.trace(power))),
             purity(state, backend=backend),
         )
+
+
+def test_singular_value_decomposition(backend):
+    zero = np.array([1, 0], dtype=complex)
+    one = np.array([0, 1], dtype=complex)
+    plus = (zero + one) / np.sqrt(2)
+    minus = (zero - one) / np.sqrt(2)
+    plus = backend.cast(plus, dtype=plus.dtype)
+    minus = backend.cast(minus, dtype=minus.dtype)
+    base = [plus, minus]
+
+    coeffs = np.random.rand(4)
+    coeffs /= np.sum(coeffs)
+    coeffs = backend.cast(coeffs, dtype=coeffs.dtype)
+
+    state = np.zeros((4, 4), dtype=complex)
+    state = backend.cast(state, dtype=state.dtype)
+    for k, coeff in enumerate(coeffs):
+        bitstring = f"{k:0{2}b}"
+        a, b = int(bitstring[0]), int(bitstring[1])
+        ket = backend.np.kron(base[a], base[b])
+        state = state + coeff * backend.np.outer(ket, ket.T)
+
+    _, S, _ = singular_value_decomposition(state, backend=backend)
+
+    S_sorted = backend.np.sort(S)
+    coeffs_sorted = backend.np.sort(coeffs)
+    if backend.name == "pytorch":
+        S_sorted, coeffs_sorted = S_sorted[0], coeffs_sorted[0]
+
+    backend.assert_allclose(S_sorted, coeffs_sorted)

tests/test_quantum_info_superoperator_transformations.py

@@ -1,5 +1,5 @@
 import numpy as np
-import pytest
+import pytest  # type: ignore
 
 from qibo import matrices
 from qibo.config import PRECISION_TOL
@@ -377,14 +377,22 @@ def test_choi_to_pauli(backend, normalize, order, pauli_order, test_superop):
     backend.assert_allclose(test_pauli, pauli_op, atol=PRECISION_TOL)
 
 
+@pytest.mark.parametrize("test_non_CP", [test_non_CP])
 @pytest.mark.parametrize("test_kraus_right", [test_kraus_right])
 @pytest.mark.parametrize("test_kraus_left", [test_kraus_left])
 @pytest.mark.parametrize("test_a1", [test_a1])
 @pytest.mark.parametrize("test_a0", [test_a0])
 @pytest.mark.parametrize("validate_cp", [False, True])
 @pytest.mark.parametrize("order", ["row", "column"])
 def test_choi_to_kraus(
-    backend, order, validate_cp, test_a0, test_a1, test_kraus_left, test_kraus_right
+    backend,
+    order,
+    validate_cp,
+    test_a0,
+    test_a1,
+    test_kraus_left,
+    test_kraus_right,
+    test_non_CP,
 ):
     axes = [1, 2] if order == "row" else [0, 3]
     test_choi = backend.cast(
@@ -425,6 +433,7 @@ def test_choi_to_kraus(
     backend.assert_allclose(evolution_a1, test_evolution_a1, atol=2 * PRECISION_TOL)
 
     if validate_cp and order == "row":
+        test_non_CP = backend.cast(test_non_CP, dtype=test_non_CP.dtype)
         (kraus_left, kraus_right), _ = choi_to_kraus(
             test_non_CP, order=order, validate_cp=validate_cp, backend=backend
         )