Python cvxpy.trace() Examples

def tracenorm(A):
    """
    Compute the trace norm of matrix A given by:

      Tr( sqrt{ A^dagger * A } )

    Parameters
    ----------
    A : numpy array
        The matrix to compute the trace norm of.
    """
    if _np.linalg.norm(A - _np.conjugate(A.T)) < 1e-8:
        #Hermitian, so just sum eigenvalue magnitudes
        return _np.sum(_np.abs(_np.linalg.eigvals(A)))
    else:
        #Sum of singular values (positive by construction)
        return _np.sum(_np.linalg.svd(A, compute_uv=False)) 

def smith_fidelity(rho: np.ndarray, sigma: np.ndarray, power) -> float:
    r"""
    Computes the Smith fidelity :math:`F_S(\rho, \sigma, power)` between two quantum states rho and
    sigma.

    The Smith fidelity is  defined as :math:`F_S = \sqrt{F^{power}}`, where F is  standard fidelity
    :math:`F = fidelity(\rho, \sigma)`. Since the power is only defined for values less than 2,
    it is always true that :math:`F_S > F`.

    At present there is no known operational interpretation of the Smith fidelity for an arbitrary
    power.

    :param rho: Is a dim by dim positive matrix with unit trace.
    :param sigma: Is a dim by dim positive matrix with unit trace.
    :param power: Is a positive scalar less than 2.
    :return: Smith Fidelity which is a scalar.
    """
    if power < 0:
        raise ValueError("Power must be positive")
    if power >= 2:
        raise ValueError("Power must be less than 2")
    return np.sqrt(fidelity(rho, sigma)) ** power 

def jtracedist(A, B, mxBasis='pp'):  # Jamiolkowski trace distance:  Tr(|J(A)-J(B)|)
    """
    Compute the Jamiolkowski trace distance between operation matrices A and B,
    given by:

      D = 0.5 * Tr( sqrt{ (J(A)-J(B))^2 } )

    where J(.) is the Jamiolkowski isomorphism map that maps a operation matrix
    to it's corresponding Choi Matrix.

    Parameters
    ----------
    A, B : numpy array
        The matrices to compute the distance between.

    mxBasis : {'std', 'gm', 'pp', 'qt'} or Basis object
        The source and destination basis, respectively.  Allowed
        values are Matrix-unit (std), Gell-Mann (gm), Pauli-product (pp),
        and Qutrit (qt) (or a custom basis object).
    """
    JA = _jam.fast_jamiolkowski_iso_std(A, mxBasis)
    JB = _jam.fast_jamiolkowski_iso_std(B, mxBasis)
    return tracedist(JA, JB) 

def hilbert_schmidt_ip(A: np.ndarray, B: np.ndarray, tol: float = 1000) -> float:
    r"""
    Computes the Hilbert-Schmidt (HS) inner product between two operators A and B.

    This inner product is defined as

    .. math::

        HS = (A|B) = Tr[A^\dagger B]

    where :math:`|B) = vec(B)` and :math:`(A|` is the dual vector to :math:`|A)`.

    :param A: Is a dim by dim positive matrix with unit trace.
    :param B: Is a dim by dim positive matrix with unit trace.
    :param tol: Tolerance in machine epsilons for np.real_if_close.
    :return: HS inner product which is a scalar.
    """
    hs_ip = np.trace(np.matmul(np.transpose(np.conj(A)), B))
    return np.ndarray.item(np.real_if_close(hs_ip, tol)) 

def bures_angle(rho: np.ndarray, sigma: np.ndarray) -> float:
    r"""
    Computes the Bures angle (AKA Bures arc or Bures length) between two states rho and sigma:

    .. math::

        D_A(\rho, \sigma) = \arccos(\sqrt{F(\rho, \sigma)})

    where :math:`F(\rho, \sigma)` is the fidelity.

    The Bures angle is a measure of statistical distance between quantum states.

    :param rho: Is a dim by dim positive matrix with unit trace.
    :param sigma: Is a dim by dim positive matrix with unit trace.
    :return: Bures angle which is a scalar.
    """
    return np.arccos(np.sqrt(fidelity(rho, sigma))) 

