Merge pull request #785 from qiboteam/liouville

Choi representation of channels + transformations between representations

doc/source/api-reference/qibo.rst

@@ -1018,18 +1018,6 @@ Basis
 Set of functions related to basis and basis transformations.
 
 
-Vectorization
-"""""""""""""
-
-.. autofunction:: qibo.quantum_info.vectorization
-
-
-Unvectorization
-"""""""""""""""
-
-.. autofunction:: qibo.quantum_info.unvectorization
-
-
 Pauli basis
 """""""""""
 
@@ -1168,6 +1156,98 @@ Random stochastic matrix
 .. autofunction:: qibo.quantum_info.random_stochastic_matrix
 
 
+Superoperator Transformations
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Functions used to convert superoperators among their possible representations.
+For more in-depth theoretical description of the representations and transformations,
+we direct the reader to `Wood, Biamonte, and Cory, Quant. Inf. Comp. 15, 0579-0811 (2015) <https://arxiv.org/abs/1111.6950>`_.
+
+
+Vectorization
+"""""""""""""
+
+.. autofunction:: qibo.quantum_info.vectorization
+
+.. note::
+    Due to `numpy` limitations on handling transposition of tensors,
+    this function will not work when the number of qubits :math:`n`
+    is such that :math:`n > 16`.
+
+
+Unvectorization
+"""""""""""""""
+
+.. autofunction:: qibo.quantum_info.unvectorization
+
+.. note::
+    Due to `numpy` limitations on handling transposition of tensors,
+    this function will not work when the number of qubits :math:`n`
+    is such that :math:`n > 16`.
+
+
+Liouville to Choi
+"""""""""""""""""
+
+.. autofunction:: qibo.quantum_info.liouville_to_choi
+
+
+Choi to Liouville
+"""""""""""""""""
+
+.. autofunction:: qibo.quantum_info.choi_to_liouville
+
+
+Choi to Kraus
+"""""""""""""
+
+.. autofunction:: qibo.quantum_info.superoperator_transformations.choi_to_kraus
+
+.. note::
+    Due to the spectral decomposition subroutine in this function, the resulting Kraus
+    operators :math:`\{K_{\alpha}\}_{\alpha}` might contain global phases. That
+    implies these operators are not exactly equal to the "true" Kraus operators
+    :math:`\{K_{\alpha}^{(\text{ideal})}\}_{\alpha}`. However, since these are
+    global phases, the operators' actions are the same, i.e.
+
+    .. math::
+        K_{\alpha} \, \rho \, K_{\alpha}^{\dagger} = K_{\alpha}^{\text{(ideal)}} \, \rho \,\, (K_{\alpha}^{\text{(ideal)}})^{\dagger} \,\,\,\,\, , \,\, \forall \, \alpha
+
+.. note::
+    User can set ``validate_CP=False`` in order to speed up execution by not checking if
+    input map ``choi_super_op`` is completely positive (CP) and Hermitian. However, that may
+    lead to erroneous outputs if ``choi_super_op`` is not guaranteed to be CP. We advise users
+    to either set this flag carefully or leave it in its default setting (``validate_CP=True``).
+
+
+Kraus to Choi
+"""""""""""""
+
+.. autofunction:: qibo.quantum_info.kraus_to_choi
+
+
+Kraus to Liouville
+""""""""""""""""""
+
+.. autofunction:: qibo.quantum_info.kraus_to_liouville
+
+
+Liouville to Kraus
+""""""""""""""""""
+
+.. autofunction:: qibo.quantum_info.liouville_to_kraus
+
+.. note::
+    Due to the spectral decomposition subroutine in this function, the resulting Kraus
+    operators :math:`\{K_{\alpha}\}_{\alpha}` might contain global phases. That
+    implies these operators are not exactly equal to the "true" Kraus operators
+    :math:`\{K_{\alpha}^{(\text{ideal})}\}_{\alpha}`. However, since these are
+    global phases, the operators' actions are the same, i.e.
+
+    .. math::
+        K_{\alpha} \, \rho \, K_{\alpha}^{\dagger} = K_{\alpha}^{\text{(ideal)}} \, \rho \,\, (K_{\alpha}^{\text{(ideal)}})^{\dagger} \,\,\,\,\, , \,\, \forall \, \alpha
+
+
 Utility Functions
 ^^^^^^^^^^^^^^^^^
 

src/qibo/gates/channels.py

@@ -101,19 +101,21 @@ def __init__(self, ops):
         self.coefficients = len(self.gates) * (1,)
         self.coefficient_sum = 1
 
-    def to_superop(self, backend=None):
-        """Returns the Liouville representation of the Kraus channel.
+    def to_choi(self, backend=None):
+        """Returns the Choi representation of the Kraus channel.
 
         Args:
-            backend (``qibo.backends.abstract.Backend``, optional): backend
-                to be used in the execution. If ``None``, it uses
-                ``GlobalBackend()``. Defaults to ``None``.
+            backend (``qibo.backends.abstract.Backend``, optional): backend to be
+                used in the execution. If ``None``, it uses ``GlobalBackend()``.
+                Defaults to ``None``.
 
         Returns:
-            Liouville representation of the channel.
+            Choi representation of the Kraus channel.
         """
         import numpy as np
 
