Phys. Rev. Lett 121, 010401 (2018)
Download Matlab files:
The computation of the quantities $T^{(m)}_{\mathcal{I}}$ for pure states reduces to a semi-analytical formula, as described in the manuscript. For general mixed states, each measure can be computed with a semidefinite program. The below Matlab implementation takes as input either pure state vectors or general density matrices, and for states of rank 2 or more also returns the optimal witness W. The semidefinite program is implemented using CVX.
function [val,W] = TI(rho,m)
% [val,W] TI(rho,m) computes the family of coherence measures as defined in PRL 121, 010401 (2018)
% rho: dxd density matrix or normalised state vector
% m: non-negative integer coefficient indicating the chosen measure
% the optimal witness W s.t. val=trace(rho*W) is returned only for non-pure rho
dim = size(rho);
d = max(dim);
if m>=d; m = d-1; end
if rank(rho) == 1 % pure state case
if min(dim) ~= 1
[~,~,v] = svd(rho);
rho = v(:,1);
end
m = m+1; % to compute the (m+1)-distillation norm
psi = sort(abs(rho),'descend');
k = 1;
minsum = sum(psi(m:d).^2);
for j=2:m
jsum = sum(psi(m-j+1:d).^2)/j;
if jsum < minsum
minsum = jsum;
k = j;
end
end
val = sum(psi(1:m-k))+k*sqrt(minsum);
val = val^2-1;
W = 0;
else % general case
rho = (rho+rho')/2; % to avoid numerical problems
cvx_begin quiet sdp
variable W(d,d) hermitian
-eye(d) <= W <= m*eye(d)
diag(W) <= 0
maximise trace(rho*W);
cvx_end
val = cvx_optval;
end
end
Notice that a simple change of the constraint diag(W) <= 0 into diag(W) == 0 gives instead the family of measures $T^{(m)}_{\mathcal{J}}$, which characterise the one-shot fidelity of coherence distillation.
Download the Matlab file:
The one-shot distillable coherence under MIO and DIO (and for pure states also SIO and IO) is computed through a semidefinite program. An example implementation with CVX is given below.
function dis = distillable_coherence(rho,tol)
% distillable_coherence(rho,tol) computes the one-shot distillable coherence under MIO/DIO
% rho: dxd density matrix or normalised state vector
% tol: error tolerance 0<=tol<=1
dim = size(rho);
d = max(dim);
if dim(1)==1; rho=rho'*rho;
elseif dim(2)==1; rho=rho*rho'; end
rho = (rho+rho')/2; % to avoid numerical problems
cvx_begin sdp quiet
variable G(d,d) hermitian
variable k nonnegative
0 <= G <= eye(d)
trace(rho*G) >= 1-tol;
diag(G) == k*ones(d,1)
minimise k
cvx_end
dis = log2(floor(1/cvx_optval));
end