Merge pull request #1451 from qiboteam/entanglement

Extend `meyer_wallach_entanglement` to `ndarrays`

src/qibo/quantum_info/entanglement.py

@@ -193,58 +193,87 @@ def entanglement_fidelity(
     return entang_fidelity
 
 
-def meyer_wallach_entanglement(circuit, backend=None):
-    """Computes the Meyer-Wallach entanglement Q of the `circuit`,
+def meyer_wallach_entanglement(state, backend=None):
+    """Compute the Meyer-Wallach entanglement :math:`Q` of a ``state``,
 
     .. math::
-        Q(\\theta) = 1 - \\frac{1}{N} \\, \\sum_{k} \\,
-            \\text{tr}\\left(\\rho_{k^{2}}(\\theta)\\right) \\, .
+        Q(\\rho) = 2\\left(1 - \\frac{1}{N} \\, \\sum_{k} \\,
+            \\text{tr}\\left(\\rho_{k}^{2}\\right)\\right) \\, ,
+
+    where :math:`\\rho_{k}^{2}` is the reduced density matrix of qubit :math:`k`,
+    and :math:`N` is the total number of qubits in ``state``.
+    We use the definition of the Meyer-Wallach entanglement as the average purity
+    proposed in `Brennen (2003) <https://dl.acm.org/doi/10.5555/2011556.2011561>`_,
+    which is equivalent to the definition introduced in `Meyer and Wallach (2002)
+    <https://doi.org/10.1063/1.1497700>`_.
 
     Args:
-        circuit (:class:`qibo.models.Circuit`): Parametrized circuit.
+        state (ndarray): statevector or density matrix.
         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: Meyer-Wallach entanglement.
+        float: Meyer-Wallach entanglement :math:`Q`.
+
+
+    References:
+        1. G. K. Brennen, *An observable measure of entanglement for pure states of
+        multi-qubit systems*, `Quantum Information and Computation, vol. 3 (6), 619-626
+        <https://dl.acm.org/doi/10.5555/2011556.2011561>`_ (2003).
+
+        2. D. A. Meyer and N. R. Wallach, *Global entanglement in multiparticle systems*,
+        `J. Math. Phys. 43, 4273–4278 <https://doi.org/10.1063/1.1497700>`_ (2002).
     """
 
     backend = _check_backend(backend)
 
-    circuit.density_matrix = True
-    nqubits = circuit.nqubits
+    if (
+        (len(state.shape) not in [1, 2])
+        or (len(state) == 0)
+        or (len(state.shape) == 2 and state.shape[0] != state.shape[1])
+    ):
+        raise_error(
+            TypeError,
+            f"state must have dims either (k,) or (k,k), but have dims {state.shape}.",
+        )
 
-    rho = backend.execute_circuit(circuit).state()
+    nqubits = int(np.log2(state.shape[-1]))
 
-    ent = 0
+    entanglement = 0
     for j in range(nqubits):
         trace_q = list(range(nqubits))
         trace_q.pop(j)
 
-        rho_r = partial_trace(rho, trace_q, backend=backend)
+        rho_r = partial_trace(state, trace_q, backend=backend)
 
         trace = purity(rho_r, backend=backend)
 
-        ent += trace
+        entanglement += trace
 
-    entanglement = 1 - ent / nqubits
-
-    return entanglement
+    return 2 * (1 - entanglement / nqubits)
 
 
 def entangling_capability(circuit, samples: int, seed=None, backend=None):
-    """Returns the entangling capability :math:`\\text{Ent}` of a parametrized
-    circuit, which is average Meyer-Wallach entanglement Q of the circuit, i.e.
+    """Return the entangling capability :math:`\\text{Ent}` of a parametrized circuit.
+
+    It is defined as the average Meyer-Wallach entanglement :math:`Q`
+    (:func:`qibo.quantum_info.meyer_wallach_entanglement`) of the ``circuit``, i.e.
 
     .. math::
-        \\text{Ent} = \\frac{2}{S}\\sum_{k}Q_k \\, ,
+        \\text{Ent} = \\frac{2}{|\\mathcal{S}|}\\sum_{\\theta_{k} \\in \\mathcal{S}}
+            \\, Q(\\rho_{k}) \\, ,
 
-    where :math:`S` is the number of samples.
+    where :math:`\\mathcal{S}` is the set of sampled circuit parameters,
+    and :math:`\\rho_{k}` is the state prepared by the circuit with uniformily-sampled
+    parameters :math:`\\theta_{k}`.
+
+    .. note::
+        Currently, function does not work with ``circuit`` that contains noisy channels.
 
     Args:
         circuit (:class:`qibo.models.Circuit`): Parametrized circuit.
-        samples (int): number of samples to estimate the integral.
+        samples (int): number of sampled circuit parameter vectors :math:`|S|`
         seed (int or :class:`numpy.random.Generator`, optional): Either a generator of random
             numbers or a fixed seed to initialize a generator. If ``None``, initializes
             a generator with a random seed. Default: ``None``.
@@ -253,7 +282,7 @@ def entangling_capability(circuit, samples: int, seed=None, backend=None):
             Defaults to ``None``.
 
     Returns:
-        float: Entangling capability.
+        float: Entangling capability :math:`\\text{Ent}`.
     """
 
     if not isinstance(samples, int):
@@ -276,13 +305,12 @@ def entangling_capability(circuit, samples: int, seed=None, backend=None):
         np.random.default_rng(seed) if seed is None or isinstance(seed, int) else seed
     )
 
