Implement and test process fidelity

qutip/core/metrics.py

@@ -43,7 +43,8 @@
 import numpy as np
 from scipy import linalg as la
 import scipy.sparse as sp
-from .superop_reps import to_kraus, to_choi, _to_superpauli, to_super
+from .superop_reps import (to_kraus, to_choi, _to_superpauli, to_super,
+                           kraus_to_choi)
 from .superoperator import operator_to_vector, vector_to_operator
 from .operators import qeye
 from .semidefinite import dnorm_problem
@@ -105,12 +106,104 @@ def fidelity(A, B):
     return float(np.real(np.sqrt(eig_vals[eig_vals > 0]).sum()))
 
 
-def process_fidelity(U1, U2, normalize=True):
+def _process_fidelity_to_id(oper):
     """
-    Calculate the process fidelity given two process operators.
+    Internal function returning the process fidelity of a quantum channel
+    to the identity quantum channel.
+    Parameters
+    ----------
+    oper : :class:`qutip.Qobj`/list
+        A unitary operator, or a superoperator in supermatrix, Choi or
+        chi-matrix form, or a list of Kraus operators
+    Returns
+    -------
+    fid : float
+    """
+    msg = 'The process fidelity to identity is only defined for ' \
+          'dimension preserving channels.'
+    if isinstance(oper, list):  # oper is a list of Kraus operators
+        d = oper[0].shape[0]
+        if oper[0].shape[1] != d:
+            raise TypeError(msg)
+        return np.sum([np.abs(k.tr()) ** 2 for k in oper]) / d ** 2
+    elif oper.type == 'oper':  # interpret as unitary
+        d = oper.shape[0]
+        if oper.shape[1] != d:
+            raise TypeError(msg)
+        return np.abs(oper.tr()) ** 2 / d ** 2
+    elif oper.type == 'super':
+        d = np.prod(oper.dims[0][0])
+        if np.prod(oper.dims[1][0]) != d:
+            raise TypeError(msg)
+        if oper.superrep == 'chi':
+            return oper[0, 0].real / d ** 2
+        else:  # oper.superrep is either 'super' or 'choi':
+            return to_super(oper).tr().real / d ** 2
+
+
+def _kraus_or_qobj_to_choi(oper):
+    if isinstance(oper, list):
+        return kraus_to_choi(oper)
+    else:
+        return to_choi(oper)
+
+
+def process_fidelity(oper, target=None):
     """
-    out = (U1 * U2).tr()
-    return out / (U1.tr() * U2.tr()) if normalize else out
+    Returns the process fidelity of a quantum channel to the target
+    channel, or to the identity channel if no target is given.
+    The process fidelity between two channels is defined as the state
+    fidelity between their normalized Choi matrices.
+    Parameters
+    ----------
+    oper : :class:`qutip.Qobj`/list
+        A unitary operator, or a superoperator in supermatrix, Choi or
+        chi-matrix form, or a list of Kraus operators
+    target : :class:`qutip.Qobj`/list
+        A unitary operator, or a superoperator in supermatrix, Choi or
+        chi-matrix form, or a list of Kraus operators
+    Returns
+    -------
+    fid : float
+        Process fidelity between oper and target,
+        or between oper and identity.
+    Notes
+    -----
+    See, for example: A. Gilchrist, N.K. Langford, M.A. Nielsen,
+    Phys. Rev. A 71, 062310 (2005).
+    The definition of state fidelity that the process fidelity is based on
+    is the one from R. Jozsa, Journal of Modern Optics, 41:12, 2315 (1994).
+    It is the square of the one implemented in
+    :func:`qutip.metrics.fidelity` which follows Nielsen & Chuang,
+    "Quantum Computation and Quantum Information"
+    """
+    if target is None:
+        return _process_fidelity_to_id(oper)
+    elif not isinstance(target, list) and target.type == 'oper':
+        # interpret target as unitary.
+        if isinstance(oper, list):  # oper is a list of Kraus operators
+            if oper[0].dims != target.dims:
+                raise TypeError('Dimensions of oper and target do not match')
+            return _process_fidelity_to_id([k * target.dag() for k in oper])
+        elif oper.type == 'oper':
+            if oper.dims != target.dims:
+                raise TypeError('Dimensions of oper and target do not match')
+            return _process_fidelity_to_id(oper*target.dag())
+        elif oper.type == 'super':
+            oper_super = to_super(oper)
+            target_dag_super = to_super(target.dag())
+            if oper_super.dims != target_dag_super.dims:
+                raise TypeError('Dimensions of oper and target do not match')
+            return _process_fidelity_to_id(oper_super * target_dag_super)
+    else:  # target is a list of Kraus operators or a superoperator
+        if not isinstance(oper, list) and oper.type == 'oper':
+            return process_fidelity(target, oper)  # reverse order
+        oper_choi = _kraus_or_qobj_to_choi(oper)
+        target_choi = _kraus_or_qobj_to_choi(target)
+        if oper_choi.dims != target_choi.dims:
+            raise TypeError('Dimensions of oper and target do not match')
+        d = np.prod(oper_choi.dims[0][0])
+        return fidelity(oper_choi / d, target_choi / d)**2
 
 
 def average_gate_fidelity(oper, target=None):

