Merge branch 'master' into inner_product

.pre-commit-config.yaml

@@ -2,14 +2,15 @@ ci:
   autofix_prs: true
 repos:
   - repo: https://github.com/pre-commit/pre-commit-hooks
-    rev: v4.6.0
+    rev: v5.0.0
     hooks:
       - id: trailing-whitespace
       - id: end-of-file-fixer
       - id: check-yaml
+      - id: check-toml
       - id: debug-statements
   - repo: https://github.com/psf/black
-    rev: 24.8.0
+    rev: 24.10.0
     hooks:
       - id: black
   - repo: https://github.com/pycqa/isort

doc/source/api-reference/qibo.rst

@@ -1704,6 +1704,12 @@ Entanglement of formation
 .. autofunction:: qibo.quantum_info.entanglement_of_formation
 
 
+Negativity
+""""""""""
+
+.. autofunction:: qibo.quantum_info.negativity
+
+
 Entanglement fidelity
 """""""""""""""""""""
 
@@ -1987,6 +1993,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/basis.py

@@ -1,4 +1,3 @@
-from functools import reduce
 from itertools import product
 from typing import Optional
 
@@ -92,43 +91,41 @@ def pauli_basis(
     backend = _check_backend(backend)
 
     pauli_labels = {"I": matrices.I, "X": matrices.X, "Y": matrices.Y, "Z": matrices.Z}
-    basis_single = [pauli_labels[label] for label in pauli_order]
+    dim = 2**nqubits
+    basis_single = backend.cast([pauli_labels[label] for label in pauli_order])
+    einsum = np.einsum if backend.name == "tensorflow" else backend.np.einsum
 
     if nqubits > 1:
-        basis_full = list(product(basis_single, repeat=nqubits))
-        basis_full = [reduce(np.kron, row) for row in basis_full]
+        input_indices = [range(3 * i, 3 * (i + 1)) for i in range(nqubits)]
+        output_indices = (i for indices in zip(*input_indices) for i in indices)
+        operands = [basis_single for _ in range(nqubits)]
+        inputs = [item for pair in zip(operands, input_indices) for item in pair]
+        basis_full = einsum(*inputs, output_indices).reshape(4**nqubits, dim, dim)
     else:
         basis_full = basis_single
 
-    basis_full = backend.cast(basis_full, dtype=basis_full[0].dtype)
-
     if vectorize and sparse:
-        basis, indexes = [], []
-        for row in basis_full:
-            row = vectorization(row, order=order, backend=backend)
-            row_indexes = backend.np.flatnonzero(row)
-            indexes.append(row_indexes)
-            basis.append(row[row_indexes])
-            del row
+        if backend.name == "tensorflow":
+            nonzero = np.nonzero
+        elif backend.name == "pytorch":
+            nonzero = lambda x: backend.np.nonzero(x, as_tuple=True)
+        else:
+            nonzero = backend.np.nonzero
+        basis = vectorization(basis_full, order=order, backend=backend)
+        indices = nonzero(basis)
+        basis = basis[indices].reshape(-1, dim)
+        indices = indices[1].reshape(-1, dim)
+
     elif vectorize and not sparse:
-        basis = [
-            vectorization(
-                backend.cast(matrix, dtype=matrix.dtype), order=order, backend=backend
-            )
-            for matrix in basis_full
-        ]
+        basis = vectorization(basis_full, order=order, backend=backend)
     else:
         basis = basis_full
 
-    basis = backend.cast(basis, dtype=basis[0].dtype)
-
     if normalize:
         basis = basis / np.sqrt(2**nqubits)
 
     if vectorize and sparse:
-        indexes = backend.cast(indexes, dtype=indexes[0][0].dtype)
-
-        return basis, indexes
+        return basis, indices
 
     return basis
 

src/qibo/quantum_info/entanglement.py

@@ -4,7 +4,11 @@
 
 from qibo.backends import _check_backend
 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,
+    partial_transpose,
+)
 from qibo.quantum_info.metrics import fidelity, purity
 
 
