quantum-shadow-maps_v2/scripts/dps_hierarchy/03_dps_level2_demo.py

55 lines
2.1 KiB
Python
Raw Permalink Normal View History

2026-07-26 14:09:49 +02:00
"""
03_dps_level2_demo.py
Demonstrates the level-2 DPS SDP (via dps_hierarchy.build_dps_problem) on
two test states:
A) the qutrit Werner state -- sanity check. PPT is already exactly
tight for this family (p_c=1/4, see 02_werner_qutrit_symbolic.py), so
DPS-2 cannot improve on it; well away from the boundary both should
agree.
B) the Tiles UPB bound-entangled state -- the interesting case. Plain
PPT is blind (min eigenvalue ~0, "PPT to machine precision" as noted
in the paper's own Tiles example), but DPS level 2 correctly detects
the entanglement.
Expected runtime: well under a minute with SCS.
"""
import cvxpy as cp
from common import RHO_TILES, werner_qutrit, plain_ppt_feasible, plain_ppt_min_eig
from dps_hierarchy import build_dps_problem, dps_feasible
SOLVER = cp.SCS # swap to cp.MOSEK if you have a license -- likely much faster
prob, rho_param, sigma = build_dps_problem(d=3, k=2)
print(f"DPS level-2 SDP built: sigma shape {sigma.shape}\n")
print("=" * 70)
print("Sanity check: qutrit Werner state (known exact threshold p=1/4)")
print("=" * 70)
for p in [0.20, 0.40]:
rho = werner_qutrit(p)
ppt_ok = plain_ppt_feasible(rho)
eig = plain_ppt_min_eig(rho)
ok = dps_feasible(prob, rho_param, rho, solver=SOLVER, eps=1e-6)
print(f"p={p:.2f}: PPT feasible={ppt_ok} (min eig {eig:+.5f}) "
f"DPS-2 feasible={ok}")
print("\n" + "=" * 70)
print("Tiles UPB bound-entangled state (the interesting case)")
print("=" * 70)
ppt_ok = plain_ppt_feasible(RHO_TILES)
eig = plain_ppt_min_eig(RHO_TILES)
print(f"Plain PPT feasible: {ppt_ok} (min eigenvalue: {eig:.10f})")
ok = dps_feasible(prob, rho_param, RHO_TILES, solver=SOLVER, eps=1e-6)
print(f"DPS level-2 feasible: {ok} "
f"({'NOT detected' if ok else 'ENTANGLEMENT DETECTED'})")
print("\nRobustness across solver tolerances (guards against SDP numerical "
"artifacts near a threshold -- see the chat for a case where this "
"mattered):")
for eps in [1e-5, 1e-6, 1e-7, 1e-8]:
dps_feasible(prob, rho_param, RHO_TILES, solver=SOLVER, eps=eps, max_iters=50000)
print(f" eps={eps}: status={prob.status}")