DEPARTED provides partial trace and partial transpose operations for
matrices with a Kronecker-product structure. 

Installation

  As of release 0.3.0, it is possible to retrieve the latest release of the
  package from the PyPi repository by running 'pip install --user departed'

Examples

  See the examples directory.

Documentation

  - ptrace (matrix, component_dims, component_mask)

    Computes the partial trace of a matrix with a Kronecker-product structure.
    
    : matrix 
      defines the matrix to be partially traced

    : component_dims
      defines a list of component dimensions

    : component_mask
      defines a mask specifying whether a component should be 
      traced out (1, True) or kept (0, False).

    Example

      Consider a density matrix R of a three-qubit quantum state. The partial
      trace over the 1st and 3rd qubit is obtained with 

        ptrace(R, [ 2, 2, 2 ], [ 1, 0, 1 ]),

      where the mask [ 1, 0, 1 ] uses 1 for components that will be traced out.

  - ptranspose (matrix, component_dims, component_mask)

    Computes the partial transpose of a matrix with a Kronecker-product structure.

    : matrix 
      defines the matrix to be partially transposed

    : component_dims
      defines a list of component dimensions

    : component_mask
      defines a mask specifying whether a component should be 
      transposed (1, True) or left intact (0, False).

    Example

      Consider a density matrix R of a three-qubit quantum state. The partial
      transpose over the 1st and 3rd qubit is obtained with

        ptranspose(R, [ 2, 2, 2 ], [ 1, 0, 1 ]),

      where the mask [ 1, 0, 1 ] uses 1 for components that will be transposed.

  - mask_from_component_list (component_list, mask_width, invert = False)

    Constructs a component_mask from a list of component indices.

    By default, 1 is set for components listed in index_list and 0 for those
    unlisted. These values can be inverted by setting invert = True. 

    : component_list
      defines a list of component indices (starting from 0)

    : mask_width
      defines the total number of components (mask width)

    : invert
      defines how component_list is translated into component_mask

    Example

      Consider a system with 5 components. Suppose the partial operation of
      choice is performed over the 2nd and 4th components. The
      corresponding mask [ 0, 1, 0, 1, 0 ] can be constructed with

        mask_from_component_list([ 1, 3 ], 5),

      where we note that the indices start from 0.

      Suppose that the partial operation should NOT affect the 3rd and 5th
      component. The respective mask [ 1, 1, 0, 1, 0 ] can be constructed with

        mask_from_component_list([ 2, 4 ], 5, invert = True),

      where setting the optional parameter (invert = True) inverts the mask.

Operating principles of Kronecker product and its tensor representation

  An understanding of the internal memory layout of matrices with a
  Kronecker-product structure is crucial for both the partial trace and partial
  transpose operations. In both algorithms, the matrix is first reshaped into a
  tensor with the axes of its constituent matrices arranged according to their
  direction, with row axes preceding column axes.

  In particular, consider a matrix with a Kronecker-product structure
  comprising N matrices and its tensor representation. The first N axes of the
  tensor are the row axes of its matrix constituents and the last N axes are
  their column axes.

  This ordering appears counter-intuitive, however, it follows from the
  internal memory layout of the underlying numpy.ndarray object.

  Example
  
    To illustrate this concept in greater detail, we consider a tensor product
    of three matrix spaces, where the first matrix space comprises (m x m)
    matrices, the second one (n x n) matrices, and finally, the third
    one (o x o) matrices. The product space contains (mno x mno) matrices.

    Let Q be a matrix (with a Kronecker-product structure) from this product
    space. Its tensor representation can be obtained through

                      T = Q.reshape([ m, n, o, m, n, o ]),

    where the first three axes of the tensor T each correspond to the row axes
    of matrices from the constituent matrix spaces and the last three axes of
    the tensor each correspond to the column axes of these matrices. 

    For example, should one have three matrices - A, B, and C - and their 
    Kronecker product, kron(kron(A, B), C), the element 

                              T[:, 0, 2, :, 1, 3 ]

    of its tensor representation equals to

                               A * B[0, 1], * C[2, 3].
  Remarks
  
    Some software libraries use tensors structurally compatible with
    Kronecker-product matrices to represent quantum states and operations.

Operating principles of partial transpose

  Consider an arbitrary bipartite quantum state defined by its density operator

               rho = sum(i, j, k, l) R(i, k, j, l) |i><j| |k><l|,

  where the rank four tensor R(i, k, j, l) determines its coefficients. The
  order of its indices (axes) follows the tensor representation of matrices
  with a Kronecker-product structure. Its row indices precede column indices.

  We define the partial transpose over the second system as

             pt(rho) = sum(i, j, k, l) R(i, k, j, l) |i><j| |l><k|
                     = sum(i, j, k, l) R(i, l, j, k) |i><j| |k><l|.

  If the density operator is encoded in a matrix with a Kronecker-product
  structure, rather than a rank four tensor, the partial transpose no longer
  amounts to a simple exchange of axes.

  It can be realized by determining the affected indices of the matrix
  and exchanging their values. Not only is this approach cumbersome, but also
  extremely inelegant.

  A better alternative is to reshape the matrix into a tensor of appropriate
  rank. The reshaping operation is computationally cheap as it only affects the
  interpretation of the underlying numpy.ndarray object.

  Suppose Q is the matrix representation of the state rho. Suppose m and n are
  the dimensions of the individual systems comprising the state. With

                         T = reshape(Q, [ m, n, m, n ])

  we obtain a rank four tensor T. Consequently, the partial transpose
  corresponds to the exchange of the 2nd and the 4th index,

                       pt(T)(i, k, j, l) = T(i, l, j, k).

  We can reconstruct the Kronecker representation by reshaping the tensor back
  into the original matrix shape.

  This idea can be generalized to an arbitrary number of components. It is
  practical to use the standard numpy.transpose procedure and supply the
  correct axial permutation. Our ptranspose does exactly that.

Operating principles of partial trace
  
  We build on the same ideas we have explored for partial transpose and then
  add a figurative cherry on top: by reordering the tensor, we can perform only
  a single trace operation even if multiple components are to be traced out.

  Let us begin with a bipartite system again and a density matrix

               rho = sum(i, j, k, l) R(i, k, j, l) |i><j| |k><l|.

  We define the partial trace performed over its second component as

         rho' = sum(u) sum(i, j, k, l) R(i, k, j, l) |i><j| <u|k> <l|u>
              = sum(i, j, k) R(i, k, j, k) |i><j|
              = sum(i, j) [ sum(k) R(i, k, j, k) ] |i><j|,

  where we have presumed that the individual { |k>, |l> } kets are orthonormal.

  If the density operator is encoded in a matrix Q with a Kronecker-product
  structure, the partial trace can be obtained from its tensor representation

                         T = reshape(Q, [ m, n, m, n ]).

  To obtain the partial trace, we must evaluate

                         M(i, j) = sum(k) T(i, k, j, k),

  which defines the elements M(i, j) of the marginal tensor. This approach can
  be extended for multi-partite systems and further optimized. We can permute
  the axes to gather the parts to be traced out into a single block and then
  compute the trace in one go.

Notes and acknowledgments

  This project was inspired by QuTiP (https://github.com/qutip/qutip) and their
  essentially identical implementation. This project aims to explain the
  underlying principles behind these functions.