def povm_jtracedist(model, targetModel, povmlbl):
    """
    Computes the Jamiolkowski trace distance between POVM maps using :func:`jtracedist`.

    Parameters
    ----------
    model, targetModel : Model
        Models containing the two POVMs to compare.

    povmlbl : str
        The POVM label

    Returns
    -------
    float
    """
    povm_mx = get_povm_map(model, povmlbl)
    target_povm_mx = get_povm_map(targetModel, povmlbl)
    return jtracedist(povm_mx, target_povm_mx, targetModel.basis) 

def fidelity(rho: np.ndarray, sigma: np.ndarray, tol: float = 1000) -> float:
    r"""
    Computes the fidelity :math:`F(\rho, \sigma)` between two quantum states rho and sigma.

    If the states are pure the expression reduces to

    .. math::

        F(|psi>,|phi>) = |<psi|phi>|^2

    The fidelity obeys :math:`0 ≤ F(\rho, \sigma) ≤ 1`, where
    :math:`F(\rho, \sigma)=1 iff \rho = \sigma`.

    :param rho: Is a dim by dim positive matrix with unit trace.
    :param sigma: Is a dim by dim positive matrix with unit trace.
    :param tol: Tolerance in machine epsilons for np.real_if_close.
    :return: Fidelity which is a scalar.
    """
    sqrt_rho = sqrtm_psd(rho)
    fid = (np.trace(sqrtm_psd(sqrt_rho @ sigma @ sqrt_rho))) ** 2
    return np.ndarray.item(np.real_if_close(fid, tol)) 

def impurity(rho: np.ndarray, dim_renorm=False, tol: float = 1000) -> float:
    """
    Calculates the impurity (or linear entropy) :math:`L = 1 - tr[ρ^2]` of a quantum state ρ.

    As stated above the lower value of the impurity depends on the dimension of ρ's Hilbert space.
    For some applications this can be undesirable. For this reason we introduce an optional
    dimensional renormalization flag with the following behavior

    If the dimensional renormalization flag is FALSE (default) then  0 ≤ L ≤ 1/dim.
    If the dimensional renormalization flag is TRUE then 0 ≤ L ≤ 1.

    where dim is the dimension of ρ's Hilbert space.

    :param rho: Is a dim by dim positive matrix with unit trace.
    :param dim_renorm: Boolean, default False.
    :param tol: Tolerance in machine epsilons for np.real_if_close.
    :return: L the impurity of the state.
    """
    imp = 1 - np.trace(rho @ rho)
    if dim_renorm:
        dim = rho.shape[0]
        imp = (dim / (dim - 1.0)) * imp
    return np.ndarray.item(np.real_if_close(imp, tol)) 

def partial_trace_super(dim1: int, dim2: int) -> np.array:
    """
    Return the partial trace superoperator in the column-major basis.

    This returns the superoperator S_TrB such that:
        S_TrB * vec(rho_AB) = vec(rho_A)
    for rho_AB = kron(rho_A, rho_B)

    Args:
        dim1: the dimension of the system not being traced
        dim2: the dimension of the system being traced over

    Returns:
        A Numpy array of the partial trace superoperator S_TrB.
    """

    iden = sps.identity(dim1)
    ptr = sps.csr_matrix((dim1 * dim1, dim1 * dim2 * dim1 * dim2))

    for j in range(dim2):
        v_j = sps.coo_matrix(([1], ([0], [j])), shape=(1, dim2))
        tmp = sps.kron(iden, v_j.tocsr())
        ptr += sps.kron(tmp, tmp)

    return ptr 

def purity(rho: np.ndarray, dim_renorm=False, tol: float = 1000) -> float:
    """
    Calculates the purity :math:`P = tr[ρ^2]` of a quantum state ρ.

    As stated above lower value of the purity depends on the dimension of ρ's Hilbert space. For
    some applications this can be undesirable. For this reason we introduce an optional dimensional
    renormalization flag with the following behavior

    If the dimensional renormalization flag is FALSE (default) then  1/dim ≤ P ≤ 1.
    If the dimensional renormalization flag is TRUE then 0 ≤ P ≤ 1.

    where dim is the dimension of ρ's Hilbert space.

    :param rho: Is a dim by dim positive matrix with unit trace.
    :param dim_renorm: Boolean, default False.
    :param tol: Tolerance in machine epsilons for np.real_if_close.
    :return: P the purity of the state.
    """
    p = np.trace(rho @ rho)
    if dim_renorm:
        dim = rho.shape[0]
        p = (dim / (dim - 1.0)) * (p - 1.0 / dim)
    return np.ndarray.item(np.real_if_close(p, tol)) 

