Using entropy_mutual function in QuTiP

Question

I am trying to calculate mutual entropies using QuTiP, but I am being unsuccessful so far. More specifically, I consider a 2^n x 2^n matrix representing the density operator of a n-qubit bipartite system AB made of system A (first m < n qubits) and B (remaining n-m qubits). No tutorial nor material on the internet addressed this specific task.

For simplicity, let us consider a 1-qubit system A and a 2-qubit system B and a density operator of dimension 8x8 representing AB in computational basis.

More practically in python, let rhoAB = Qobj=(np.random.rand(8,8)), and assume that this is a valid density operator.

How should I call entropy_mutual so that I can get this measure between A and B, in particular, regarding the arguments selA and selB? Ideally, I would call something like entopy_mutual(rhoAB, selA=[1], selB=[2,3]) but this not the approach how the function interprets the subsystems and their respective dimensions.

Answer 1

Best to look at the source code when the documentation isn't helpful enough. The definition of entropy_mutual is

def entropy_mutual(rho, selA, selB, base=e, sparse=False):
    """
    Calculates the mutual information S(A:B) between selection
    components of a system density matrix.

    Parameters
    ----------
    rho : qobj
        Density matrix for composite quantum systems
    selA : int/list
        `int` or `list` of first selected density matrix components.
    selB : int/list
        `int` or `list` of second selected density matrix components.
    base : {e,2}
        Base of logarithm.
    sparse : {False,True}
        Use sparse eigensolver.

    Returns
    -------
    ent_mut : float
       Mutual information between selected components.

    """
    if isinstance(selA, int):
        selA = [selA]
    if isinstance(selB, int):
        selB = [selB]
    if rho.type != 'oper':
        raise TypeError("Input must be a density matrix.")
    if (len(selA) + len(selB)) != len(rho.dims[0]):
        raise TypeError("Number of selected components must match " +
                        "total number.")

    rhoA = ptrace(rho, selA)
    rhoB = ptrace(rho, selB)
    out = (entropy_vn(rhoA, base, sparse=sparse) +
           entropy_vn(rhoB, base, sparse=sparse) -
           entropy_vn(rho, base, sparse=sparse))
    return out

So we see selA and selB are passed as arguments to compute the partial trace. I am not too familiar with qutip but here is an example computing $S(A:B)$ for $\rho_{AB}$ where $A$ is a qubit system and $B$ is a two-qubit system.

import qutip as qtp
# note there is a rand_dm function
# We should also let qutip know how are systems are partitioned
# This is so it knows how to correctly compute the partial trace
rho = qtp.rand_dm(8, dims=[[2,4],[2,4]])
qtp.entropy_mutual(rho,0,1)

With the above example we could also specify the second system as two-qubits instead of a four dimensional system i.e.

rho = qtp.rand_dm(8, dims=[[2,2,2],[2,2,2]])
qtp.entropy_mutual(rho,0,[1,2])