+        from qibo.quantum_info.superoperator_transformations import vectorization
+
         if backend is None:  # pragma: no cover
             from qibo.backends import GlobalBackend
 
@@ -129,13 +131,41 @@ def to_superop(self, backend=None):
             self.gates += (I(*self.target_qubits),)
 
         super_op = np.zeros((4**self.nqubits, 4**self.nqubits), dtype="complex")
-        super_op = backend.cast(super_op, dtype=super_op.dtype)
         for coeff, gate in zip(self.coefficients, self.gates):
             kraus_op = FusedGate(*range(self.nqubits))
             kraus_op.append(gate)
             kraus_op = kraus_op.asmatrix(backend)
-            kraus_op = coeff * np.kron(np.conj(kraus_op), kraus_op)
-            super_op += backend.cast(kraus_op, dtype=kraus_op.dtype)
+            kraus_op = vectorization(coeff * kraus_op)
+            super_op += np.outer(kraus_op, np.conj(kraus_op))
+            del kraus_op
+
+        super_op = backend.cast(super_op, dtype=super_op.dtype)
+
+        return super_op
+
+    def to_superop(self, backend=None):
+        """Returns the Liouville representation of the Kraus channel.
+
+        Args:
+            backend (``qibo.backends.abstract.Backend``, optional): backend
+                to be used in the execution. If ``None``, it uses
+                ``GlobalBackend()``. Defaults to ``None``.
+
+        Returns:
+            Liouville representation of the channel.
+        """
+        import numpy as np
+
+        from qibo.quantum_info.superoperator_transformations import choi_to_liouville
+
+        if backend is None:  # pragma: no cover
+            from qibo.backends import GlobalBackend
+
+            backend = GlobalBackend()
+
+        super_op = self.to_choi(backend=backend)
+        super_op = choi_to_liouville(super_op)
+        super_op = backend.cast(super_op, dtype=super_op.dtype)
 
         return super_op
 
@@ -162,7 +192,23 @@ def to_pauli_liouville(self, normalize: bool = False, backend=None):
 
             backend = GlobalBackend()
 
-        super_op = self.to_superop(backend=backend)
+        self.nqubits = 1 + max(self.target_qubits)
+
+        if self.name != "KrausChannel":
+            p0 = 1
+            for coeff in self.coefficients:
+                p0 = p0 - coeff
+            self.coefficients += (p0,)
+            self.gates += (I(*self.target_qubits),)
+
+        super_op = np.zeros((4**self.nqubits, 4**self.nqubits), dtype="complex")
+        super_op = backend.cast(super_op, dtype=super_op.dtype)
+        for coeff, gate in zip(self.coefficients, self.gates):
+            kraus_op = FusedGate(*range(self.nqubits))
+            kraus_op.append(gate)
+            kraus_op = kraus_op.asmatrix(backend)
+            kraus_op = coeff * np.kron(np.conj(kraus_op), kraus_op)
+            super_op += backend.cast(kraus_op, dtype=kraus_op.dtype)
 
         # unitary that transforms from comp basis to pauli basis
         U = backend.cast(comp_basis_to_pauli(self.nqubits, normalize))