def main(n=3, p=3):
    # Generate problem data
    C = randn_psd(n)
    As = [randn_symm(n) for _ in range(p)]
    Bs = np.random.randn(p)

    # Extract problem data using cvxpy
    X = cp.Variable((n, n), PSD=True)
    objective = cp.trace(C@X)
    constraints = [cp.trace(As[i]@X) == Bs[i] for i in range(p)]
    prob = cp.Problem(cp.Minimize(objective), constraints)
    A, b, c, cone_dims = scs_data_from_cvxpy_problem(prob)

    # Print problem size
    mn_plus_m_plus_n = A.size + b.size + c.size
    n_plus_2n = c.size + 2 * b.size
    entries_in_derivative = mn_plus_m_plus_n * n_plus_2n
    print(f"""n={n}, p={p}, A.shape={A.shape}, nnz in A={A.nnz}, derivative={mn_plus_m_plus_n}x{n_plus_2n} ({entries_in_derivative} entries)""")

    # Compute solution and derivative maps
    start = time.perf_counter()
    x, y, s, derivative, adjoint_derivative = diffcp.solve_and_derivative(
        A, b, c, cone_dims, eps=1e-5)
    end = time.perf_counter()
    print("Compute solution and set up derivative: %.2f s." % (end - start))

    # Derivative
    lsqr_args = dict(atol=1e-5, btol=1e-5)
    start = time.perf_counter()
    dA, db, dc = adjoint_derivative(diffcp.cones.vec_symm(
        C), np.zeros(y.size), np.zeros(s.size), **lsqr_args)
    end = time.perf_counter()
    print("Evaluate derivative: %.2f s." % (end - start))

    # Adjoint of derivative
    start = time.perf_counter()
    dx, dy, ds = derivative(A, b, c, **lsqr_args)
    end = time.perf_counter()
    print("Evaluate adjoint of derivative: %.2f s." % (end - start)) 

def _project(self, matrix):
        """The largest value in each row is set to 1.
        """
        result = np.zeros(self.shape)
        for i in range(self.shape[0]):
            idx = np.argmax(matrix[i,:])
            result[i, idx] = 1
        return result
        # ordering = self.a.T.dot(result)
        # indices = np.argsort(ordering)
        # return result[:,indices]
        # import cvxpy as cvx
        # X = cvx.Bool(self.shape)
        # constr = [cvx.sum(X, axis=1) == 1, cvx.diff((self.a.T*X).T) >= 0]
        # prob = cvx.Problem(cvx.Maximize(cvx.trace(matrix.T*X)), constr)
        # prob.solve(solver=cvx.GUROBI, timeLimit=10)
        # return X.value

    # def _neighbors(self, matrix):
    #     neighbors_list = []
    #     idxs = np.argmax(matrix, axis=1)
    #     for i in range(self.shape[0]):
    #         for j in range(self.shape[1]):
    #             if j != idxs[i]:
    #                 new_mat = matrix.copy()
    #                 new_mat[i,j] = 1
    #                 new_mat[i,idxs[i]] = 0
    #                 neighbors_list += [new_mat]
    #     return neighbors_list 

def trace_distance(rho: np.ndarray, sigma: np.ndarray) -> float:
    r"""
    Computes the trace distance between two states rho and sigma:

    .. math::

        T(\rho, \sigma) = (1/2)||\rho-\sigma||_1

    where :math:`||X||_1` denotes the 1 norm of X.

    :param rho: Is a dim by dim positive matrix with unit trace.
    :param sigma: Is a dim by dim positive matrix with unit trace.
    :return: Trace distance which is a scalar.
    """
    return (0.5) * np.linalg.norm(rho - sigma, 1) 

def infidelity(rho: np.ndarray, sigma: np.ndarray, tol: float = 1000) -> float:
    r"""
    Computes the infidelity, :math:`1 - F(\rho, \sigma)`, between two quantum states rho and sigma
    where :math:`F(\rho, \sigma)` is the fidelity.

    :param rho: Is a dim by dim positive matrix with unit trace.
    :param sigma: Is a dim by dim positive matrix with unit trace.
    :param tol: Tolerance in machine epsilons for np.real_if_close.
    :return: Infidelity which is a scalar.
    """
    return 1 - fidelity(rho, sigma, tol) 

