OperatorSinkhorn

From QETLAB
Jump to: navigation, search
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 [1]).

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[2] 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 on "expand" to the right to view the MATLAB source code for this function.

  1. %%  OPERATORSINKHORN    Performs the operator Sinkhorn iteration
  2. %   This function has one required argument:
  3. %     RHO: a density matrix
  4. %
  5. %   SIGMA = OperatorSinkhorn(RHO) is a density matrix that is locally
  6. %   equivalent to RHO, but has both of its partial traces proportional to
  7. %   the identity (RHO must be bipartite; if it is multipartite, see the
  8. %   optional arguments below).
  9. %
  10. %   Such a density matrix SIGMA does not always exist if RHO is low-rank.
  11. %   An error is returned in these cases.
  12. %
  13. %   This function has two optional input arguments:
  14. %     DIM (default has two subsystems of equal dimension)
  15. %     TOL (default sqrt(eps))
  16. %
  17. %   This function has one optional output argument:
  18. %     F: a cell containing local matrices
  19. %
  20. %   [SIGMA,F] = OperatorSinkhorn(RHO,DIM,TOL) returns SIGMA and F such that
  21. %   SIGMA has all of its (single-party) reduced density matrices
  22. %   proportional to the identity, and SIGMA = Tensor(F)*RHO*Tensor(F)'. In
  23. %   other words, F contains invertible local operations that demonstrate
  24. %   that RHO and SIGMA are locally equivalent.
  25. %
  26. %   DIM is a 1-by-2 vector containing the dimensions of the subsystems on
  27. %   which RHO acts (RHO can act on any number of parties). TOL is the
  28. %   numerical tolerance used when determining when the operator Sinkhorn
  29. %   iteration has converged.
  30. %
  31. %   URL: http://www.qetlab.com/OperatorSinkhorn
  32.  
  33. %   requires: opt_args.m, PartialTrace.m, PermuteSystems.m
  34. %             
  35. %   author: Nathaniel Johnston (nathaniel@njohnston.ca)
  36. %   package: QETLAB
  37. %   last updated: October 3, 2014
  38.  
  39. function [sigma,F] = OperatorSinkhorn(rho,varargin)
  40.  
  41. dX = length(rho);
  42. sdX = round(sqrt(dX));
  43. tr_rho = trace(rho);
  44.  
  45. % set optional argument defaults: dim=sqrt(length(rho)), tol=sqrt(eps)
  46. [dim,tol] = opt_args({ [sdX, sdX], sqrt(eps) },varargin{:});
  47. num_sys = length(dim);
  48.  
  49. % allow the user to enter a single number for dim
  50. if(num_sys == 1)
  51.     dim = [dim,dX/dim];
  52.     if abs(dim(2) - round(dim(2))) >= 2*dX*eps
  53.         error('OperatorSinkhorn:InvalidDim','If DIM is a scalar, X must be square and DIM must evenly divide length(X); please provide the DIM array containing the dimensions of the subsystems.');
  54.     end
  55.     dim(2) = round(dim(2));
  56.     num_sys = 2;
  57. end
  58. tr_rho_p = tr_rho^(1/(2*num_sys));
  59.  
  60. % Prepare the iteration.
  61. for j = num_sys:-1:1
  62.     Prho{j} = eye(dim(j))/dim(j);
  63.     Prho_tmp{j} = Prho{j};
  64.     F{j} = eye(dim(j))*tr_rho_p;
  65.     ldim(j) = prod(dim(1:j-1));
  66.     rdim(j) = prod(dim(j+1:end));
  67. end
  68.  
  69. % Perform the operator Sinkhorn iteration.
  70. lastwarn(''); % clears any previous warnings
  71. warning('off','MATLAB:singularMatrix'); % we want to catch invertibility warnings, but not display them
  72. warning('off','MATLAB:nearlySingularMatrix'); % we want to catch invertibility warnings, but not display them
  73. it_err = 1;
  74. while it_err > tol
  75.     it_err = 0;
  76.     max_cond = 0;
  77.  
  78.     % Loop over each of the systems and apply a filter on each one.
  79.     try
  80.         for j = 1:num_sys
  81.             % Compute the reduced density matrix on the j-th system.
  82.             Prho_tmp{j} = PartialTrace(rho,setdiff(1:num_sys,j),dim);
  83.             Prho_tmp{j} = (Prho_tmp{j}+Prho_tmp{j}')/2; % for numerical stability
  84.             it_err = it_err + norm(Prho{j}-Prho_tmp{j});
  85.             Prho{j} = Prho_tmp{j};
  86.  
  87.             % Apply the filter.
  88.             T = sqrtm(inv(Prho{j}))/sqrt(dim(j));
  89.             Tk = kron(speye(ldim(j)),kron(T,speye(rdim(j))));
  90.             rho = Tk*rho*Tk';
  91.             F{j} = T*F{j};
  92.  
  93.             max_cond = max(max_cond,cond(F{j}));
  94.         end
  95.     catch err
  96.         it_err = 1;
  97.     end
  98.  
  99.     % Make sure that the local transformation performed is invertible --
  100.     % otherwise, the iteration will typically not converge and the results
  101.     % will be meaningless.
  102.     [~,warnid] = lastwarn;
  103.     if(it_err == 1 || max_cond >= 1/tol || strcmpi('MATLAB:nearlySingularMatrix',warnid) || strcmpi('MATLAB:singularMatrix',warnid))
  104.         error('OperatorSinkhorn:LowRank','The operator Sinkhorn iteration does not converge for RHO. This is often the case if RHO is not of full rank.');
  105.     end
  106. end
  107. warning('on','MATLAB:singularMatrix')
  108. warning('on','MATLAB:nearlySingularMatrix')
  109.  
  110. sigma = (rho+rho')/2; % done for numerical stability reasons
  111. sigma = tr_rho*sigma/trace(sigma); % correct the scaling of the output

References

  1. L. Gurvits. Classical complexity and quantum entanglement. J. Comput. System Sci., 69:448–484, 2004.
  2. J. M. Leinaas, J. Myrheim, and E. Ovrum. Geometrical aspects of entanglement. Phys. Rev. A, 74:012313, 2006. E-print: arXiv:quant-ph/0605079