Computational supplement

def distance(N,name): import cvxpy as cp J = herm(choi(N)-idchoi()) W = cp.Variable((4,4),hermitian=True) rho = cp.Variable((2,2),hermitian=True) primal = cp.Problem(cp.Maximize(cp.real(cp.trace(J@W))), [W>>0,rho>>0,cp.trace(rho)==1,cp.kron(I,rho)-W>>0]) lo = cp_solve(primal,name+"_diamond_primal") D = cp.Variable((4,4),hermitian=True) lam = cp.Variable() dual = cp.Problem(cp.Minimize(lam), [D>>0,D-J>>0,lam*I-(D[:2,:2]+D[2:,2:])>>0]) hi = cp_solve(dual,name+"_diamond_dual") rv = herm(D.value) shift = max(0.,-mineig(rv),-mineig(rv-J))+1e-9 rv += shift*np.eye(4) upper = float(np.linalg.eigvalsh( herm(partial_output(rv,2,2))).max()) resid = max(0.,-mineig(W.value),-mineig(rho.value), -mineig(np.kron(I,rho.value)-W.value), abs(np.trace(rho.value)-1)) need(resid<=STOL and abs(hi-lo)<=STOL,"diamond numerical gap "+name) save(name+"_distance",J=J,W=W.value,rho=rho.value, Z_raw=D.value,Z_repaired=rv) return "primal":lo,"dual":hi,"repaired_dual":upper, "primal_residual":float(resid),"gap":hi-lo, "certified":False

def recovery(N,name): import cvxpy as cp b = len(N[0,0]) R = cp.Variable((2*b,2*b),hermitian=True) mats = {} for i,j,E in units(): mats[i,j] = cp.bmat([[cp.reshape(cp.trace( R[a*b:(a+1)*b,z*b:(z+1)*b]@N[i,j].T),(1,1),order="C") for z in range(2)] for a in range(2)]) J = sum(cp.kron(mats[i,j],E) for i,j,E in units())-idchoi() D = cp.Variable((4,4),hermitian=True) lam = cp.Variable() prob = cp.Problem(cp.Minimize(lam), [R>>0,R[:b,:b]+R[b:,b:]==np.eye(b), D>>0,D-J>>0,lam*I-(D[:2,:2]+D[2:,2:])>>0]) estimate = cp_solve(prob,name+"_recovery") raw = R.value rr = max(0.,-mineig(raw), float(np.max(np.abs(partial_output(raw,2,b)-np.eye(b))))) need(rr<=STOL,"raw recovery residual "+name) fixed = repaired_decoder(raw) check(name+"_repaired_TP",partial_output(fixed,2,b),np.eye(b),STOL) need(mineig(fixed)>=-STOL,"repaired decoder positivity") achieved = apply_decoder(fixed,N) achieved_distance = distance(achieved,name+"_recovered")

C = sum((np.kron(E.T,N[i,j].T)/4 for i,j,E in units()), np.zeros((2*b,2*b),complex)) Yv = cp.Variable((b,b),hermitian=True) dual = cp.Problem(cp.Minimize(cp.real(cp.trace(Yv))), [cp.kron(I,Yv)-C>>0]) cp_solve(dual,name+"_fidelity_dual") Ynum = herm(Yv.value) shift = max(0.,-mineig(np.kron(I,Ynum)-C))+1e-9 Ynum += shift*np.eye(b) lower = max(0.,1-float(np.trace(Ynum).real)) need(lower<=achieved_distance["repaired_dual"]+STOL, "inconsistent recovery bounds") save(name+"_decoder",R_raw=raw,R_repaired=fixed,Y=Ynum, **{"unit_ return "estimate":estimate,"raw_residual":rr, "achieved_distance":achieved_distance, "fidelity_lower_candidate":lower,"certified":False