qutip/tests/core/test_metrics.py

@@ -52,8 +52,8 @@
     rand_ket_haar, rand_dm_ginibre, rand_unitary_haar
 )
 from qutip import (
-    Qobj, to_super, to_choi, tensor, create, destroy, jmat, identity, qdiags,
-    sigmax, sigmay, sigmaz, qeye, fock_dm, basis,
+    Qobj, to_super, to_choi, to_chi, to_kraus, tensor, create, destroy, jmat,
+    identity, qdiags, sigmax, sigmay, sigmaz, qeye, fock_dm, basis,
 )
 from qutip.core.metrics import *
 from qutip.qip.operations.gates import hadamard_transform, swap
@@ -204,6 +204,67 @@ def test_fidelity_overlap():
             np.abs(psi.dag() * phi)
         )
 
+
+def test_process_fidelity_of_identity():
+    """
+    Metrics: process fidelity of identity map is 1
+    """
+    num_qubits = 3
+    oper = qeye(num_qubits*[2])
+    for superrep_conversion in [
+            lambda x:x, to_super, to_choi, to_chi, to_kraus]:
+        f = process_fidelity(superrep_conversion(oper))
+        assert_(np.isrealobj(f))
+        assert_almost_equal(f, 1)
+
+
+def test_process_fidelity_identical_maps():
+    """
+    Metrics: process fidelity of a map to itself is 1
+    """
+    num_qubits = 2
+    for k in range(10):
+        oper = rand_unitary(2**num_qubits, dims=2*[num_qubits*[2]])
+        f = process_fidelity(oper, oper)
+        assert_almost_equal(f, 1)
+        oper = rand_super_bcsz(2**num_qubits, dims=2*[2*[num_qubits*[2]]])
+        for superrep_conversion in [
+                lambda x: x, to_choi, to_chi, to_kraus]:
+            oper_converted = superrep_conversion(oper)
+            f = process_fidelity(oper_converted, oper_converted)
+            assert_almost_equal(f, 1)
+
+
+def test_process_fidelity_consistency():
+    """
+    Metrics: process fidelity independent of channel representation
+    """
+    num_qubits = 2
+    for k in range(10):
+        fidelities_u_to_u = []
+        fidelities_u_to_id = []
+        u1 = rand_unitary(2**num_qubits, dims=2*[num_qubits*[2]])
+        u2 = rand_unitary(2**num_qubits, dims=2*[num_qubits*[2]])
+        for map1 in [lambda x:x, to_super, to_choi, to_chi, to_kraus]:
+            fidelities_u_to_id.append(process_fidelity(map1(u1)))
+            for map2 in [lambda x:x, to_super, to_choi, to_chi, to_kraus]:
+                fidelities_u_to_u.append(process_fidelity(map1(u1), map2(u2)))
+        assert_almost_equal(fidelities_u_to_id, fidelities_u_to_id[0])
+        assert_almost_equal(fidelities_u_to_u, fidelities_u_to_u[0])
+
+
+def test_process_fidelity_unitary_invariance():
+    """
+    Metrics: process fidelity, invariance under unitary trans.
+    """
+    for k in range(10):
+        op1 = rand_super_bcsz(10)
+        op2 = rand_super_bcsz(10)
+        u = to_super(rand_unitary(10))
+        assert_almost_equal(process_fidelity(op1, op2),
+                            process_fidelity(u*op1, u*op2))
+
+
 def test_tracedist1():
     """
     Metrics: Trace dist., invariance under unitary trans.