diff --git a/doc/source/api-reference/qibo.rst b/doc/source/api-reference/qibo.rst index 191ea01a19..b2b1e9075d 100644 --- a/doc/source/api-reference/qibo.rst +++ b/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 ^^^^^^^^^^^^^^^^ diff --git a/src/qibo/backends/abstract.py b/src/qibo/backends/abstract.py index 468fa89008..4b168cf0ef 100644 --- a/src/qibo/backends/abstract.py +++ b/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 diff --git a/src/qibo/backends/numpy.py b/src/qibo/backends/numpy.py index fe1441b7c6..e9e57ca41d 100644 --- a/src/qibo/backends/numpy.py +++ b/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 diff --git a/src/qibo/backends/pytorch.py b/src/qibo/backends/pytorch.py index 22fbbc5640..1a8051ab4e 100644 --- a/src/qibo/backends/pytorch.py +++ b/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 diff --git a/src/qibo/backends/tensorflow.py b/src/qibo/backends/tensorflow.py index 2088cd69be..13a4a5bd25 100644 --- a/src/qibo/backends/tensorflow.py +++ b/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( diff --git a/src/qibo/quantum_info/linalg_operations.py b/src/qibo/quantum_info/linalg_operations.py index 0521af97f0..ed23e89b2e 100644 --- a/src/qibo/quantum_info/linalg_operations.py +++ b/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) diff --git a/src/qibo/quantum_info/superoperator_transformations.py b/src/qibo/quantum_info/superoperator_transformations.py index 158dd73ac4..7d74a0ef76 100644 --- a/src/qibo/quantum_info/superoperator_transformations.py +++ b/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): diff --git a/tests/test_quantum_info_operations.py b/tests/test_quantum_info_operations.py index 4d25e5c4a6..9d9957f205 100644 --- a/tests/test_quantum_info_operations.py +++ b/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) diff --git a/tests/test_quantum_info_superoperator_transformations.py b/tests/test_quantum_info_superoperator_transformations.py index d93db2e5e6..8ac14aa1fd 100644 --- a/tests/test_quantum_info_superoperator_transformations.py +++ b/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,6 +377,7 @@ 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]) @@ -384,7 +385,14 @@ def test_choi_to_pauli(backend, normalize, order, pauli_order, test_superop): @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 )