-    res = []
+    capability = []
     for _ in range(samples):
         params = local_state.uniform(-np.pi, np.pi, circuit.trainable_gates.nparams)
         circuit.set_parameters(params)
-        entanglement = meyer_wallach_entanglement(circuit, backend=backend)
-        res.append(entanglement)
-
-    capability = 2 * np.real(np.sum(res)) / samples
+        state = backend.execute_circuit(circuit).state()
+        entanglement = meyer_wallach_entanglement(state, backend=backend)
+        capability.append(entanglement)
 
-    return capability
+    return 2 * np.real(np.sum(capability)) / samples

src/qibo/quantum_info/linalg_operations.py

@@ -149,8 +149,8 @@ def partial_trace(state, traced_qubits: Union[List[int], Tuple[int]], backend=No
         return backend.np.reshape(rho, shape)
 
     order = tuple(sorted(traced_qubits))
-    order += tuple(i for i in range(nqubits) if i not in traced_qubits)
-    order += tuple(i + nqubits for i in order)
+    order += tuple(set(list(range(nqubits))) ^ set(traced_qubits))
+    order += tuple(k + nqubits for k in order)
     shape = 2 * (2 ** len(traced_qubits), 2 ** (nqubits - len(traced_qubits)))
 
     state = backend.np.transpose(state, order)

tests/test_quantum_info_entanglement.py

@@ -126,21 +126,28 @@ def test_entanglement_fidelity(backend, channel, nqubits, check_hermitian):
 
 
 def test_meyer_wallach_entanglement(backend):
+    with pytest.raises(TypeError):
+        state = np.random.rand(2, 3, 2)
+        state = backend.cast(state, dtype=state.dtype)
+        test = meyer_wallach_entanglement(state, backend=backend)
+
     nqubits = 2
 
     circuit1 = Circuit(nqubits)
     circuit1.add([gates.RX(0, np.pi / 4)] for _ in range(nqubits))
+    state1 = backend.execute_circuit(circuit1).state()
 
     circuit2 = Circuit(nqubits)
     circuit2.add([gates.RX(0, np.pi / 4)] for _ in range(nqubits))
     circuit2.add(gates.CNOT(0, 1))
+    state2 = backend.execute_circuit(circuit2).state()
 
     backend.assert_allclose(
-        meyer_wallach_entanglement(circuit1, backend=backend), 0.0, atol=PRECISION_TOL
+        meyer_wallach_entanglement(state1, backend=backend), 0.0, atol=PRECISION_TOL
     )
 
     backend.assert_allclose(
-        meyer_wallach_entanglement(circuit2, backend=backend), 0.5, atol=PRECISION_TOL
+        meyer_wallach_entanglement(state2, backend=backend), 1, atol=PRECISION_TOL
     )