Code accompanying the paper arXiv:1906.01648.
Download Matlab files:
The computation of the maximal achievable fidelity of distillation as well as the one-shot $\varepsilon$-error distillable entanglement under several classes of operations can be cast as semidefinite programs. These include $PPT_+$-preserving operations, Rains-preserving operations, and PPT operations; for pure states, the characterisaton includes even LOCC and one-way LOCC (see the manuscript for more details).
The implementation below requires CVX and QETLAB. The code takes a simplified approach when it encounters a pure state (either provided as a density matrix or state vector), and employs analytical (for the fidelity) or second-order conic programming (for the distillable entanglement) methods. For general states, it evaluates a semidefinite program. The notation is consistent with the paper where relevant.
Tested on Matlab R2018b with CVX 2.1 and QETLAB 0.9.
[fid, W] = dist_fid(rho, m, op, dim)
Computes the fidelity of entanglement distillation under a given class of operations
Required arguments:
rho: $n \times n$ density matrix or $n \times 1$ state vectorm: desired output size of maximally entangled stateOptional arguments:
op: class of operations (default ppt)
Allowed choices:
ppt (PPT operations),
pptp ($PPT_+$-preserving),
rains (Rains-preserving)
Disregarded when input is pure (all classes including LOCC give the same result)
dim: a vector [dA, dB] of the dimensions of the subsystems (default both systems of equal dimension)
Output:
fid: optimal fidelity of distillationW: optimal witness such that fid = trace(rho*W)
If the input is provided as a vector, W is also returned as one.
function [fid, W] = dist_fid(rho, m, op, dim)
dimrho = size(rho);
drho = max(dimrho);
if nargin < 4 % assume dimensions of subsystems are equal
dim = [sqrt(drho),sqrt(drho)];
d = sqrt(drho);
else
d = min(dim);
end
frho = full(rho);
if rank(frho) == 1 % pure state case
% obtain Schmidt decomposition
if min(dimrho) ~= 1
[~,~,v] = svd(frho);
rho = v(:,1);
end
[xi,U,V] = SchmidtDecomposition(rho, dim, -1);
% compute m-distillation norm of Schmidt vector
xi = sort(xi,'descend');
k = 1;
minsum = sum(xi(m:d).^2);
for j=2:m
jsum = sum(xi(m-j+1:d).^2)/j;
if jsum < minsum
minsum = jsum;
k = j;
end
end
fid = sum(xi(1:m-k))+k*sqrt(minsum);
fid = fid^2/m;
if nargout > 1
% construct optimal witness
wxi = xi./sqrt(minsum);
[~,wp] = maxk(xi,m-k);
wxi(wp) = 1;
W = zeros(drho,1);
for j=1:d
W = W + wxi(j) * kron(U(:,j), V(:,j));
end
W = W/sqrt(m);
if min(dimrho) ~= 1
W = W*W';
end
end
else % mixed state case
I = eye(drho);
rho = (rho+rho')/2; % avoids numerical problems
cvx_begin sdp quiet
variable W(drho,drho) hermitian semidefinite
W <= I
if nargin<3 || lower(op)=="ppt" || op==""
-I/m <= PartialTranspose(W,2,dim) <= I/m
elseif lower(op)=="pptp"
variable X(drho,drho) hermitian semidefinite
W <= X
PartialTranspose(X,2,dim) <= I/m
elseif lower(op)=="rains"
variable X(drho,drho) hermitian semidefinite
W <= X
-I/m <= PartialTranspose(X,2,dim) <= I/m
else
error('Unsupported type of operations. Allowed choices: ppt, pptp, rains.');
end
maximize trace(rho*W)
cvx_end
fid = cvx_optval;
end
[ent, ent0, W] = dist_ent (rho, tol, op, dim)
Computes the one-shot distillable entanglement under a given class of operations.
Required arguments:
rho: $n \times n$ density matrix or $n \times 1$ state vectortol: desired error tolerance in distillation fidelity, 0$\leq$tol$<$ 1Optional arguments:
op: class of operations (default: ppt)
Allowed choices:
ppt (PPT operations),
pptp ($PPT_+$-preserving operations),
rains (Rains-preserving operations)
Disregarded when input is pure (all classes including LOCC give the same result)
dim: a vector [dA, dB] of the dimensions of the subsystems (default: both systems of equal dimension)
Output:
ent: maximum size of maximally entangled state which can be obtained from rho under operations op up to tolerance tol, i.e. log2(ent) is the usual one-shot distillable entanglement.
The function makes an attempt to correct for CVX’s numerical accuracy (the default accuracy being within $\sqrt{\epsilon}$ of the optimal value, where $\epsilon$ is the machine precision); one should consult the exact value ent0 if needed.
ent0: the non-floored value of the optimisation problem such that log2(floor(1/ent0)) gives distillable entanglementW: optimal witness in the optimisation problem (see Prop. 4 of the manuscript)
If the input is provided as a vector, W is also returned as one.
function [ent, ent0, W] = dist_ent(rho, tol, op, dim)
dimrho = size(rho);
drho = max(dimrho);
if nargin < 4 % assume dimensions of subsystems are equal
dim = [sqrt(drho),sqrt(drho)];
d = sqrt(drho);
else
d = min(dim);
end
frho = full(rho);
if rank(frho) == 1 % pure state case
if min(dimrho) ~= 1
[~,~,v] = svd(frho);
rho = v(:,1);
end
[xi,U,V] = SchmidtDecomposition(rho, dim, -1);
cvx_begin quiet
variable x(d) nonnegative
norm(x) <= 1
x'*xi >= sqrt(1-tol);
minimize max(x)
cvx_end
if nargout > 2
% construct optimal witness
W = zeros(drho,1);
for j=1:d
W = W + x(j) * kron(U(:,j), V(:,j));
end
if min(dimrho) ~= 1
W = W*W';
end
end
ent0 = cvx_optval^2;
else % mixed state case
I = eye(drho);
rho = (rho+rho')/2; % avoids numerical problems
cvx_begin sdp quiet
variable W(drho,drho) hermitian semidefinite
W <= I
if nargin<3 || lower(op)=="ppt" || op==""
t = norm(PartialTranspose(W,2,dim));
elseif lower(op)=="pptp"
variable X(drho,drho) hermitian semidefinite
W <= X
t = lambda_max(PartialTranspose(X,2,dim));
elseif lower(op)=="rains"
variable X(drho,drho) hermitian semidefinite
W <= X
t = norm(PartialTranspose(X,2,dim));
else
error('Unsupported type of operations. Allowed choices: ppt, pptp, rains.');
end
trace(rho*W) >= 1-tol;
minimize t
cvx_end
ent0 = cvx_optval;
end
ent = floor(1/ent0);
prec = cvx_precision;
if abs(ent0 - 1/(ent+1)) <= prec(1) % to account for CVX precision
ent = ent+1;
end
end