initial commit

scripts/README.md

@@ -0,0 +1,81 @@
+# Reproducibility Scripts
+
+This directory contains small numerical scripts used to reproduce values quoted in the paper. Scripts should be executable from the repository root unless noted otherwise.
+
+## Requirements
+
+- Python 3
+- NumPy
+
+Install the only current dependency with:
+
+```bash
+python3 -m pip install numpy
+```
+
+## Current Scripts
+
+### `tiles_upb.py`
+
+Computes the two-qutrit Tiles unextendible-product-basis benchmark used in `paper/symmetric_shadow_maps_formal.tex`. The script uses Gell-Mann generators scaled so that `Tr(sigma_i sigma_j) = 3 delta_ij`, matching the paper's generator normalization.
+
+Run from the repository root:
+
+```bash
+python3 scripts/tiles_upb.py
+```
+
+The script constructs the five Tiles UPB product vectors, forms the normalized projector onto their four-dimensional orthogonal complement, and prints:
+
+- the UPB Gram matrix;
+- the density-matrix trace, Hermiticity check, and spectrum;
+- the partial-transpose spectrum and minimum eigenvalue;
+- the shadow-map value `||M_A(rho)||_*`, normalized by `sqrt((3-1)(3-1)) = 2`;
+- the CCNR/realignment trace norm for comparison.
+
+The values used in `paper/symmetric_shadow_maps_formal.tex` are:
+
+| quantity | value |
+| --- | ---: |
+| minimum eigenvalue of partial transpose | `-1.5922869149236308e-16` |
+| unnormalized correlation nuclear norm | `2.1068432645403345` |
+| normalized shadow-map value | `1.0534216322701673` |
+| CCNR/realignment trace norm | `1.087412464837521` |
+
+The tiny negative partial-transpose eigenvalue is numerical roundoff; the state is the standard PPT-entangled Tiles UPB state. The paper compares the normalized shadow-map value to the separability bound `<= 1`.
+
+### `grraph_state_cuts.py`
+
+Computes the cut-resolved bigraduated shadow-map values for the four-qubit ring graph state with edges `(1,2)`, `(2,3)`, `(3,4)`, `(4,1)`. Internally the script uses zero-based qubit labels.
+
+Run from the repository root:
+
+```bash
+python3 scripts/grraph_state_cuts.py
+```
+
+The script constructs the graph state, computes its Pauli correlation tensor, forms the normalized `2|2` source-target unfoldings, and prints:
+
+- the full normalized nuclear norm `||M_S||_*` for each cut;
+- the normalized full-sector block, obtained by keeping only nonidentity Pauli labels on both source qubits and both target qubits;
+- the eigenvalues of the two-qubit source marginal and its distance from the maximally mixed state.
+
+The values used in `paper/symmetric_shadow_maps_formal.tex` are:
+
+| cut | full normalized norm | normalized full-sector block | source marginal |
+| --- | ---: | ---: | --- |
+| adjacent `{1,2}|{3,4}` | `5` | `5/3` | maximally mixed |
+| diagonal `{1,3}|{2,4}` | `7/3` | `1` | eigenvalues `1/2, 1/2, 0, 0` |
+| adjacent `{1,4}|{2,3}` | `5` | `5/3` | maximally mixed |
+
+The script also prints raw, unnormalized nuclear norms. The paper compares only the normalized values to the separability bound `<= 1`.
+
+## Adding New Scripts
+
+When adding another script, add a short entry above with:
+
+- the purpose of the script;
+- the command needed to run it from the repository root;
+- required dependencies beyond NumPy, if any;
+- the paper values, table, or figure it reproduces;
+- notes about normalization conventions if the output includes both raw and normalized quantities.

scripts/grraph_state_cuts.py