def test_sdp(self):
        set_seed(2)

        n = 3
        p = 3
        C = cp.Parameter((n, n))
        A = [cp.Parameter((n, n)) for _ in range(p)]
        b = [cp.Parameter((1, 1)) for _ in range(p)]

        C_tch = torch.randn(n, n, requires_grad=True).double()
        A_tch = [torch.randn(n, n, requires_grad=True).double()
                 for _ in range(p)]
        b_tch = [torch.randn(1, 1, requires_grad=True).double()
                 for _ in range(p)]

        X = cp.Variable((n, n), symmetric=True)
        constraints = [X >> 0]
        constraints += [
            cp.trace(A[i]@X) == b[i] for i in range(p)
        ]
        prob = cp.Problem(cp.Minimize(cp.trace(C@X) + cp.sum_squares(X)),
                          constraints)
        layer = CvxpyLayer(prob, [C] + A + b, [X])
        torch.autograd.gradcheck(lambda *x: layer(*x,
                                                  solver_args={'eps': 1e-12}),
                                 [C_tch] + A_tch + b_tch,
                                 eps=1e-6,
                                 atol=1e-3,
                                 rtol=1e-3) 

def test_sdp(self):
        tf.random.set_seed(5)

        n = 3
        p = 3
        C = cp.Parameter((n, n))
        A = [cp.Parameter((n, n)) for _ in range(p)]
        b = [cp.Parameter((1, 1)) for _ in range(p)]

        C_tf = tf.Variable(tf.random.normal((n, n), dtype=tf.float64))
        A_tf = [tf.Variable(tf.random.normal((n, n), dtype=tf.float64))
                for _ in range(p)]
        b_tf = [tf.Variable(tf.random.normal((1, 1), dtype=tf.float64))
                for _ in range(p)]

        X = cp.Variable((n, n), symmetric=True)
        constraints = [X >> 0]
        constraints += [
            cp.trace(A[i]@X) == b[i] for i in range(p)
        ]
        problem = cp.Problem(cp.Minimize(
            cp.trace(C @ X) - cp.log_det(X) + cp.sum_squares(X)),
            constraints)
        layer = CvxpyLayer(problem, [C] + A + b, [X])
        values = [C_tf] + A_tf + b_tf
        with tf.GradientTape() as tape:
            soln = layer(*values,
                         solver_args={'eps': 1e-10, 'max_iters': 10000})[0]
            summed = tf.math.reduce_sum(soln)
        grads = tape.gradient(summed, values)

        def f():
            problem.solve(cp.SCS, eps=1e-10, max_iters=10000)
            return np.sum(X.value)

        numgrads = numerical_grad(f, [C] + A + b, values, delta=1e-4)
        for g, ng in zip(grads, numgrads):
            np.testing.assert_allclose(g, ng, atol=1e-1) 

def gate_error(channel, target=None, require_cp=True, require_tp=False):
    r"""Return the gate error of a noisy quantum channel.

    The gate error :math:`E` is given by the average gate infidelity

    .. math::
        E(\mathcal{E}, U) = 1 - F_{\text{ave}}(\mathcal{E}, U)

    where :math:`F_{\text{ave}}(\mathcal{E}, U)` is the
    :meth:`~qiskit.quantum_info.average_gate_fidelity` of the input
    quantum *channel* :math:`\mathcal{E}` with a *target* unitary
    :math:`U`.

    Args:
        channel (QuantumChannel): noisy quantum channel.
        target (Operator or None): target unitary operator.
            If `None` target is the identity operator [Default: None].
        require_cp (bool): require channel to be completely-positive
            [Default: True].
        require_tp (bool): require channel to be trace-preserving
            [Default: False].

    Returns:
        float: The average gate error :math:`E`.

    Raises:
        QiskitError: if the channel and target do not have the same dimensions,
                     or have different input and output dimensions.
        QiskitError: if the channel and target or are not completely-positive
                     (with ``require_cp=True``) or not trace-preserving
                     (with ``require_tp=True``).
    """
    return 1 - average_gate_fidelity(
        channel, target=target, require_cp=require_cp, require_tp=require_tp) 