@@ -116,6 +120,39 @@ def entanglement_of_formation(
     return ent_of_form
 
 
+def negativity(state, bipartition, backend=None):
+    """Calculates the negativity of a bipartite quantum state.
+
+    Given a bipartite state :math:`\\rho \\in \\mathcal{H}_{A} \\otimes \\mathcal{H}_{B}`,
+    the negativity :math:`\\operatorname{Neg}(\\rho)` is given by
+
+    .. math::
+        \\operatorname{Neg}(\\rho) = \\frac{1}{2} \\,
+            \\left( \\norm{\\rho_{B}}_{1} - 1 \\right) \\, ,
+
+    where :math:`\\rho_{B}` is the reduced density matrix after tracing out qubits in
+    partition :math:`A`, and :math:`\\norm{\\cdot}_{1}` is the Schatten :math:`1`-norm
+    (also known as nuclear norm or trace norm).
+
+    Args:
+        state (ndarray): statevector or density matrix.
+        bipartition (list or tuple or ndarray): qubits in the subsystem to be traced out.
+        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:
+        float: Negativity :math:`\\operatorname{Neg}(\\rho)` of state :math:`\\rho`.
+    """
+    backend = _check_backend(backend)
+
+    reduced = partial_transpose(state, bipartition, backend)
+    reduced = backend.np.conj(reduced.T) @ reduced
+    norm = backend.np.trace(matrix_power(reduced, 1 / 2, backend))
+
+    return float(backend.np.real((norm - 1) / 2))
+
+
 def entanglement_fidelity(
     channel, nqubits: int, state=None, check_hermitian: bool = False, backend=None
 ):

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):
@@ -26,8 +27,9 @@ def vectorization(state, order: str = "row", backend=None):
     .. math::
         |\\rho) = \\sum_{k, l} \\, \\rho_{kl} \\, \\ket{l} \\otimes \\ket{k}
 
+    If ``state`` is a 3-dimensional tensor it is interpreted as a batch of states.
     Args:
-        state: state vector or density matrix.
+        state: statevector, density matrix, an array of statevectors, or an array of density matrices.
         order (str, optional): If ``"row"``, vectorization is performed
             row-wise. If ``"column"``, vectorization is performed
             column-wise. If ``"system"``, a block-vectorization is
@@ -40,13 +42,13 @@ def vectorization(state, order: str = "row", backend=None):
         ndarray: Liouville representation of ``state``.
     """
     if (
-        (len(state.shape) >= 3)
+        (len(state.shape) > 3)
         or (len(state) == 0)
         or (len(state.shape) == 2 and state.shape[0] != state.shape[1])
     ):
         raise_error(
             TypeError,
-            f"Object must have dims either (k,) or (k,k), but have dims {state.shape}.",
+            f"Object must have dims either (k,), (k, k), (N, 1, k) or (N, k, k), but have dims {state.shape}.",
         )
 
     if not isinstance(order, str):
@@ -62,25 +64,36 @@ def vectorization(state, order: str = "row", backend=None):
 
     backend = _check_backend(backend)
 
+    dims = state.shape[-1]
+
     if len(state.shape) == 1:
         state = backend.np.outer(state, backend.np.conj(state))
+    elif len(state.shape) == 3 and state.shape[1] == 1:
+        state = backend.np.einsum(
+            "aij,akl->aijkl", state, backend.np.conj(state)
+        ).reshape(state.shape[0], dims, dims)
 
     if order == "row":
-        state = backend.np.reshape(state, (1, -1))[0]
+        state = backend.np.reshape(state, (-1, dims**2))
     elif order == "column":
-        state = state.T
-        state = backend.np.reshape(state, (1, -1))[0]
+        indices = list(range(len(state.shape)))
+        indices[-2:] = reversed(indices[-2:])
+        state = backend.np.transpose(state, indices)
+        state = backend.np.reshape(state, (-1, dims**2))
     else:
-        dim = len(state)
-        nqubits = int(np.log2(dim))
+        nqubits = int(np.log2(state.shape[-1]))
 
-        new_axis = []
+        new_axis = [0]
         for qubit in range(nqubits):
-            new_axis += [qubit + nqubits, qubit]
+            new_axis.extend([qubit + nqubits + 1, qubit + 1])
 
-        state = backend.np.reshape(state, [2] * 2 * nqubits)
+        state = backend.np.reshape(state, [-1] + [2] * 2 * nqubits)
         state = backend.np.transpose(state, new_axis)
-        state = backend.np.reshape(state, (-1,))
+        state = backend.np.reshape(state, (-1, 2 ** (2 * nqubits)))
+
+    state = backend.np.squeeze(
+        state, axis=tuple(i for i, ax in enumerate(state.shape) if ax == 1)
+    )
 
     return state
 
@@ -483,10 +496,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):