OperatorSinkhorn

OperatorSinkhorn
Performs the operator Sinkhorn iteration, making all single-party reduced states proportional to the identity
Other toolboxes required	none
Related functions	FilterNormalForm OperatorSchmidtDecomposition
Function category	Miscellaneous

OperatorSinkhorn is a function that performs operator Sinkhorn iterative scaling on a (potentially multipartite) operator, resulting in a locally equivalent operator that has all of its single-party reduced states proportional to the identity (see Section 4 of ).

Syntax

• SIGMA = OperatorSinkhorn(RHO)
• SIGMA = OperatorSinkhorn(RHO,DIM)
• SIGMA = OperatorSinkhorn(RHO,DIM,TOL)
• [SIGMA,F] = OperatorSinkhorn(RHO,DIM,TOL)

Argument descriptions

Input arguments

• RHO: An positive semidefinite operator (typically a density matrix) that acts on a bipartite or multipartite Hilbert space. This operator will have the operator Sinkhorn iteration applied to it.
• DIM (optional, default has RHO acting on two subsystems of equal size): A vector containing the dimensions of the two or more subsystems on which RHO acts.
• TOL (optional, default sqrt(eps)): The numerical tolerance used to determine when the iteration has converged.

Output arguments

• SIGMA: A positive semidefinite operator with three properties: (1) its trace is the same trace as RHO, (2) it is locally equivalent to RHO (i.e., there exist invertible local operators that transform RHO into SIGMA), and (3) all of its single-party reduced density matrices are proportional to the identity operator.
• F (optional): A cell containing the local filtering operations used to convert RHO into SIGMA. That is, a cell containing invertible operators so that Tensor(F)*RHO*Tensor(F)' equals SIGMA.

Examples

Verifying the partial traces

The following code generates a random density matrix, performs operator Sinkhorn iterative scaling, and then verifies that the result has the claimed properties (both partial traces are proportional to the identity and it is locally equivalent to the input).

>> rho = RandomDensityMatrix(9);
>> [sigma,F] = OperatorSinkhorn(rho);
>> PartialTrace(sigma,1)

ans =

   0.3333 + 0.0000i   0.0000 - 0.0000i  -0.0000 - 0.0000i
   0.0000 + 0.0000i   0.3333 + 0.0000i   0.0000 - 0.0000i
  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.3333 + 0.0000i

>> PartialTrace(sigma,2)

ans =

   0.3333 + 0.0000i  -0.0000 + 0.0000i   0.0000 - 0.0000i
  -0.0000 - 0.0000i   0.3333 + 0.0000i  -0.0000 + 0.0000i
   0.0000 + 0.0000i  -0.0000 - 0.0000i   0.3333 + 0.0000i

>> norm(Tensor(F)*rho*Tensor(F)' - sigma)

ans =

   1.8731e-16

A multipartite example

The following example is similar to the previous one, except the density matrix that is used lives on 3-qubit space (i.e., the tripartite space where each local dimension is 2).

>> rho = RandomDensityMatrix(8);
>> [sigma,F] = OperatorSinkhorn(rho,[2,2,2]);
>> PartialTrace(sigma,[1,2],[2,2,2])

ans =

   0.5000 + 0.0000i  -0.0000 - 0.0000i
  -0.0000 + 0.0000i   0.5000 + 0.0000i

>> PartialTrace(sigma,[1,3],[2,2,2])

ans =

   0.5000 + 0.0000i  -0.0000 + 0.0000i
  -0.0000 - 0.0000i   0.5000 + 0.0000i

>> PartialTrace(sigma,[2,3],[2,2,2])

ans =

   0.5000 + 0.0000i  -0.0000 + 0.0000i
  -0.0000 - 0.0000i   0.5000 + 0.0000i

>> norm(Tensor(F)*rho*Tensor(F)' - sigma)

ans =

   3.4349e-16

Iteration may not converge for low-rank states

It is known that the operator Sinkhorn iteration converges whenever RHO has full rank. However, the iteration may fail to converge for low-rank density matrices – an error will be produced by this script in these cases. The following code demonstrates an example.

>> rho = kron(RandomDensityMatrix(3,0,2),eye(3)); % clearly can not have its left partial trace scaled to be proportional to the identity
>> OperatorSinkhorn(rho)
Error using OperatorSinkhorn (line 100)
The operator Sinkhorn iteration does not converge for RHO. This is often the case if RHO is not of full
rank.

Source code

Click here to view this function's source code on github.