def entanglement_fidelity(A, B, mxBasis='pp'):
    """
    Returns the "entanglement" process fidelity between gate
    matrices A and B given by :

      F = Tr( sqrt{ sqrt(J(A)) * J(B) * sqrt(J(A)) } )^2

    where J(.) is the Jamiolkowski isomorphism map that maps a operation matrix
    to it's corresponding Choi Matrix.

    Parameters
    ----------
    A, B : numpy array
        The matrices to compute the fidelity between.

    mxBasis : {'std', 'gm', 'pp', 'qt'} or Basis object
        The basis of the matrices.  Allowed values are Matrix-unit (std),
        Gell-Mann (gm), Pauli-product (pp), and Qutrit (qt)
        (or a custom basis object).
    """
    d2 = A.shape[0]
    def isTP(x): return _np.isclose(x[0, 0], 1.0) and all(
        [_np.isclose(x[0, i], 0) for i in range(d2)])

    def isUnitary(x): return _np.allclose(_np.identity(d2, 'd'), _np.dot(x, x.conjugate().T))

    if isTP(A) and isTP(B) and isUnitary(B):  # then assume TP-like gates & use simpler formula
        TrLambda = _np.trace(_np.dot(A, B.conjugate().T))  # same as using _np.linalg.inv(B)
        d2 = A.shape[0]
        return TrLambda / d2

    JA = _jam.jamiolkowski_iso(A, mxBasis, mxBasis)
    JB = _jam.jamiolkowski_iso(B, mxBasis, mxBasis)
    return fidelity(JA, JB) 

def tracedist(A, B):
    """
    Compute the trace distance between matrices A and B,
    given by:

      D = 0.5 * Tr( sqrt{ (A-B)^dagger * (A-B) } )

    Parameters
    ----------
    A, B : numpy array
        The matrices to compute the distance between.
    """
    return 0.5 * tracenorm(A - B) 

def average_gate_fidelity(channel,
                          target=None,
                          require_cp=True,
                          require_tp=False):
    r"""Return the average gate fidelity of a noisy quantum channel.

    The average gate fidelity :math:`F_{\text{ave}}` is given by

    .. math::
        F_{\text{ave}}(\mathcal{E}, U)
            &= \int d\psi \langle\psi|U^\dagger
                \mathcal{E}(|\psi\rangle\!\langle\psi|)U|\psi\rangle \\
            &= \frac{d F_{\text{pro}}(\mathcal{E}, U) + 1}{d + 1}

    where :math:`F_{\text{pro}}(\mathcal{E}, U)` is the
    :meth:`~qiskit.quantum_info.process_fidelity` of the input quantum
    *channel* :math:`\mathcal{E}` with a *target* unitary :math:`U`, and
    :math:`d` is the dimension of the *channel*.

    Args:
        channel (QuantumChannel): noisy quantum channel.
        target (Operator or None): target unitary operator.
            If `None` target is the identity operator [Default: None].
        require_cp (bool): require channel to be completely-positive
            [Default: True].
        require_tp (bool): require channel to be trace-preserving
            [Default: False].

    Returns:
        float: The average gate fidelity :math:`F_{\text{ave}}`.

    Raises:
        QiskitError: if the channel and target do not have the same dimensions,
                     or have different input and output dimensions.
        QiskitError: if the channel and target or are not completely-positive
                     (with ``require_cp=True``) or not trace-preserving
                     (with ``require_tp=True``).
    """
    if isinstance(channel, (list, np.ndarray, Operator, Pauli)):
        channel = Operator(channel)
    else:
        channel = SuperOp(channel)
    dim, _ = channel.dim
    f_pro = process_fidelity(channel,
                             target=target,
                             require_cp=require_cp,
                             require_tp=require_tp)
    return (dim * f_pro + 1) / (dim + 1) 