src/qibo/quantum_info/__init__.py

@@ -1,4 +1,5 @@
 from qibo.quantum_info.basis import *
 from qibo.quantum_info.metrics import *
 from qibo.quantum_info.random_ensembles import *
+from qibo.quantum_info.superoperator_transformations import *
 from qibo.quantum_info.utils import *

src/qibo/quantum_info/basis.py

@@ -5,114 +5,7 @@
 
 from qibo import matrices
 from qibo.config import raise_error
-
-
-def vectorization(state, order: str = "row"):
-    """Returns state :math:`\\rho` in its Liouville
-    representation :math:`\\ket{\\rho}`.
-
-    Args:
-        state: state vector or density matrix.
-        order (str, optional): If ``"row"``, vectorization is performed
-            row-wise. If ``"column"``, vectorization is performed
-            column-wise. If ``"system"``, a block-vectorization is
-            performed. Default is ``"row"``.
-
-    Returns:
-        Liouville representation of ``state``.
-    """
-
-    if (
-        (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}.",
-        )
-
-    if not isinstance(order, str):
-        raise_error(
-            TypeError, f"order must be type str, but it is type {type(order)} instead."
-        )
-    else:
-        if (order != "row") and (order != "column") and (order != "system"):
-            raise_error(
-                ValueError,
-                f"order must be either 'row' or 'column' or 'system', but it is {order}.",
-            )
-
-    if len(state.shape) == 1:
-        state = np.outer(state, np.conj(state))
-
-    if order == "row":
-        state = np.reshape(state, (1, -1), order="C")[0]
-    elif order == "column":
-        state = np.reshape(state, (1, -1), order="F")[0]
-    else:
-        d = len(state)
-        nqubits = int(np.log2(d))
-
-        new_axis = []
-        for x in range(nqubits):
-            new_axis += [x + nqubits, x]
-        state = np.reshape(
-            np.transpose(np.reshape(state, [2] * 2 * nqubits), axes=new_axis), -1
-        )
-
-    return state
-
-
-def unvectorization(state, order: str = "row"):
-    """Returns state :math:`\\rho` from its Liouville
-    representation :math:`\\ket{\\rho}`.
-
-    Args:
-        state: :func:`vectorization` of a quantum state.
-        order (str, optional): If ``"row"``, unvectorization is performed
-            row-wise. If ``"column"``, unvectorization is performed
-            column-wise. If ``"system"``, system-wise vectorization is
-            performed. Default is ``"row"``.
-
-    Returns:
-        Density matrix of ``state``.
-    """
-
-    if len(state.shape) != 1:
-        raise_error(
-            TypeError,
-            f"Object must have dims (k,), but have dims {state.shape}.",
-        )
-
-    if not isinstance(order, str):
-        raise_error(
-            TypeError, f"order must be type str, but it is type {type(order)} instead."
-        )
-    else:
-        if (order != "row") and (order != "column") and (order != "system"):
-            raise_error(
-                ValueError,
-                f"order must be either 'row' or 'column' or 'system', but it is {order}.",
-            )
-
-    d = int(np.sqrt(len(state)))
-
-    if (order == "row") or (order == "column"):
-        order = "C" if order == "row" else "F"
-        state = np.reshape(state, (d, d), order=order)
-    else:
-        nqubits = int(np.log2(d))
-        axes_old = list(np.arange(0, 2 * nqubits))
-        state = np.reshape(
-            np.transpose(
-                np.reshape(state, [2] * 2 * nqubits),
-                axes=axes_old[1::2] + axes_old[0::2],
-            ),
-            [2**nqubits] * 2,
-        )
-
-    return state
+from qibo.quantum_info.superoperator_transformations import vectorization
 
 
 def pauli_basis(