@@ -0,0 +1,83 @@
+import numpy as np
+from itertools import product, combinations
+
+I=np.array([[1,0],[0,1]],complex)
+X=np.array([[0,1],[1,0]],complex)
+Y=np.array([[0,-1j],[1j,0]],complex)
+Z=np.array([[1,0],[0,-1]],complex)
+paulis=[I,X,Y,Z]
+labels=['I','X','Y','Z']
+
+def kron_all(ops):
+    out=ops[0]
+    for op in ops[1:]: out=np.kron(out,op)
+    return out
+
+def graph_state(n, edges):
+    psi=np.ones(2**n,complex)/np.sqrt(2**n)
+    # basis index bits qubit 0..n-1 as MSB? phase invariant consistent
+    for idx in range(2**n):
+        bits=[(idx>>(n-1-q))&1 for q in range(n)]
+        phase=1
+        for a,b in edges:
+            if bits[a]*bits[b]: phase*=-1
+        psi[idx]*=phase
+    return psi
+
+def corr_tensor(psi,n):
+    coeff={}
+    for inds in product(range(4), repeat=n):
+        op=kron_all([paulis[i] for i in inds])
+        val=np.vdot(psi, op@psi)
+        if abs(val)>1e-9:
+            coeff[inds]=float(np.real_if_close(val))
+    return coeff
+
+def matrix_for_cut(coeff,S,n, source_full=False, target_full=False):
+    Sc=[i for i in range(n) if i not in S]
+    row=[]; col=[]
+    for inds in product(range(4), repeat=len(S)):
+        if all(i==0 for i in inds): continue
+        if source_full and any(i==0 for i in inds): continue
+        row.append(inds)
+    for inds in product(range(4), repeat=len(Sc)):
+        if all(i==0 for i in inds): continue
+        if target_full and any(i==0 for i in inds): continue
+        col.append(inds)
+    M=np.zeros((len(col),len(row))) # target rows, source cols
+    for r,tinds in enumerate(col):
+        for c,sinds in enumerate(row):
+            full=[0]*n
+            for q,ind in zip(S,sinds): full[q]=ind
+            for q,ind in zip(Sc,tinds): full[q]=ind
+            M[r,c]=coeff.get(tuple(full),0.0)
+    return M,row,col
+
+def partial_rho(psi, keep, n):
+    rho=np.outer(psi, psi.conj()).reshape([2]*n*2)
+    trace=[i for i in range(n) if i not in keep]
+    # trace out from high to low axes
+    for q in sorted(trace, reverse=True):
+        rho=np.trace(rho, axis1=q, axis2=q+rho.ndim//2)
+    d=2**len(keep)
+    return rho.reshape(d,d)
+
+n=4
+edges=[(0,1),(1,2),(2,3),(3,0)]
+psi=graph_state(n,edges)
+coeff=corr_tensor(psi,n)
+print('nonzero coeffs', len(coeff))
+for cutname,S in [('adjacent12|34',[0,1]),('diagonal13|24',[0,2]),('adjacent14|23',[0,3])]:
+    M,_,_=matrix_for_cut(coeff,S,n)
+    G,_,_=matrix_for_cut(coeff,S,n,source_full=True,target_full=True)
+    norm=np.linalg.svd(M/3, compute_uv=False).sum()
+    rawn=np.linalg.svd(M, compute_uv=False).sum()
+    gnorm=np.linalg.svd(G/3, compute_uv=False).sum()
+    graw=np.linalg.svd(G, compute_uv=False).sum()
+    print('\n',cutname)
+    print('full normalized',norm,'raw',rawn,'sing',np.linalg.svd(M/3,compute_uv=False))
+    print('genuine normalized',gnorm,'raw',graw,'sing',np.linalg.svd(G/3,compute_uv=False))
+    rho=partial_rho(psi,S,n)
+    print('rho eigen',np.linalg.eigvalsh(rho),'max mixed dist',np.linalg.norm(rho-np.eye(4)/4))
+    # list nonzero stabilizers labels
+    # print rows/cols maybe

scripts/tiles_upb.py

@@ -0,0 +1,105 @@
+import numpy as np
+
+np.set_printoptions(precision=5, suppress=True)
+
+# ---------- Gell-Mann generators for d=3, normalized so Tr(sigma_i sigma_j) = 3*delta_ij ----------
+i_ = 1j
+lam = [None]*9
+lam[1] = np.array([[0,1,0],[1,0,0],[0,0,0]], dtype=complex)
+lam[2] = np.array([[0,-i_,0],[i_,0,0],[0,0,0]], dtype=complex)
+lam[3] = np.array([[1,0,0],[0,-1,0],[0,0,0]], dtype=complex)
+lam[4] = np.array([[0,0,1],[0,0,0],[1,0,0]], dtype=complex)
+lam[5] = np.array([[0,0,-i_],[0,0,0],[i_,0,0]], dtype=complex)
+lam[6] = np.array([[0,0,0],[0,0,1],[0,1,0]], dtype=complex)
+lam[7] = np.array([[0,0,0],[0,0,-i_],[0,i_,0]], dtype=complex)
+lam[8] = (1/np.sqrt(3))*np.array([[1,0,0],[0,1,0],[0,0,-2]], dtype=complex)
+
+# check standard normalization Tr(lam_a lam_b) = 2 delta_ab
+for a in range(1,9):
+    for b in range(1,9):
+        val = np.trace(lam[a]@lam[b])
+        if a==b and not np.isclose(val,2):
+            print("WARN std norm", a,b,val)
+        if a!=b and not np.isclose(val,0):
+            print("WARN std orth", a,b,val)
+
+sigma = [None] + [np.sqrt(3/2)*lam[k] for k in range(1,9)]  # d=3 -> Tr(sigma_i sigma_j)=3 delta_ij
+
+# sanity check
+for a in range(1,9):
+    for b in range(1,9):
+        val = np.trace(sigma[a]@sigma[b]).real
+        expected = 3.0 if a==b else 0.0
+        assert abs(val-expected) < 1e-9, (a,b,val)
+print("Generator normalization OK: Tr(sigma_i sigma_j) = 3 delta_ij")
+
+# ---------- Tiles UPB (Bennett, DiVincenzo, Mor, Shor, Smolin, Terhal 1999) ----------
+e0 = np.array([1,0,0], dtype=complex)
+e1 = np.array([0,1,0], dtype=complex)
+e2 = np.array([0,0,1], dtype=complex)
+
+def nrm(v):
+    return v/np.linalg.norm(v)
+
+psi = []
+psi.append(np.kron(e0, nrm(e0-e1)))
+psi.append(np.kron(e2, nrm(e1-e2)))
+psi.append(np.kron(nrm(e0-e1), e2))
+psi.append(np.kron(nrm(e1-e2), e0))
+psi.append(np.kron(nrm(e0+e1+e2), nrm(e0+e1+e2)))
+
+# check orthonormality
+G = np.array([[np.vdot(p,q) for q in psi] for p in psi])
+print("\nGram matrix of the 5 UPB vectors (should be identity):")
+print(np.round(G,6))
+
+P = sum(np.outer(p, p.conj()) for p in psi)
+I9 = np.eye(9, dtype=complex)
+rho = (I9 - P)/4.0
+
+print("\nTr(rho) =", np.trace(rho).real, " (should be 1)")
+print("rho is Hermitian:", np.allclose(rho, rho.conj().T))
+eigvals_rho = np.linalg.eigvalsh(rho)
+print("eigenvalues of rho (should be >=0, rank 4 nonzero):", np.round(eigvals_rho,5))
+
+# ---------- PPT check ----------
+def partial_transpose_B(rho, dA=3, dB=3):
+    r = rho.reshape(dA,dB,dA,dB)
+    rpt = r.transpose(0,3,2,1)
+    return rpt.reshape(dA*dB, dA*dB)
+
+rho_pt = partial_transpose_B(rho)
+eig_pt = np.linalg.eigvalsh(rho_pt)
+print("\nEigenvalues of partial transpose (PPT check):")
+print(np.round(eig_pt,6))
+print("min eigenvalue of PT:", eig_pt.min(), " -> PPT" if eig_pt.min() > -1e-9 else " -> NPT (entangled via ordinary PPT already)")
+
+# ---------- correlation tensor / shadow map ----------
+T = np.zeros((8,8))
+for a in range(1,9):
+    for b in range(1,9):
+        op = np.kron(sigma[a], sigma[b])
+        T[a-1,b-1] = np.trace(rho @ op).real
+
+s = np.linalg.svd(T, compute_uv=False)
+nuclear_T = s.sum()
+dA, dB = 3, 3
+norm_const = np.sqrt((dA-1)*(dB-1))
+M_norm = nuclear_T / norm_const
+
+print("\nSingular values of correlation tensor T:", np.round(s,5))
+print("Nuclear norm ||T||_* =", nuclear_T)
+print("Normalization constant sqrt((dA-1)(dB-1)) =", norm_const)
+print("Shadow-map value ||M_A(rho)||_* =", M_norm, " (separable bound: <= 1)")
+
+# ---------- CCNR / realignment criterion for comparison ----------
+def realign(rho, dA=3, dB=3):
+    r = rho.reshape(dA,dB,dA,dB)
+    # standard realignment: R_{(i mu),(j nu)} = rho_{ij,mu nu}
+    R = r.transpose(0,2,1,3).reshape(dA*dA, dB*dB)
+    return R
+
+R = realign(rho)
+s_R = np.linalg.svd(R, compute_uv=False)
+ccnr = s_R.sum()
+print("\nCCNR (realignment) trace norm:", ccnr, " (separable bound: <= 1)")