def fidelity_upper_bound(operationMx):
    """
    Get an upper bound on the fidelity of the given
      operation matrix with any unitary operation matrix.

    The closeness of the result to one tells
     how "unitary" the action of operationMx is.

    Parameters
    ----------
    operationMx : numpy array
        The operation matrix to act on.

    Returns
    -------
    float
        The resulting upper bound on fidelity(operationMx, anyUnitaryGateMx)
    """
    choi = _jam.jamiolkowski_iso(operationMx, choiMxBasis="std")
    choi_evals, choi_evecs = _np.linalg.eig(choi)
    maxF_direct = max([_np.sqrt(max(ev.real, 0.0)) for ev in choi_evals]) ** 2

    iMax = _np.argmax([ev.real for ev in choi_evals])  # index of maximum eigenval
    closestVec = choi_evecs[:, iMax:(iMax + 1)]

    # #print "DEBUG: closest evec = ", closestUnitaryVec
    # new_evals = _np.zeros( len(closestUnitaryVec) ); new_evals[iClosestU] = 1.0
    # # gives same result:
    # closestUnitaryJmx = _np.dot(choi_evecs, _np.dot( _np.diag(new_evals), _np.linalg.inv(choi_evecs) ) )
    closestJmx = _np.kron(closestVec, _np.transpose(_np.conjugate(closestVec)))  # closest rank-1 Jmx
    closestJmx /= _mt.trace(closestJmx)  # normalize so trace of Jmx == 1.0

    maxF = fidelity(choi, closestJmx)

    if not _np.isnan(maxF):

        #Uncomment for debugging
        #if abs(maxF - maxF_direct) >= 1e-6:
        #    print "DEBUG: operationMx:\n",operationMx
        #    print "DEBUG: choiMx:\n",choi
        #    print "DEBUG choi_evals = ",choi_evals, " iMax = ",iMax
        #    #print "DEBUG: J = \n", closestUnitaryJmx
        #    print "DEBUG: eigvals(J) = ", _np.linalg.eigvals(closestJmx)
        #    print "DEBUG: trace(J) = ", _mt.trace(closestJmx)
        #    print "DEBUG: maxF = %f,  maxF_direct = %f" % (maxF, maxF_direct)
        #    raise ValueError("ERROR: maxF - maxF_direct = %f" % (maxF -maxF_direct))
        assert(abs(maxF - maxF_direct) < 1e-6)
    else:
        maxF = maxF_direct  # case when maxF is nan, due to scipy sqrtm function being buggy - just use direct F

    closestOpMx = _jam.jamiolkowski_iso_inv(closestJmx, choiMxBasis="std")
    return maxF, closestOpMx

    #closestU_evals, closestU_evecs = _np.linalg.eig(closestUnitaryGateMx)
    #print "DEBUG: U = \n", closestUnitaryGateMx
    #print "DEBUG: closest U evals = ",closestU_evals
    #print "DEBUG:  evecs = \n",closestU_evecs 

def quantum_chernoff_bound(rho: np.ndarray,
                           sigma: np.ndarray,
                           tol: float = 1000) -> Tuple[float, float]:
    r"""
    Computes the quantum Chernoff bound between rho and sigma.

    It is defined as

    .. math::

        ξ_{QCB}(\rho, \sigma) = - \log[ \min_{0\le s\le 1} tr(\rho^s \sigma^{1-s}) ]

    It is also common to study the non-logarithmic variety of the quantum Chernoff bound

    .. math::

        Q_{QCB}(\rho, \sigma) = \min_{0\le s\le 1} tr(\rho^s \sigma^{1-s})

    The quantum Chernoff bound has many nice properties, see [QCB]_. Importantly it is
    operationally important in the following context. Given n copies of rho or sigma the minimum
    error probability for discriminating rho from sigma is :math:`P_{e,min,n} ~ exp[-n ξ_{QCB}]`.

    .. [QCB] The Quantum Chernoff Bound.
          Audenaert et al.
          Phys. Rev. Lett. 98, 160501 (2007).
          https://dx.doi.org/10.1103/PhysRevLett.98.160501
          https://arxiv.org/abs/quant-ph/0610027

    :param rho: Is a dim by dim positive matrix with unit trace.
    :param sigma: Is a dim by dim positive matrix with unit trace.
    :param tol: Tolerance in machine epsilons for np.real_if_close.
    :return: The non-logarithmic quantum Chernoff bound and the s achieving the minimum.
    """

    def f(s):
        s = np.real_if_close(s)
        return np.trace(
            np.matmul(fractional_matrix_power(rho, s), fractional_matrix_power(sigma, 1 - s)))

    f_min = minimize_scalar(f, bounds=(0, 1), method='bounded')
    s_opt = np.real_if_close(f_min.x, tol)
    qcb = np.real_if_close(f_min.fun, tol)
    return qcb, s_opt 

def entanglement_fidelity(pauli_lio0: np.ndarray,
                          pauli_lio1: np.ndarray,
                          tol: float = 1000) -> float:
    r"""
    Returns the entanglement fidelity (F_e) between two channels, E and F, represented as Pauli
    Liouville matrix.

    The expression is

    .. math::

            F_e(E,F) = Tr[E^\dagger F] / (dim^2),

    where dim is the dimension of the Hilbert space associated with E and F.

    See the following references for more information:

    [GRAPTN]_ referenced in the superoperator_tools module. In particular section V subsection G.

    .. [H**3] General teleportation channel, singlet fraction and quasi-distillation.
           Horodecki et al.
           PRA 60, 1888 (1999).
           https://doi.org/10.1103/PhysRevA.60.1888
           https://arxiv.org/abs/quant-ph/9807091

    .. [GFID] A simple formula for the average gate fidelity of a quantum dynamical operation.
           M. Nielsen.
           Physics Letters A 303, 249 (2002).
           https://doi.org/10.1016/S0375-9601(02)01272-0
           https://arxiv.org/abs/quant-ph/0205035

    :param pauli_lio0: A dim**2 by dim**2 Pauli-Liouville matrix
    :param pauli_lio1: A dim**2 by dim**2 Pauli-Liouville matrix
    :param tol: Tolerance in machine epsilons for np.real_if_close.
    :return: Returns the entanglement fidelity between pauli_lio0 and pauli_lio1 which is a scalar.
    """
    assert pauli_lio0.shape == pauli_lio1.shape
    assert pauli_lio0.shape[0] == pauli_lio1.shape[1]
    dim_squared = pauli_lio0.shape[0]
    dim = int(np.sqrt(dim_squared))
    Fe = np.trace(np.matmul(np.transpose(np.conj(pauli_lio0)), pauli_lio1)) / (dim ** 2)
    return np.ndarray.item(np.real_if_close(Fe, tol)) 

def fidelity(A, B):
    """
    Returns the quantum state fidelity between density
      matrices A and B given by :

      F = Tr( sqrt{ sqrt(A) * B * sqrt(A) } )^2

    To compute process fidelity, pass this function the
    Choi matrices of the two processes, or just call
    :function:`entanglement_fidelity` with the operation matrices.

    Parameters
    ----------
    A : numpy array
        First density matrix.

    B : numpy array
        Second density matrix.

    Returns
    -------
    float
        The resulting fidelity.
    """
    evals, U = _np.linalg.eig(A)
    if len([ev for ev in evals if abs(ev) > 1e-8]) == 1:
        # special case when A is rank 1, A = vec * vec^T and sqrt(A) = A
        ivec = _np.argmax(evals)
        vec = U[:, ivec:(ivec + 1)]
        F = evals[ivec].real * _np.dot(_np.conjugate(_np.transpose(vec)), _np.dot(B, vec)).real  # vec^T * B * vec
        return float(F)

    evals, U = _np.linalg.eig(B)
    if len([ev for ev in evals if abs(ev) > 1e-8]) == 1:
        # special case when B is rank 1 (recally fidelity is sym in args)
        ivec = _np.argmax(evals)
        vec = U[:, ivec:(ivec + 1)]
        F = evals[ivec].real * _np.dot(_np.conjugate(_np.transpose(vec)), _np.dot(A, vec)).real  # vec^T * A * vec
        return float(F)

    #if _np.array_equal(A, B): return 1.0  # HACK - some cases when A and B are perfecty equal sqrtm(A) fails...
    sqrtA = _hack_sqrtm(A)  # _spl.sqrtm(A)
    # test the scipy sqrtm function - sometimes fails when rank defficient
    #assert(_np.linalg.norm(_np.dot(sqrtA, sqrtA) - A) < 1e-8)
    if _np.linalg.norm(_np.dot(sqrtA, sqrtA) - A) > 1e-8:
        evals = _np.linalg.eigvals(A)
        _warnings.warn(("sqrtm(A) failure when computing fidelity - beware result. "
                        "Maybe due to rank defficiency - eigenvalues of A are: %s") % evals)
    F = (_mt.trace(_hack_sqrtm(_np.dot(sqrtA, _np.dot(B, sqrtA)))).real)**2  # Tr( sqrt{ sqrt(A) * B * sqrt(A) } )^2
    return float(F)