Designing Quantum Error-Correcting Codes with Linear Programming

October 16, 2026

The Quantum Delsarte Bound and a Stabilizer Code Search

How many logical qubits can we protect, and how well, using a given number of physical qubits? An $[[n,k,d]]$ quantum code encodes $k$ logical qubits into $n$ physical qubits and corrects any error acting on up to $\lfloor (d-1)/2 \rfloor$ qubits. Two questions drive code design:

Upper bounds. Which triples $(n,k,d)$ are impossible?
Construction. For the triples that are not ruled out, can we actually find a code?

In this article we answer both with a single concrete example. First, we use the quantum Delsarte (quantum linear programming) bound to compute, for every $n \le 16$ and $k < n$, the largest distance $d$ that is not excluded. Second, we run a random stabilizer-code search for $[[5,1,d]]$ and check that the code it finds is consistent with the LP constraints.

1. Mathematical Background

1.1 Weight enumerators of a stabilizer code

Let $\mathcal{S}$ be the stabilizer group of an $[[n,k]]$ stabilizer code, so $|\mathcal{S}| = 2^{n-k}$. Let $\mathcal{N}$ be its normalizer, so $|\mathcal{N}| = 2^{n+k}$. Define

where $\mathrm{wt}(E)$ is the number of qubits on which the Pauli operator $E$ acts non-trivially. Then $A_0 = B_0 = 1$ and

$$
\sum_{j=0}^{n} A_j = 2^{n-k}, \qquad \sum_{j=0}^{n} B_j = 2^{n+k}.
$$

1.2 Quantum MacWilliams identity

With $K = 2^k$, the two enumerators are linked by

$$
B_j = \frac{K}{2^n} \sum_{i=0}^{n} Q_j(i), A_i, \qquad
Q_j(i) = \sum_{s=0}^{j} (-1)^s, 3^{,j-s} \binom{i}{s}\binom{n-i}{j-s},
$$

where $Q_j(i)$ is the quaternary Krawtchouk polynomial.

1.3 The shadow enumerator

Rains showed that the following quantity must also be non-negative:

$$
S_j = \frac{K}{2^n} \sum_{i=0}^{n} (-1)^i, Q_j(i), A_i ;\ge; 0 .
$$

1.4 The linear program

If an $[[n,k,d]]$ stabilizer code exists, the vector $A = (A_0,\dots,A_n)$ must satisfy:

$$
\begin{aligned}
& A_0 = 1, \qquad \sum_{j} A_j = 2^{n-k}, \qquad A_j \ge 0,\
& B_j = A_j \quad (1 \le j \le d-1) \qquad \text{(no undetectable error of weight} < d \text{ outside } \mathcal{S}),\
& B_j \ge A_j \quad (j \ge d),\
& S_j \ge 0 \quad (0 \le j \le n).
\end{aligned}
$$

All of these constraints are linear in $A$. If the system is infeasible, no $[[n,k,d]]$ code exists. This is the quantum Delsarte (quantum LP) bound. The condition is necessary but not sufficient, so feasibility does not guarantee existence.

Feasibility is monotone in $d$: raising $d$ only adds constraints. For each $(n,k)$ we can therefore increase $d$ until the LP first becomes infeasible. The last feasible value is $d_{\mathrm{LP}}(n,k)$.

We compare it with the quantum Singleton bound:

$$
n - k ;\ge; 2(d-1) \quad\Longleftrightarrow\quad d \le \left\lfloor \frac{n-k}{2} \right\rfloor + 1 .
$$

1.5 The code search

Every $n$-qubit Pauli operator (up to phase) is a vector $(x,|,z) \in \mathbb{F}_2^{2n}$. Two Paulis commute exactly when their symplectic product vanishes:

$$
\langle (x|z),(x’|z’) \rangle = x\cdot z’ + z\cdot x’ \equiv 0 \pmod 2 .
$$

A stabilizer code is an $(n-k)$-dimensional isotropic subspace of $\mathbb{F}_2^{2n}$. We build one at random as follows:

Pick a Pauli that commutes with all earlier generators and is not already in their span.
Repeat until we have $n-k$ generators.
Compute $A_j$ and $B_j$ by direct enumeration.
Read off the distance as the minimum weight in $\mathcal{N}\setminus\mathcal{S}$.

For $[[5,1]]$ the famous perfect code $[[5,1,3]]$ should appear. Its enumerators are

$$
A(z) = 1 + 15 z^4, \qquad B(z) = 1 + 30 z^3 + 15 z^4 + 18 z^5 .
$$

2. The Concrete Problem

Task A (bound). For every $2 \le n \le 16$ and $0 \le k \le n-1$, determine

$$
d_{\mathrm{LP}}(n,k) = \max{, d : \text{the LP above is feasible} ,}
$$

and compare it with the Singleton bound.

Task B (search). Sample 4000 random $[[5,1]]$ stabilizer codes, record their distances, and take the best one. Then verify that its enumerators satisfy the quantum MacWilliams identity and the shadow inequalities, and that the LP is feasible for $d=3$ and infeasible for $d=4$.

3. Full Source Code

import numpy as np
import matplotlib.pyplot as plt
from math import comb
from scipy.optimize import linprog
from mpl_toolkits.mplot3d import Axes3D  # noqa: F401

plt.style.use('dark_background')

# ----------------------------------------------------------------------
# 1. Krawtchouk matrices (computed once per n with exact integer arithmetic)
# ----------------------------------------------------------------------
NMAX = 16

def krawtchouk_matrix(n):
    Q = np.zeros((n + 1, n + 1), dtype=float)
    for j in range(n + 1):
        for i in range(n + 1):
            total = 0
            for s in range(j + 1):
                total += (-1) ** s * 3 ** (j - s) * comb(i, s) * comb(n - i, j - s)
            Q[j, i] = float(total)
    return Q

KRAW = {n: krawtchouk_matrix(n) for n in range(2, NMAX + 1)}

# ----------------------------------------------------------------------
# 2. Quantum LP feasibility test
# ----------------------------------------------------------------------
def lp_feasible(n, k, d):
    Q = KRAW[n]
    scale = 2.0 ** (k - n)
    sign = (-1.0) ** np.arange(n + 1)
    Bmat = scale * Q
    Smat = scale * Q * sign[None, :]
    I = np.eye(n + 1)

    eq_rows = [I[0], np.ones(n + 1)]
    eq_rhs = [1.0, 2.0 ** (n - k)]
    for j in range(1, d):
        eq_rows.append(Bmat[j] - I[j])
        eq_rhs.append(0.0)

    ub_rows = []
    for j in range(d, n + 1):
        ub_rows.append(I[j] - Bmat[j])
    for j in range(n + 1):
        ub_rows.append(-Smat[j])

    A_eq = np.array(eq_rows)
    b_eq = np.array(eq_rhs)
    A_ub = np.array(ub_rows)
    b_ub = np.zeros(len(ub_rows))

    m = np.abs(A_eq).max(axis=1)
    m[m == 0] = 1.0
    A_eq = A_eq / m[:, None]
    b_eq = b_eq / m

    m = np.abs(A_ub).max(axis=1)
    m[m == 0] = 1.0
    A_ub = A_ub / m[:, None]

    res = linprog(np.zeros(n + 1), A_ub=A_ub, b_ub=b_ub,
                  A_eq=A_eq, b_eq=b_eq,
                  bounds=[(0, None)] * (n + 1), method='highs')
    ok = (res.status == 0)
    return ok, (res.x if ok else None)

def singleton(n, k):
    return min(n, (n - k) // 2 + 1)

# ----------------------------------------------------------------------
# 3. Task A: LP bound for every (n, k)
# ----------------------------------------------------------------------
dLP = {}
for n in range(2, NMAX + 1):
    for k in range(0, n):
        best = 0
        for d in range(1, n + 1):
            ok, _ = lp_feasible(n, k, d)
            if ok:
                best = d
            else:
                break
        dLP[(n, k)] = best

print("Quantum LP upper bound d_LP(n, k)")
header = "  n \\ k | " + " ".join(f"{k:2d}" for k in range(0, NMAX))
print(header)
print("-" * len(header))
for n in range(2, NMAX + 1):
    row = []
    for k in range(0, NMAX):
        row.append(f"{dLP[(n, k)]:2d}" if k < n else " .")
    print(f"  {n:4d}  | " + " ".join(row))

print()
ok3, x3 = lp_feasible(5, 1, 3)
ok4, _ = lp_feasible(5, 1, 4)
print("LP feasible for [[5,1,3]] :", ok3)
print("LP feasible for [[5,1,4]] :", ok4)
if ok3:
    print("LP solution A_j for [[5,1,3]] :", np.round(x3, 6))

# ----------------------------------------------------------------------
# 4. Task B: random stabilizer-code search for [[5,1,d]]
# ----------------------------------------------------------------------
def build_pauli_tables(n):
    idx = np.arange(4 ** n)
    digits = (idx[:, None] // (4 ** np.arange(n))[None, :]) % 4   # 0:I 1:X 2:Z 3:Y
    X = (digits & 1).astype(np.int64)
    Z = (digits >> 1).astype(np.int64)
    pw = 2 ** np.arange(n)
    codes = X @ pw + ((Z @ pw) << n)
    wt = (digits > 0).sum(axis=1)
    return X, Z, codes, wt

def random_stabilizer_code(n, k, rng, tables):
    X, Z, codes, wt = tables
    elems = np.array([0], dtype=np.int64)
    allowed = np.ones(len(codes), dtype=bool)
    for _ in range(n - k):
        cand = np.flatnonzero(allowed & ~np.isin(codes, elems))
        g = rng.choice(cand)
        comm = (X @ Z[g] + Z @ X[g]) % 2
        allowed &= (comm == 0)
        elems = np.concatenate([elems, elems ^ codes[g]])
    in_S = np.isin(codes, elems)
    in_N = allowed
    A = np.bincount(wt[in_S], minlength=n + 1)
    B = np.bincount(wt[in_N], minlength=n + 1)
    diff = in_N & ~in_S
    d = int(wt[diff].min()) if diff.any() else n
    return A, B, d

n_s, k_s, trials = 5, 1, 4000
rng = np.random.default_rng(2024)
tables = build_pauli_tables(n_s)

ds = []
bestA, bestB, bestd = None, None, -1
for _ in range(trials):
    A, B, d = random_stabilizer_code(n_s, k_s, rng, tables)
    ds.append(d)
    if d > bestd:
        bestA, bestB, bestd = A, B, d

vals, counts = np.unique(ds, return_counts=True)
print()
print(f"Random search over {trials} stabilizer codes with n={n_s}, k={k_s}")
for v, c in zip(vals, counts):
    print(f"  distance d = {v} : {c} codes ({100.0 * c / trials:.2f} %)")
print("Best code distance      :", bestd)
print("A_j (stabilizer)        :", bestA)
print("B_j (normalizer)        :", bestB)

Q5 = KRAW[n_s]
B_pred = 2.0 ** (k_s - n_s) * (Q5 @ bestA)
S_found = 2.0 ** (k_s - n_s) * (Q5 @ (((-1.0) ** np.arange(n_s + 1)) * bestA))
print("MacWilliams max error   :", float(np.abs(B_pred - bestB).max()))
print("min (B_j - A_j)         :", int((bestB - bestA).min()))
print("min shadow S_j          :", float(S_found.min()))
okb, _ = lp_feasible(n_s, k_s, bestd)
okn, _ = lp_feasible(n_s, k_s, bestd + 1)
print(f"LP feasible at d={bestd}   :", okb)
print(f"LP feasible at d={bestd + 1}   :", okn)

# ----------------------------------------------------------------------
# 5. One combined figure (six panels, two of them 3D)
# ----------------------------------------------------------------------
keys = sorted(dLP.keys())
xs = np.array([n for (n, k) in keys], dtype=float)
ys = np.array([k for (n, k) in keys], dtype=float)
zs = np.array([dLP[key] for key in keys], dtype=float)
gaps = np.array([max(singleton(n, k) - dLP[(n, k)], 0) for (n, k) in keys], dtype=float)

fig = plt.figure(figsize=(21, 12))
fig.suptitle("Quantum Delsarte (LP) Bound and Stabilizer Code Search", fontsize=18)

# (1) 3D bars: LP bound
ax1 = fig.add_subplot(2, 3, 1, projection='3d')
cols = plt.cm.plasma(zs / max(zs.max(), 1.0))
ax1.bar3d(xs - 0.4, ys - 0.4, np.zeros_like(zs), 0.8, 0.8, zs, color=cols, shade=True)
ax1.set_xlabel("n (physical qubits)")
ax1.set_ylabel("k (logical qubits)")
ax1.set_zlabel(r"$d_{\mathrm{LP}}$")
ax1.set_title(r"LP upper bound $d_{\mathrm{LP}}(n,k)$")
ax1.view_init(elev=28, azim=-60)

# (2) 3D bars: gap to Singleton
ax2 = fig.add_subplot(2, 3, 2, projection='3d')
cols2 = plt.cm.cool(gaps / max(gaps.max(), 1.0))
ax2.bar3d(xs - 0.4, ys - 0.4, np.zeros_like(gaps), 0.8, 0.8, gaps, color=cols2, shade=True)
ax2.set_xlabel("n")
ax2.set_ylabel("k")
ax2.set_zlabel("Singleton - LP")
ax2.set_title("Improvement of the LP bound over Singleton")
ax2.view_init(elev=28, azim=-60)

# (3) Heatmap
ax3 = fig.add_subplot(2, 3, 3)
M = np.full((NMAX, NMAX - 1), np.nan)
for (n, k), v in dLP.items():
    M[k, n - 2] = v
im = ax3.imshow(np.ma.masked_invalid(M), origin='lower', aspect='auto',
                cmap='viridis', extent=(1.5, NMAX + 0.5, -0.5, NMAX - 0.5))
for (n, k), v in dLP.items():
    ax3.text(n, k, str(v), ha='center', va='center', fontsize=7, color='white')
fig.colorbar(im, ax=ax3, fraction=0.046)
ax3.set_xlabel("n")
ax3.set_ylabel("k")
ax3.set_title(r"Heatmap of $d_{\mathrm{LP}}(n,k)$")

# (4) Lines
ax4 = fig.add_subplot(2, 3, 4)
palette = ['#ff6b6b', '#4ecdc4', '#ffe66d', '#a29bfe']
for color, k in zip(palette, [0, 1, 2, 3]):
    ns = [n for n in range(2, NMAX + 1) if k < n]
    ax4.plot(ns, [dLP[(n, k)] for n in ns], 'o-', color=color, label=f"LP, k={k}")
    ax4.plot(ns, [singleton(n, k) for n in ns], '--', color=color, alpha=0.5)
ax4.set_xlabel("n")
ax4.set_ylabel("distance d")
ax4.set_title("LP bound (solid) vs Singleton bound (dashed)")
ax4.grid(alpha=0.3)
ax4.legend()

# (5) Enumerators of the best code found
ax5 = fig.add_subplot(2, 3, 5)
js = np.arange(n_s + 1)
ax5.bar(js - 0.2, bestA, 0.4, color='#ff9f43', label=r"$A_j$ (stabilizer)")
ax5.bar(js + 0.2, bestB, 0.4, color='#54a0ff', label=r"$B_j$ (normalizer)")
ax5.set_xlabel("weight j")
ax5.set_ylabel("count")
ax5.set_title(r"Enumerators of the best code found: $[[5,1,%d]]$" % bestd)
ax5.legend()
ax5.grid(alpha=0.3)

# (6) Distance histogram
ax6 = fig.add_subplot(2, 3, 6)
ax6.bar(vals, counts, color='#1dd1a1')
for v, c in zip(vals, counts):
    ax6.text(v, c, str(c), ha='center', va='bottom')
ax6.set_xticks(vals)
ax6.set_xlabel("distance d")
ax6.set_ylabel("number of codes")
ax6.set_title("Distances of 4000 random [[5,1]] stabilizer codes")
ax6.grid(alpha=0.3)

plt.tight_layout(rect=[0, 0, 1, 0.95])
plt.show()

4. Code Walkthrough

4.1 Krawtchouk matrices

krawtchouk_matrix(n) evaluates $Q_j(i)$ directly from its definition, using Python’s exact integers and math.comb. comb returns $0$ when the lower index exceeds the upper, so the boundary cases of the double binomial sum need no special handling. The matrices are computed once for each $n$ and cached in KRAW, so the thousands of LP calls never recompute them.

4.2 The LP feasibility test

lp_feasible(n, k, d) builds exactly the constraints of Section 1.4 in matrix form.

Bmat is $\frac{K}{2^n} Q$, so $B = \texttt{Bmat},A$. Smat is the same matrix with the alternating sign $(-1)^i$ applied column-wise, which gives the shadow enumerator.
The equality block contains $A_0 = 1$, $\sum_j A_j = 2^{n-k}$, and $B_j - A_j = 0$ for $1 \le j \le d-1$.
The inequality block contains $A_j - B_j \le 0$ for $j \ge d$, and $-S_j \le 0$ for all $j$.
Each row is divided by its largest absolute coefficient. Without this, entries range from $1$ to about $10^{9}$ and the solver’s tolerances become unreliable.
The objective is identically zero, because we only ask whether a feasible point exists. HiGHS reports status $0$ for a feasible problem, and any other status (in particular status $2$, infeasible) is treated as “no code”.

4.3 Speed

The expensive-looking part is the double loop over $(n,k)$ together with an inner loop over $d$. Three choices keep it fast:

The Krawtchouk matrices are cached, as above.
The inner loop stops at the first infeasible $d$, which is valid because feasibility is monotone in $d$.
Each LP has at most $17$ variables, and HiGHS solves it in milliseconds.

The code search is vectorized with NumPy. All $4^n$ Pauli operators are tabulated once as bit vectors. A Pauli is encoded as the integer X | (Z << n), so multiplication of Paulis (ignoring phase) is a single XOR, and commutation with a generator is one matrix–vector product modulo $2$.

4.4 Random stabilizer code

random_stabilizer_code adds generators one at a time:

allowed marks all Paulis that commute with every generator chosen so far. np.isin(codes, elems) removes those already in the generated group.
A random candidate $g$ is drawn from what remains. A candidate always exists while fewer than $n$ generators have been chosen, because the normalizer is strictly larger than the stabilizer.
After choosing $g$, allowed is intersected with the commutant of $g$, and the group elems is doubled by XOR with the code of $g$.
At the end, in_S marks the stabilizer and in_N (the final allowed mask) marks the normalizer. np.bincount over the weights gives $A_j$ and $B_j$.
The distance is the minimum weight over $\mathcal{N}\setminus\mathcal{S}$.

No rejection sampling is needed, so every trial succeeds.

4.5 Consistency checks

After the search, the best code is tested against the theory:

the quantum MacWilliams identity, $B = \frac{K}{2^n} Q A$ (the maximum error should be at rounding level),
the inequality $B_j \ge A_j$,
non-negativity of the shadow enumerator,
LP feasibility at the found distance $d$ and, as a sanity check, infeasibility at $d+1$.

4.6 The figure

All plots are drawn in one figure with a $2 \times 3$ grid, so the output is a single image.

5. Execution Results

Quantum LP upper bound d_LP(n, k)
  n \ k |  0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15
---------------------------------------------------------
     2  |  2  1  .  .  .  .  .  .  .  .  .  .  .  .  .  .
     3  |  3  1  1  .  .  .  .  .  .  .  .  .  .  .  .  .
     4  |  4  2  2  1  .  .  .  .  .  .  .  .  .  .  .  .
     5  |  5  3  2  1  1  .  .  .  .  .  .  .  .  .  .  .
     6  |  6  3  2  2  2  1  .  .  .  .  .  .  .  .  .  .
     7  |  7  3  2  2  2  1  1  .  .  .  .  .  .  .  .  .
     8  |  8  3  3  3  2  2  2  1  .  .  .  .  .  .  .  .
     9  |  9  3  3  3  2  2  2  1  1  .  .  .  .  .  .  .
    10  | 10  4  4  3  3  2  2  2  2  1  .  .  .  .  .  .
    11  | 11  5  4  3  3  3  2  2  2  1  1  .  .  .  .  .
    12  | 12  5  4  4  4  3  3  2  2  2  2  1  .  .  .  .
    13  | 13  5  4  4  4  4  3  3  2  2  2  1  1  .  .  .
    14  | 14  5  5  5  4  4  4  3  3  2  2  2  2  1  .  .
    15  | 15  5  5  5  5  4  4  4  3  3  2  2  2  1  1  .
    16  | 16  6  6  5  5  5  4  4  4  3  3  2  2  2  2  1

LP feasible for [[5,1,3]] : True
LP feasible for [[5,1,4]] : False
LP solution A_j for [[5,1,3]] : [ 1.  0.  0.  0. 15. -0.]

Random search over 4000 stabilizer codes with n=5, k=1
  distance d = 1 : 2321 codes (58.02 %)
  distance d = 2 : 1663 codes (41.58 %)
  distance d = 3 : 16 codes (0.40 %)
Best code distance      : 3
A_j (stabilizer)        : [ 1  0  0  0 15  0]
B_j (normalizer)        : [ 1  0  0 30 15 18]
MacWilliams max error   : 0.0
min (B_j - A_j)         : 0
min shadow S_j          : 0.0
LP feasible at d=3   : True
LP feasible at d=4   : False

6. Interpreting the Results

6.1 The 3D bar chart of $d_{\mathrm{LP}}(n,k)$ (top left)

The horizontal axes are the number of physical qubits $n$ and logical qubits $k$, and the height of each bar is $d_{\mathrm{LP}}(n,k)$. Two trends are visible:

Along the $n$ axis the bars rise: more physical qubits allow a larger distance.
Along the $k$ axis the bars fall: encoding more logical qubits into the same $n$ qubits costs protection.

This is the basic rate–distance trade-off of quantum coding. The LP bound shows exactly where the wall of impossibility lies.

6.2 The gap to the Singleton bound (top middle)

Here each bar is the amount by which the LP bound beats the quantum Singleton bound, $\lfloor (n-k)/2 \rfloor + 1 - d_{\mathrm{LP}}$. Where a bar is zero, the LP adds nothing beyond Singleton. Where a bar is positive, the LP has excluded distances that Singleton still allowed. The Singleton bound is a counting argument that ignores the finer structure of the weight distribution. The LP uses the full MacWilliams and shadow structure, so it is never weaker, and it becomes noticeably stronger as $n$ grows.

6.3 Heatmap and line plots (top right and bottom left)

The heatmap gives the exact integer values of $d_{\mathrm{LP}}$ and shows the staircase shape of the feasible region. In the line plot, each color is one value of $k$. The solid curve is the LP bound and the dashed curve of the same color is Singleton. The solid curve lies on or below the dashed one everywhere, and the distance between them is the extra information the LP provides.

6.4 The best code found (bottom middle)

The orange bars are the stabilizer enumerator $A_j$ and the blue bars are the normalizer enumerator $B_j$. For the perfect $[[5,1,3]]$ code we expect

$$
A = (1,0,0,0,15,0), \qquad B = (1,0,0,30,15,18).
$$

The total of $A$ is $16 = 2^{4}$ and the total of $B$ is $64 = 2^{6}$, as required. Because $B_j = A_j$ for $j = 1, 2$ and $B_3 = 30 > 0 = A_3$, there are weight-3 logical operators. This is exactly the statement that the distance is $3$. The console check that the MacWilliams error is at rounding level confirms that these enumerators are exactly the ones the LP constraints describe.

6.5 The distance histogram (bottom right)

Most random $[[5,1]]$ stabilizer codes have a small distance, and only a small fraction reach $d = 3$. The LP bound tells us that $d = 4$ is impossible, so $d = 3$ is optimal. Without that certificate, a random search could never tell whether a better code might still be hiding.

7. Conclusion

Linear programming gives a systematic way to explore the landscape of quantum codes:

The constraints $B_j = A_j$ for $j < d$, $B_j \ge A_j$, and $S_j \ge 0$, together with the quantum MacWilliams identity, define a small LP for each triple $(n,k,d)$.
Infeasibility of this LP proves that no $[[n,k,d]]$ code exists, and the table $d_{\mathrm{LP}}(n,k)$ summarizes all such certificates for $n \le 16$.
The random stabilizer search complements the bound by constructing actual codes, and for $[[5,1]]$ it meets the bound with the perfect five-qubit code.

The same pipeline extends naturally to larger parameters, to non-additive codes by dropping the integrality of $A_j$, and to sharper bounds by adding further constraints to the LP.

Lower-Bounding the Secret Key Rate of Quantum Key Distribution as a Convex Optimization Problem

October 15, 2026

Quantum key distribution (QKD) promises information-theoretic security, but a promise is not a number. To run a real BB84 link we must answer a concrete question: given the error rates Alice and Bob observe, how many secret bits per sifted bit can they safely distill? The answer is a convex optimization problem over everything an eavesdropper could possibly have done. In this article we formulate it, reduce it to a one-dimensional convex program, solve it numerically in Python, and verify the result against the celebrated Shor–Preskill formula.

1. The problem

Consider the entanglement-based picture of BB84. Alice and Bob share an unknown two-qubit state $\rho_{AB}$, and Eve holds a purification of it. Alice and Bob measure in the $Z$ basis to generate the key, and in the $X$ basis to estimate the phase errors. The observed statistics are the bit error rate $e_z$ and the phase error rate $e_x$:

$$
\mathrm{Tr}!\left[(Z\otimes Z),\rho_{AB}\right] = 1-2e_z, \qquad
\mathrm{Tr}!\left[(X\otimes X),\rho_{AB}\right] = 1-2e_x .
$$

With one-way reverse reconciliation, the asymptotic key rate is given by the Devetak–Winter formula

$$
r = H(Z_B,|,E) - f,H(Z_B,|,Z_A) ,
$$

where $f \ge 1$ is the efficiency of the error-correcting code and $H(Z_B|Z_A)=h(e_z)$, with the binary entropy $h(x)=-x\log_2 x-(1-x)\log_2(1-x)$.

The quantity $H(Z_B|E)$ depends on Eve’s attack, which we do not know. We only know that the state must be consistent with the observed statistics. For a purification held by Eve, the conditional entropy can be written as a quantum relative entropy:

where $D(\rho|\sigma)=\mathrm{Tr},\rho(\log_2\rho-\log_2\sigma)$ and $\mathcal{Z}_B$ is the pinching map that measures Bob’s qubit in the $Z$ basis. The worst case over all attacks gives the key rate lower bound, which is the following convex program:

The objective is convex because $\rho\mapsto D(\rho|\mathcal{Z}_B(\rho))$ is the composition of the jointly convex relative entropy with a linear map. All constraints are linear, so any local minimum is a global minimum.

2. Reducing the problem to one variable

A general two-qubit state has 15 real parameters, but symmetry collapses the problem dramatically. Twirling $\rho_{AB}$ over the group ${I\otimes I,\ X\otimes X,\ Y\otimes Y,\ Z\otimes Z}$ leaves both constraints unchanged and maps $\rho_{AB}$ to a Bell-diagonal state. Because $D(\rho|\mathcal{Z}_B(\rho))$ is convex and invariant under these local Pauli unitaries, the twirled state never has a larger objective value. The minimizer can therefore be searched for among Bell-diagonal states:

$$
\rho_{AB}=\lambda_{\Phi^+}|\Phi^+\rangle\langle\Phi^+|+\lambda_{\Phi^-}|\Phi^-\rangle\langle\Phi^-|+\lambda_{\Psi^+}|\Psi^+\rangle\langle\Psi^+|+\lambda_{\Psi^-}|\Psi^-\rangle\langle\Psi^-| .
$$

The two observed error rates translate into linear constraints on the weights:

$$
e_z=\lambda_{\Psi^+}+\lambda_{\Psi^-},\qquad e_x=\lambda_{\Phi^-}+\lambda_{\Psi^-},\qquad \textstyle\sum_i\lambda_i=1 .
$$

Four unknowns and three equations leave a single free parameter $t=\lambda_{\Psi^-}$:

$$
\boldsymbol\lambda(t)=\big(1-e_x-e_z+t,; e_x-t,; e_z-t,; t\big),
\qquad t\in[\max(0,e_x+e_z-1),\ \min(e_x,e_z)] .
$$

This one parameter is Eve’s entire freedom, namely the correlation between the bit-flip and phase-flip errors she induces.

For a Bell-diagonal state the pinching map simply erases coherences, and one finds $S(\mathcal{Z}_B(\rho))=1+h(e_z)$ and $S(\rho)=H(\boldsymbol\lambda)$, where $H$ is the Shannon entropy. Hence

$$
D!\left(\rho|\mathcal{Z}_B(\rho)\right)=1+h(e_z)-H(\boldsymbol\lambda(t)),
\qquad
r(t)=1-H(\boldsymbol\lambda(t))-(f-1),h(e_z).
$$

Minimizing the key rate means maximizing the Shannon entropy of $\boldsymbol\lambda(t)$, which is a concave maximization over an interval. The maximum entropy is attained by the independent distribution $t^\ast=e_xe_z$, giving $H=h(e_x)+h(e_z)$ and the Shor–Preskill bound

$$
r_{\min}=1-h(e_x)-f,h(e_z).
$$

Our numerical solver should rediscover exactly this result without being told.

3. A concrete example

We use the following setup.

Observed bit error rate: $e_z = 0.02$
Observed phase error rate: $e_x = 0.04$
Error-correction efficiency: $f = 1.1$

The program does four things.

It builds the Bell-diagonal state explicitly as a $4\times 4$ matrix, applies the pinching map, and evaluates the quantum relative entropy numerically via eigenvalues.
It solves the convex program with a bounded scalar minimizer and compares the result with the closed-form solution.
It sweeps the $(e_x,e_z)$ plane with a vectorized ternary search, so a $60\times 60$ grid of convex programs is solved in a fraction of a second.
It computes the QBER threshold beyond which no secret key can be extracted.

4. Source code

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import gridspec
from scipy.optimize import minimize_scalar, brentq

plt.style.use("dark_background")
EPS = 1e-15

# ---------------------------------------------------------------
# Basic information-theoretic helpers
# ---------------------------------------------------------------
def h2(x):
    x = np.asarray(x, dtype=float)
    q = np.clip(x, EPS, 1.0 - EPS)
    return -q * np.log2(q) - (1.0 - q) * np.log2(1.0 - q)

def shannon(p):
    # p has shape (4, ...) ; entropy is taken along axis 0
    p = np.asarray(p, dtype=float)
    q = np.clip(p, EPS, 1.0)
    return -np.sum(np.where(p > EPS, p * np.log2(q), 0.0), axis=0)

# ---------------------------------------------------------------
# Eve's one-parameter family of Bell-diagonal states
# order: Phi+, Phi-, Psi+, Psi-
# ---------------------------------------------------------------
def eve_state(t, ex, ez):
    return np.array([1.0 - ex - ez + t, ex - t, ez - t, t])

def eve_vec(t, ex, ez):
    return np.stack([1.0 - ex - ez + t, ex - t, ez - t, t])

# ---------------------------------------------------------------
# Explicit 4x4 matrix construction of the relative entropy
# ---------------------------------------------------------------
s = 1.0 / np.sqrt(2.0)
bell = [
    np.array([s, 0.0, 0.0, s]),    # Phi+
    np.array([s, 0.0, 0.0, -s]),   # Phi-
    np.array([0.0, s, s, 0.0]),    # Psi+
    np.array([0.0, s, -s, 0.0]),   # Psi-
]
PROJ = [
    np.kron(np.eye(2), np.diag([1.0, 0.0])),
    np.kron(np.eye(2), np.diag([0.0, 1.0])),
]

def rho_from(lam):
    return sum(l * np.outer(v, v) for l, v in zip(lam, bell))

def pinch(rho):
    return sum(P @ rho @ P for P in PROJ)

def vn_entropy(rho):
    w = np.linalg.eigvalsh(rho)
    w = w[w > EPS]
    return float(-np.sum(w * np.log2(w)))

def rel_ent(lam):
    rho = rho_from(lam)
    return vn_entropy(pinch(rho)) - vn_entropy(rho)

def key_rate_numeric(t, ex, ez, f):
    return rel_ent(eve_state(t, ex, ez)) - f * float(h2(ez))

# ---------------------------------------------------------------
# Convex program: Eve minimizes the key rate over t
# ---------------------------------------------------------------
def solve_eve(ex, ez, f):
    lo, hi = max(0.0, ex + ez - 1.0), min(ex, ez)
    if hi - lo < 1e-12:
        return lo, key_rate_numeric(lo, ex, ez, f)
    res = minimize_scalar(
        lambda t: key_rate_numeric(t, ex, ez, f),
        bounds=(lo, hi), method="bounded", options={"xatol": 1e-12},
    )
    return float(res.x), float(res.fun)

# Vectorized ternary search: solves a whole grid of convex programs at once
def optimal_t_vec(EX, EZ, iters=80):
    lo = np.maximum(0.0, EX + EZ - 1.0)
    hi = np.minimum(EX, EZ)
    for _ in range(iters):
        m1 = lo + (hi - lo) / 3.0
        m2 = hi - (hi - lo) / 3.0
        f1 = -shannon(eve_vec(m1, EX, EZ))
        f2 = -shannon(eve_vec(m2, EX, EZ))
        left = f1 < f2
        hi = np.where(left, m2, hi)
        lo = np.where(left, lo, m1)
    return 0.5 * (lo + hi)

def rate_grid(EX, EZ, f):
    T = optimal_t_vec(EX, EZ)
    return 1.0 - shannon(eve_vec(T, EX, EZ)) - (f - 1.0) * h2(EZ)

def closed_form(ex, ez, f):
    return 1.0 - h2(ex) - f * h2(ez)

# ---------------------------------------------------------------
# Example problem
# ---------------------------------------------------------------
ex0, ez0, f0 = 0.04, 0.02, 1.1
t_star, r_num = solve_eve(ex0, ez0, f0)
r_closed = float(closed_form(ex0, ez0, f0))
lam_star = eve_state(t_star, ex0, ez0)
check_identity = rel_ent(lam_star) - (1.0 + float(h2(ez0)) - float(shannon(lam_star)))

print("=== Example: e_x = %.3f, e_z = %.3f, f = %.2f ===" % (ex0, ez0, f0))
print("Optimal t (numerical)        : %.8f" % t_star)
print("Optimal t (e_x * e_z)        : %.8f" % (ex0 * ez0))
print("Eve's optimal Bell weights   :", np.round(lam_star, 6))
print("Key rate lower bound (numeric): %.8f" % r_num)
print("Key rate (Shor-Preskill)     : %.8f" % r_closed)
print("Absolute difference          : %.3e" % abs(r_num - r_closed))
print("Matrix vs analytic entropy   : %.3e" % abs(check_identity))

# ---------------------------------------------------------------
# QBER thresholds for symmetric noise e_x = e_z = Q
# ---------------------------------------------------------------
print("\n=== QBER thresholds (e_x = e_z = Q) ===")
thresholds = {}
for fv in (1.0, 1.1, 1.2):
    thr = brentq(lambda q: solve_eve(q, q, fv)[1], 1e-3, 0.3, xtol=1e-10)
    thresholds[fv] = thr
    print("f = %.1f  ->  Q_threshold = %.4f (%.2f %%)" % (fv, thr, 100.0 * thr))

# ---------------------------------------------------------------
# Grid computation over (e_x, e_z)
# ---------------------------------------------------------------
n = 60
exs = np.linspace(0.0005, 0.15, n)
ezs = np.linspace(0.0005, 0.15, n)
EX, EZ = np.meshgrid(exs, ezs)
R = rate_grid(EX, EZ, f0)
err_grid = np.max(np.abs(R - closed_form(EX, EZ, f0)))
print("\nMax |numerical - closed form| over the %dx%d grid: %.3e" % (n, n, err_grid))

# ---------------------------------------------------------------
# Visualization: one figure with six panels
# ---------------------------------------------------------------
fig = plt.figure(figsize=(21, 12))
gs = gridspec.GridSpec(2, 3, figure=fig)
fig.suptitle("QKD Key Rate Lower Bound as a Convex Optimization Problem", fontsize=18, y=0.995)

# (1) Convexity along Eve's free parameter
ax1 = fig.add_subplot(gs[0, 0])
lo0, hi0 = max(0.0, ex0 + ez0 - 1.0), min(ex0, ez0)
tg = np.linspace(lo0, hi0, 300)
rg = np.array([key_rate_numeric(tt, ex0, ez0, f0) for tt in tg])
ax1.plot(tg, rg, color="cyan", lw=2.5, label="key rate $r(t)$")
ax1.scatter([t_star], [r_num], color="red", s=120, zorder=5, label="numerical optimum $t^\\ast$")
ax1.axvline(ex0 * ez0, color="yellow", ls="--", lw=1.5, label="$t=e_x e_z$")
ax1.set_xlabel("Eve's parameter $t=\\lambda_{\\Psi^-}$")
ax1.set_ylabel("key rate  [bits / sifted bit]")
ax1.set_title("Convex objective along Eve's attack family")
ax1.grid(alpha=0.25)
ax1.legend()

# (2) 3D surface of the minimal key rate
ax2 = fig.add_subplot(gs[0, 1], projection="3d")
surf = ax2.plot_surface(EX * 100, EZ * 100, R, cmap="viridis", edgecolor="none", alpha=0.95)
ax2.contourf(EX * 100, EZ * 100, R, zdir="z", offset=float(R.min()), cmap="viridis", alpha=0.6)
ax2.set_xlabel("$e_x$ [%]")
ax2.set_ylabel("$e_z$ [%]")
ax2.set_zlabel("$r_{\\min}$")
ax2.set_title("Key rate lower bound surface ($f=%.1f$)" % f0)
ax2.view_init(elev=28, azim=-125)
fig.colorbar(surf, ax=ax2, shrink=0.55, pad=0.08)

# (3) Secure region map
ax3 = fig.add_subplot(gs[0, 2])
cf = ax3.contourf(EX * 100, EZ * 100, R, levels=30, cmap="viridis")
ax3.contour(EX * 100, EZ * 100, R, levels=[0.0], colors="red", linewidths=2.5)
ax3.scatter([ex0 * 100], [ez0 * 100], color="white", s=90, marker="*", zorder=5)
ax3.set_xlabel("$e_x$ [%]")
ax3.set_ylabel("$e_z$ [%]")
ax3.set_title("Secure region (red line: $r_{\\min}=0$)")
fig.colorbar(cf, ax=ax3)

# (4) Key rate vs symmetric QBER for several code efficiencies
ax4 = fig.add_subplot(gs[1, 0])
Q = np.linspace(0.0005, 0.2, 400)
colors = {1.0: "lime", 1.1: "orange", 1.2: "magenta"}
for fv in (1.0, 1.1, 1.2):
    Tq = optimal_t_vec(Q, Q)
    rq = 1.0 - shannon(eve_vec(Tq, Q, Q)) - (fv - 1.0) * h2(Q)
    ax4.plot(Q * 100, rq, color=colors[fv], lw=2.3, label="$f=%.1f$" % fv)
    ax4.axvline(thresholds[fv] * 100, color=colors[fv], ls="--", alpha=0.7)
ax4.axhline(0.0, color="white", lw=1)
ax4.set_xlabel("QBER $Q=e_x=e_z$ [%]")
ax4.set_ylabel("key rate  [bits / sifted bit]")
ax4.set_title("Key rate and thresholds (dashed)")
ax4.grid(alpha=0.25)
ax4.legend()

# (5) Different attacks, same observed statistics
ax5 = fig.add_subplot(gs[1, 1])
t_vals = [lo0, t_star, hi0]
labels = ["$t_{\\min}$", "$t^\\ast$ (optimal)", "$t_{\\max}$"]
r_vals = [key_rate_numeric(tt, ex0, ez0, f0) for tt in t_vals]
bars = ax5.bar(labels, r_vals, color=["steelblue", "crimson", "seagreen"], edgecolor="white")
for b, v in zip(bars, r_vals):
    ax5.text(b.get_x() + b.get_width() / 2, v, "%.4f" % v, ha="center", va="bottom", fontsize=12)
ax5.set_ylim(min(r_vals) * 0.9, max(r_vals) * 1.03)
ax5.set_ylabel("key rate  [bits / sifted bit]")
ax5.set_title("Three attacks, identical error statistics")
ax5.grid(alpha=0.25, axis="y")

# (6) 3D valley: normalized Eve parameter u vs e_x
ax6 = fig.add_subplot(gs[1, 2], projection="3d")
uu = np.linspace(0.0, 1.0, 50)
exv = np.linspace(0.005, 0.15, 50)
U, EXV = np.meshgrid(uu, exv)
TT = U * np.minimum(EXV, ez0)
RV = 1.0 - shannon(eve_vec(TT, EXV, ez0)) - (f0 - 1.0) * float(h2(ez0))
ax6.plot_surface(U, EXV * 100, RV, cmap="plasma", edgecolor="none", alpha=0.92)
exl = np.linspace(0.005, 0.15, 80)
tl = optimal_t_vec(exl, np.full_like(exl, ez0))
ul = tl / np.minimum(exl, ez0)
rl = 1.0 - shannon(eve_vec(tl, exl, ez0)) - (f0 - 1.0) * float(h2(ez0))
ax6.plot(ul, exl * 100, rl + 0.003, color="cyan", lw=3.5)
ax6.set_xlabel("normalized $u=(t-t_{\\min})/(t_{\\max}-t_{\\min})$")
ax6.set_ylabel("$e_x$ [%]")
ax6.set_zlabel("$r(t)$")
ax6.set_title("Convex valley (cyan: optimal attack, $e_z=%.0f$%%)" % (ez0 * 100))
ax6.view_init(elev=30, azim=-60)

fig.tight_layout()
plt.show()

5. Code walkthrough

Helpers. h2 is the binary entropy and shannon is the Shannon entropy along the first axis, so that one call can process an entire grid of attacks. Both clip their arguments at $10^{-15}$, which makes $0\log 0$ evaluate safely to zero.

Eve’s attack family. eve_state implements $\boldsymbol\lambda(t)=(1-e_x-e_z+t,\ e_x-t,\ e_z-t,\ t)$. Every point of this family reproduces the observed $e_x$ and $e_z$ exactly, and $t$ is the only thing Eve can still tune.

Explicit quantum computation. The four Bell vectors are written in the basis $|00\rangle,|01\rangle,|10\rangle,|11\rangle$. rho_from assembles $\rho_{AB}$, pinch applies $\mathcal{Z}_B$ using the projectors $I\otimes|z\rangle\langle z|$, and vn_entropy computes the von Neumann entropy from the eigenvalues. rel_ent then evaluates $D(\rho|\mathcal{Z}_B(\rho))=S(\mathcal{Z}_B(\rho))-S(\rho)$, with no analytic shortcut involved.

The convex program. solve_eve minimizes $r(t)$ over the feasible interval with a bounded scalar minimizer. Because the objective is convex, the result is the global minimum, so there is no need for multiple restarts.

Speed-up. Solving thousands of convex programs one by one would be wasteful. optimal_t_vec performs a ternary search on all grid points simultaneously with NumPy broadcasting. Ternary search is valid precisely because the objective is convex (unimodal), and 80 iterations shrink the interval by a factor of $(2/3)^{80}\approx 10^{-14}$. The whole $60\times60$ surface is therefore computed with a few hundred array operations.

Self-verification. The program prints the difference between the numerical optimum and the closed-form $1-h(e_x)-f,h(e_z)$ for the example, the difference between the matrix-based relative entropy and its analytic expression $1+h(e_z)-H(\boldsymbol\lambda)$, and the maximum deviation over the entire grid. All three should be at the level of numerical round-off.

Thresholds. For symmetric noise $e_x=e_z=Q$, brentq finds the root of $r_{\min}(Q)=0$ for several values of $f$. For a perfect code ($f=1$) this recovers the famous BB84 limit of about 11 %.

6. Execution results

=== Example: e_x = 0.040, e_z = 0.020, f = 1.10 ===
Optimal t (numerical)        : 0.00080000
Optimal t (e_x * e_z)        : 0.00080000
Eve's optimal Bell weights   : [9.408e-01 3.920e-02 1.920e-02 8.000e-04]
Key rate lower bound (numeric): 0.60212321
Key rate (Shor-Preskill)     : 0.60212321
Absolute difference          : 3.331e-16
Matrix vs analytic entropy   : 1.110e-16

=== QBER thresholds (e_x = e_z = Q) ===
f = 1.0  ->  Q_threshold = 0.1100 (11.00 %)
f = 1.1  ->  Q_threshold = 0.1023 (10.23 %)
f = 1.2  ->  Q_threshold = 0.0955 (9.55 %)

Max |numerical - closed form| over the 60x60 grid: 6.453e-16

7. Interpreting the results

Console output. The optimal attack parameter $t^\ast$ coincides with $e_xe_z$, which means that the worst-case eavesdropper makes the bit-flip and phase-flip errors statistically independent. Independence maximizes the entropy of Eve’s error pattern, and that is exactly what minimizes the information Alice and Bob can extract. The numerical key rate agrees with $1-h(e_x)-f,h(e_z)$ to round-off precision, and the relative-entropy identity check confirms that the matrix computation and the analytic formula describe the same quantity.

Panel 1: the convex objective. Along Eve’s only free parameter, the key rate is a convex curve with a single valley. The red marker (numerical optimum) sits on top of the yellow dashed line ($t=e_xe_z$). Any other value of $t$ corresponds to an attack that is weaker from Eve’s point of view, because it yields a higher key rate for the same observed errors.

Panel 2: the 3D surface. The key rate decreases monotonically in both $e_x$ and $e_z$. The surface is nearly 1 in the corner of vanishing noise, and it slopes downward through zero and into negative values, where no secret key can be distilled. The projected contour on the floor shows the same structure seen from above.

Panel 3: the secure region. The red contour is the boundary $r_{\min}=0$. Inside it, a positive key rate is guaranteed against every possible attack; outside it, key generation must be aborted. The white star marks our example, comfortably inside the secure region. The boundary is symmetric in $e_x$ and $e_z$ only when $f=1$. For $f>1$ the bit error rate is penalized more heavily than the phase error rate, because error correction is applied to the bit errors only.

Panel 4: thresholds. Each curve reaches zero at its threshold (dashed line). A perfect code ($f=1$) tolerates the largest noise, about 11 %, and every additional percent of error-correction inefficiency lowers the tolerable QBER. This is the quantitative reason why good error-correcting codes matter for practical QKD systems.

Panel 5: different attacks, same statistics. The three bars correspond to three attacks that all produce identical $e_x$ and $e_z$ and therefore look indistinguishable to Alice and Bob. Yet they give different key rates. The red bar, the independent-error attack, is the lowest. This illustrates why a security proof must take the minimum over all attacks consistent with the data rather than trust a single plausible noise model.

Panel 6: the 3D valley. Here the normalized attack parameter $u$ runs along one axis and the phase error rate along the other. For every fixed $e_x$ the surface forms a convex valley in $u$, and the cyan curve traces the bottom of this valley, which is the optimal attack as $e_x$ varies. The valley is smooth and has no secondary minima, which is the geometric face of convexity and the reason why a simple ternary search is guaranteed to find the true worst case.

8. Conclusion

The security of BB84 reduces to a convex program: minimize a quantum relative entropy over all states compatible with the observed error rates. Symmetry reduces this program from fifteen parameters to a single one, and the solution is the maximum-entropy attack in which bit and phase errors are independent, giving the Shor–Preskill rate $1-h(e_x)-f,h(e_z)$. The same framework, with a numerical convex solver replacing the closed form, extends to protocols and channels where no analytic solution exists, such as asymmetric protocols, the six-state protocol, and finite-key or imperfect-source scenarios.

Quantum Values of Nonlocal Games

October 14, 2026

How Much Better Can Entanglement Do? A Worked CHSH Example in Python

Two players, Alice and Bob, are placed in separate rooms and cannot communicate. A referee sends each of them a question, and each returns an answer. Together they win or lose according to a fixed rule. If they are allowed to share an entangled quantum state, how high can their winning probability go?

That number is the quantum value of the nonlocal game:

$$
\omega^*(G)=\sup_{|\psi\rangle,;{M^x_a},;{N^y_b}};\sum_{x,y}\pi(x,y)\sum_{a,b}V(a,b\mid x,y),\langle\psi|,M^x_a\otimes N^y_b,|\psi\rangle
$$

Here $\pi(x,y)$ is the question distribution, $V(a,b\mid x,y)\in{0,1}$ is the winning predicate, and the supremum runs over all shared states and all local measurements.

In this post we compute it for the most famous example, the CHSH game, in three independent ways: an analytic strategy, a numerical see-saw optimisation, and a Monte Carlo simulation of the game itself.

1. The CHSH game

The referee draws two uniformly random bits $x,y\in{0,1}$ and sends $x$ to Alice and $y$ to Bob. They reply with bits $a$ and $b$. They win if and only if

$$
a\oplus b ;=; x\wedge y .
$$

In words: the answers must be equal unless both questions are $1$, in which case they must differ.

Classical value. Any shared randomness is a mixture of deterministic strategies, so it suffices to check the $2^4=16$ deterministic ones. Winning all four question pairs would require $a_0\oplus b_0=0$, $a_0\oplus b_1=0$, $a_1\oplus b_0=0$ and $a_1\oplus b_1=1$. Adding the four equations modulo 2 gives $0=1$, which is a contradiction. At most three of the four can hold, so

$$
\omega_{\text{classical}}=\frac34 .
$$

Quantum value. With $\pm1$-valued observables $A_x,B_y$, the winning probability becomes

$$
P_{\text{win}}=\frac14\sum_{x,y}\frac{1+(-1)^{xy}\langle A_x\otimes B_y\rangle}{2}
=\frac12+\frac{S}{8},
$$

where

$$
S=\langle A_0B_0\rangle+\langle A_0B_1\rangle+\langle A_1B_0\rangle-\langle A_1B_1\rangle .
$$

Classically $|S|\le 2$, which reproduces $3/4$. Quantum mechanically, Tsirelson’s bound gives $|S|\le 2\sqrt2$, so

$$
\omega^*(\mathrm{CHSH})=\frac12+\frac{2\sqrt2}{8}=\cos^2!\left(\frac{\pi}{8}\right)\approx 0.853553 .
$$

The no-signalling value is $1$, so the game separates three regimes: $0.75<0.8536<1$.

An explicit optimal strategy

Take the Bell state $|\Phi^+\rangle=\frac{1}{\sqrt2}(|00\rangle+|11\rangle)$ and the observables

$$
A(\theta)=\cos\theta,Z+\sin\theta,X .
$$

On $|\Phi^+\rangle$ we have $\langle Z\otimes Z\rangle=\langle X\otimes X\rangle=1$ and $\langle Z\otimes X\rangle=\langle X\otimes Z\rangle=0$, hence

$$
\langle\Phi^+|A(\alpha)\otimes A(\beta)|\Phi^+\rangle=\cos(\alpha-\beta).
$$

Choosing $\alpha_0=0,\ \alpha_1=\pi/2,\ \beta_0=\pi/4,\ \beta_1=-\pi/4$ gives each of the four correlators magnitude $\tfrac{\sqrt2}{2}$ with exactly the right signs, so $S=2\sqrt2$.

Without assuming the strategy: the see-saw method

The explicit strategy shows that $\cos^2(\pi/8)$ is achievable. To check that a generic optimiser also finds nothing better, we use a see-saw iteration in local dimension $d$. Define the Bell operator

$$
W=\sum_{x,y}(-1)^{xy},A_x\otimes B_y .
$$

The iteration repeats three exactly solvable steps:

State step. For fixed observables, the best state is the eigenvector of $W$ with the largest eigenvalue $\lambda_{\max}$.
Alice step. For a fixed state and fixed $B_y$, maximise $\operatorname{Tr}(A_xM_x)$ over $\pm1$-valued observables, where $M_x$ is obtained by partially tracing out Bob. The solution is $A_x=\operatorname{sign}!\big(\tfrac{M_x+M_x^\dagger}{2}\big)$.
Bob step. The same as above with the roles exchanged.

Each step can only increase the value, so the winning probability $\tfrac12+\tfrac{\lambda}{8}$ rises monotonically and converges to a local optimum.

2. The complete source code

Everything is in a single cell.

import itertools
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm

np.set_printoptions(precision=6, suppress=True)
rng = np.random.default_rng(2026)

# ------------------------------------------------------------
# 1. CHSH game: win iff a XOR b == x AND y, inputs uniform
# ------------------------------------------------------------
def win_pred(a, b, x, y):
    return (a ^ b) == (x & y)

# ------------------------------------------------------------
# 2. Classical value by brute force over deterministic strategies
# ------------------------------------------------------------
def classical_value():
    best = 0.0
    for a0, a1, b0, b1 in itertools.product([0, 1], repeat=4):
        A = (a0, a1)
        B = (b0, b1)
        wins = sum(win_pred(A[x], B[y], x, y) for x in (0, 1) for y in (0, 1))
        best = max(best, wins / 4)
    return best

# ------------------------------------------------------------
# 3. Quantum strategy with Bell state and rotated observables
# ------------------------------------------------------------
I2 = np.eye(2)
Z = np.array([[1, 0], [0, -1]], dtype=complex)
X = np.array([[0, 1], [1, 0]], dtype=complex)
phi_plus = np.array([1, 0, 0, 1], dtype=complex) / np.sqrt(2)

def obs(theta):
    return np.cos(theta) * Z + np.sin(theta) * X

def win_probability(alpha, beta):
    total = 0.0
    for x in (0, 1):
        for y in (0, 1):
            corr = np.real(phi_plus.conj() @ np.kron(obs(alpha[x]), obs(beta[y])) @ phi_plus)
            total += 0.25 * (1 + (-1) ** (x * y) * corr) / 2
    return total

def win_probability_fast(a0, a1, b0, b1):
    # <A(a) x B(b)> = cos(a - b) on the Bell state (vectorised version)
    s = np.cos(a0 - b0) + np.cos(a0 - b1) + np.cos(a1 - b0) - np.cos(a1 - b1)
    return 0.5 + s / 8

alpha_opt = (0.0, np.pi / 2)
beta_opt = (np.pi / 4, -np.pi / 4)

# ------------------------------------------------------------
# 4. See-saw optimisation in dimension d (no assumption on the strategy)
# ------------------------------------------------------------
def random_state(dim, rng):
    v = rng.normal(size=dim) + 1j * rng.normal(size=dim)
    return v / np.linalg.norm(v)

def sign_operator(M):
    w, V = np.linalg.eigh((M + M.conj().T) / 2)
    s = np.where(w >= 0, 1.0, -1.0)
    return (V * s) @ V.conj().T

def random_observable(d, rng):
    return sign_operator(rng.normal(size=(d, d)) + 1j * rng.normal(size=(d, d)))

def seesaw(d, iters=60, rng=rng):
    A = [random_observable(d, rng) for _ in range(2)]
    B = [random_observable(d, rng) for _ in range(2)]
    history = []
    for _ in range(iters):
        # (i) optimal state = top eigenvector of the Bell operator
        W = sum((-1) ** (x * y) * np.kron(A[x], B[y]) for x in (0, 1) for y in (0, 1))
        w, V = np.linalg.eigh(W)
        psi = V[:, -1]
        rho = np.outer(psi, psi.conj()).reshape(d, d, d, d)  # indices a b a' b'
        # (ii) optimal Alice observables for fixed state and Bob
        for x in (0, 1):
            Mx = sum((-1) ** (x * y) * np.einsum('abcd,db->ac', rho, B[y]) for y in (0, 1))
            A[x] = sign_operator(Mx)
        # (iii) optimal Bob observables for fixed state and Alice
        for y in (0, 1):
            My = sum((-1) ** (x * y) * np.einsum('abcd,ca->bd', rho, A[x]) for x in (0, 1))
            B[y] = sign_operator(My)
        W = sum((-1) ** (x * y) * np.kron(A[x], B[y]) for x in (0, 1) for y in (0, 1))
        value = np.real(psi.conj() @ W @ psi)
        history.append(0.5 + value / 8)
    return np.array(history)

# ------------------------------------------------------------
# 5. Monte Carlo simulation of the actual game
# ------------------------------------------------------------
def simulate_game(n_rounds, rng):
    wins = np.zeros(n_rounds, dtype=float)
    xs = rng.integers(0, 2, n_rounds)
    ys = rng.integers(0, 2, n_rounds)
    # exact outcome distributions for each input pair (Born rule)
    probs = {}
    for x in (0, 1):
        for y in (0, 1):
            P = np.zeros(4)
            for a in (0, 1):
                for b in (0, 1):
                    Pa = (I2 + (-1) ** a * obs(alpha_opt[x])) / 2
                    Pb = (I2 + (-1) ** b * obs(beta_opt[y])) / 2
                    P[2 * a + b] = np.real(phi_plus.conj() @ np.kron(Pa, Pb) @ phi_plus)
            probs[(x, y)] = P / P.sum()
    for x in (0, 1):
        for y in (0, 1):
            idx = np.where((xs == x) & (ys == y))[0]
            outcome = rng.choice(4, size=len(idx), p=probs[(x, y)])
            a = outcome // 2
            b = outcome % 2
            wins[idx] = ((a ^ b) == (x & y)).astype(float)
    return np.cumsum(wins) / np.arange(1, n_rounds + 1)

# ------------------------------------------------------------
# 6. Run everything
# ------------------------------------------------------------
c_val = classical_value()
q_val = win_probability(alpha_opt, beta_opt)
tsirelson = np.cos(np.pi / 8) ** 2

print("Classical value            :", c_val)
print("Quantum value (explicit)   :", round(q_val, 10))
print("cos^2(pi/8) (analytic)     :", round(tsirelson, 10))
print("No-signalling value        :", 1.0)

hist2 = [seesaw(2, iters=40) for _ in range(8)]
hist3 = [seesaw(3, iters=40) for _ in range(8)]
print("See-saw best (d=2)         :", round(max(h[-1] for h in hist2), 8))
print("See-saw best (d=3)         :", round(max(h[-1] for h in hist3), 8))

n_rounds = 200000
running = simulate_game(n_rounds, rng)
print("Monte Carlo win rate       :", round(running[-1], 5), f"({n_rounds} rounds)")

# ------------------------------------------------------------
# 7. One combined figure
# ------------------------------------------------------------
fig = plt.figure(figsize=(20, 12))

ax1 = fig.add_subplot(2, 3, 1)
labels = ["Classical", "Quantum", "No-signalling"]
vals = [c_val, tsirelson, 1.0]
bars = ax1.bar(labels, vals, color=["#4c72b0", "#dd8452", "#8c8c8c"])
for b, v in zip(bars, vals):
    ax1.text(b.get_x() + b.get_width() / 2, v + 0.01, f"{v:.4f}", ha="center", fontsize=12)
ax1.set_ylim(0.5, 1.05)
ax1.set_ylabel("Maximum winning probability")
ax1.set_title("CHSH game: three kinds of value")

b0 = np.linspace(-np.pi, np.pi, 160)
b1 = np.linspace(-np.pi, np.pi, 160)
B0, B1 = np.meshgrid(b0, b1)
P = win_probability_fast(0.0, np.pi / 2, B0, B1)

ax2 = fig.add_subplot(2, 3, 2, projection="3d")
surf = ax2.plot_surface(B0, B1, P, cmap=cm.viridis, linewidth=0, antialiased=True, alpha=0.95)
ax2.scatter([np.pi / 4], [-np.pi / 4], [tsirelson], color="red", s=80, label="optimum")
ax2.set_xlabel(r"$\beta_0$")
ax2.set_ylabel(r"$\beta_1$")
ax2.set_zlabel("Win probability")
ax2.set_title("3D landscape of Bob's angles")
ax2.legend(loc="upper left")
fig.colorbar(surf, ax=ax2, shrink=0.6, pad=0.1)

ax3 = fig.add_subplot(2, 3, 3)
cs = ax3.contourf(B0, B1, P, levels=30, cmap="viridis")
ax3.plot([np.pi / 4], [-np.pi / 4], "r*", markersize=18)
ax3.set_xlabel(r"$\beta_0$")
ax3.set_ylabel(r"$\beta_1$")
ax3.set_title("Contour view (red star = optimum)")
fig.colorbar(cs, ax=ax3)

ax4 = fig.add_subplot(2, 3, 4)
for i, h in enumerate(hist2):
    ax4.plot(h, color="#dd8452", alpha=0.7, label="d = 2" if i == 0 else None)
for i, h in enumerate(hist3):
    ax4.plot(h, color="#55a868", alpha=0.7, label="d = 3" if i == 0 else None)
ax4.axhline(tsirelson, color="red", linestyle="--", label=r"$\cos^2(\pi/8)$")
ax4.axhline(0.75, color="gray", linestyle=":", label="classical 3/4")
ax4.set_xlabel("See-saw iteration")
ax4.set_ylabel("Winning probability")
ax4.set_title("See-saw optimisation from random starts")
ax4.legend()

ax5 = fig.add_subplot(2, 3, 5)
ax5.plot(np.arange(1, n_rounds + 1), running, color="#4c72b0", linewidth=1)
ax5.axhline(tsirelson, color="red", linestyle="--", label=r"$\cos^2(\pi/8)\approx 0.8536$")
ax5.axhline(0.75, color="gray", linestyle=":", label="classical 3/4")
ax5.set_xscale("log")
ax5.set_xlabel("Number of rounds")
ax5.set_ylabel("Running win rate")
ax5.set_title("Monte Carlo simulation of the game")
ax5.legend()

ax6 = fig.add_subplot(2, 3, 6)
theta = np.linspace(0, 2 * np.pi, 400)
ax6.plot(np.cos(theta), np.sin(theta), color="lightgray")
colors = {"A0": "#4c72b0", "A1": "#55a868", "B0": "#dd8452", "B1": "#c44e52"}
dirs = {"A0": alpha_opt[0], "A1": alpha_opt[1], "B0": beta_opt[0], "B1": beta_opt[1]}
for name, ang in dirs.items():
    ax6.annotate("", xy=(np.cos(ang), np.sin(ang)), xytext=(0, 0),
                 arrowprops=dict(arrowstyle="->", color=colors[name], lw=2.5))
    ax6.text(1.12 * np.cos(ang), 1.12 * np.sin(ang), name, color=colors[name],
             fontsize=13, ha="center", va="center")
ax6.set_xlim(-1.4, 1.4)
ax6.set_ylim(-1.4, 1.4)
ax6.set_aspect("equal")
ax6.set_title("Optimal measurement directions (Bloch x-z plane)")

plt.tight_layout()
plt.show()

3. Code walkthrough

Section 1–2: the game and the classical value

win_pred(a, b, x, y) is the winning predicate $a\oplus b=x\wedge y$ written with Python’s bitwise operators. classical_value() enumerates all $2^4=16$ deterministic strategies $(a_0,a_1,b_0,b_1)$, counts how many of the four question pairs each one wins, and keeps the best fraction. Deterministic strategies are enough because any classical strategy using shared randomness is a convex combination of them. The result is exactly $3/4$.

Section 3: the explicit quantum strategy

obs(theta) builds the $\pm1$-valued observable $A(\theta)=\cos\theta,Z+\sin\theta,X$, a measurement along a direction in the $x$–$z$ plane of the Bloch sphere.
phi_plus is the Bell state $|\Phi^+\rangle$ as a four-component vector.
win_probability evaluates $\tfrac14\sum_{x,y}\tfrac{1+(-1)^{xy}\langle A_x\otimes B_y\rangle}{2}$ literally, using np.kron for the tensor product. It serves as the reference implementation.
win_probability_fast uses the closed form $\langle A(\alpha)\otimes A(\beta)\rangle=\cos(\alpha-\beta)$, so the entire expression is a sum of four cosines. It accepts NumPy arrays, which matters in Section 7.

alpha_opt and beta_opt store the optimal angles derived above.

Section 4: the see-saw

sign_operator(M) takes the Hermitian part of $M$, diagonalises it, and replaces each eigenvalue by $+1$ or $-1$ according to its sign. This is the matrix sign function and gives the best $\pm1$-valued observable for a linear objective $\operatorname{Tr}(AM)$.
random_observable(d, rng) produces a random starting point by applying sign_operator to a random complex matrix.
Inside seesaw, the state is reshaped into a four-index tensor rho[a, b, a', b']. The einsum calls contract this tensor with Bob’s observable (or Alice’s) to give the effective operator $M_x$ (or $M_y$) acting on the other party’s system. The result is passed to sign_operator.
After each full sweep, the Bell operator’s expectation value is converted to a winning probability $\tfrac12+\tfrac{\lambda}{8}$ and appended to history.

The state is found as the top eigenvector of a $d^2\times d^2$ matrix, so for $d=2,3$ the whole optimisation takes a fraction of a second. Using local dimension $d=3$ lets the optimiser use strategies beyond qubits. The fact that it does not exceed $\cos^2(\pi/8)$ is consistent with Tsirelson’s theorem.

Section 5: Monte Carlo simulation

This section plays the game rather than computing an expectation. First, for each of the four question pairs, the exact joint outcome distribution $p(a,b\mid x,y)=\langle\Phi^+|P^x_a\otimes P^y_b|\Phi^+\rangle$ is computed from the projectors $P=\tfrac{I\pm A}{2}$. Then 200,000 question pairs are sampled.

Speed-up. A naive implementation would loop over rounds and sample one outcome per round, which is slow in pure Python. Here the rounds are grouped by question pair, and all the outcomes for each group are drawn with a single vectorised rng.choice call. The final running average is computed with np.cumsum. Only four sampling calls are needed instead of 200,000.

Section 6: running everything

The first block prints the classical, explicit quantum, analytic and no-signalling values. The see-saw is then run from 8 random starts in each of $d=2$ and $d=3$, and the best final values are printed. Finally the Monte Carlo game is run and its win rate is printed. Random starts are used because see-saw only guarantees a local optimum, and several restarts guard against a poor one.

Section 7: the figure

Speed-up for the 3D surface. Evaluating the winning probability on a $160\times160$ grid with win_probability would require 25,600 calls with four np.kron multiplications each. Instead, win_probability_fast accepts the whole meshgrid and returns the full surface in one vectorised expression.

4. Results as graphs

The six panels are in one figure.

How to read each panel

Panel 1, three kinds of value. The bars show $0.75$ (classical), $\cos^2(\pi/8)\approx0.8536$ (quantum) and $1$ (no-signalling). Entanglement gives a real, measurable advantage, but it stops well short of what an unrestricted no-signalling resource would allow.

Panel 2, 3D landscape. With Alice fixed at $\alpha_0=0,\ \alpha_1=\pi/2$, the surface shows Bob’s winning probability as a function of his two angles $(\beta_0,\beta_1)$. It is a smooth, periodic landscape with several peaks, and the red point marks the global maximum at $(\pi/4,-\pi/4)$ with height $\cos^2(\pi/8)$. The surface never rises above that point, which is Tsirelson’s bound in geometric form.

Panel 3, contour view. The same surface seen from above. The red star sits at the center of the brightest region, and the symmetric pattern of ridges reflects the cosine structure $\cos(\alpha-\beta)$.

Panel 4, see-saw convergence. Each curve is one random start. Every curve increases monotonically, as the theory guarantees, and the curves for both $d=2$ and $d=3$ settle at the red dashed line $\cos^2(\pi/8)$ and never cross it. Some starts get there in a few iterations; others pause on a plateau first. This is numerical evidence that the explicit qubit strategy is optimal even when more general strategies are allowed.

Panel 5, Monte Carlo. The running win rate of the simulated game fluctuates strongly for the first few hundred rounds and then converges toward the red dashed line. The log-scaled horizontal axis makes the $1/\sqrt{n}$ shrinkage of the fluctuations visible. The long-run rate clearly lies above the gray classical line at $3/4$, so no classical strategy could explain it.

Panel 6, optimal directions. The four arrows show the measurement axes of the optimal strategy in the Bloch $x$–$z$ plane. Alice’s axes ($A_0,A_1$) are orthogonal, and Bob’s axes ($B_0,B_1$) are orthogonal to each other and lie exactly halfway between Alice’s two axes. This symmetric picture is the geometric heart of the strategy: each of Bob’s axes is only $45^\circ$ away from both of Alice’s, which makes three of the four correlations as strong as possible, while the fourth is tuned to carry the required minus sign.

5. Console output

Classical value            : 0.75
Quantum value (explicit)   : 0.8535533906
cos^2(pi/8) (analytic)     : 0.8535533906
No-signalling value        : 1.0
See-saw best (d=2)         : 0.85355339
See-saw best (d=3)         : 0.85355339
Monte Carlo win rate       : 0.85456 (200000 rounds)

How to read the numbers

Classical value is exactly $3/4$, in agreement with the argument that no deterministic strategy wins all four question pairs.
Quantum value (explicit) and $\cos^2(\pi/8)$ agree to ten digits, confirming the analytic derivation.
See-saw best for $d=2$ and $d=3$ both match $\cos^2(\pi/8)$ to the printed precision. Higher local dimension did not help.
Monte Carlo win rate is close to $0.8536$. With $n=200{,}000$ rounds the statistical standard deviation is about $\sqrt{p(1-p)/n}\approx 0.0008$, so deviations in the third decimal place are expected.

6. Wrap-up

The quantum value of a nonlocal game is the best winning probability over all entangled strategies. For the CHSH game it equals

$$
\omega^*(\mathrm{CHSH})=\cos^2!\left(\frac{\pi}{8}\right)=\frac12+\frac{\sqrt2}{4}\approx 0.8536 ,
$$

compared with $3/4$ classically. The explicit Bell-state strategy reaches this value, the see-saw optimisation finds nothing better, and a simulation of the game confirms it empirically. For general games the exact quantum value is hard to compute, but see-saw gives lower bounds, and semidefinite hierarchies such as the NPA hierarchy give matching upper bounds. Try changing the winning predicate or the question distribution in the code and see how the three values move.

Verifying Quantum Steering with a Semidefinite Program

October 13, 2026

Werner States, Steerable Weight, and Critical Visibility

Quantum steering sits between entanglement and Bell nonlocality. Alice measures her half of a shared state, and she can steer Bob’s conditional states in a way that no local hidden state model can reproduce. In this post we decide steerability numerically with a semidefinite program (SDP), using a two-qubit Werner state as a concrete example.

1. The problem

Alice and Bob share a state $\rho_{AB}$. Alice chooses a measurement setting $x \in {1,\dots,m}$ and obtains an outcome $a \in {0,1}$, described by projectors $P_{a|x}$. Bob’s unnormalized conditional states form the assemblage

The assemblage is unsteerable if it admits a local hidden state (LHS) model:

$$
\sigma_{a|x} = \sum_{\lambda} D(a|x,\lambda),\sigma_\lambda,
\qquad \sigma_\lambda \succeq 0 ,
$$

where $\lambda = (\lambda_1,\dots,\lambda_m)$ runs over all $2^m$ deterministic strategies and $D(a|x,\lambda)=\delta_{a,\lambda_x}$. Since there are finitely many $\lambda$ and every constraint is linear in $\sigma_\lambda$, this is a semidefinite feasibility problem.

Steerable weight

Instead of a yes/no answer, we compute a quantitative measure, the steerable weight. We split the assemblage into an unsteerable part and a remainder, and maximize the weight of the unsteerable part:

$$
\begin{aligned}
\mu^\star = \max_ \quad & \sum_\lambda \mathrm{Tr},\sigma_\lambda \
\text{s.t.} \quad & \sum_{\lambda}D(a|x,\lambda),\sigma_\lambda \preceq \sigma_{a|x}\quad \forall a,x,\
& \sigma_\lambda \succeq 0 ,
\end{aligned}
\qquad
\mathrm{SW} = 1-\mu^\star .
$$

The assemblage is unsteerable exactly when $\mathrm{SW}=0$.

Critical visibility in a single SDP

Take the Werner state

$$
\rho(p) = p,|\psi^-\rangle\langle\psi^-| + (1-p),\frac{\mathbb{1}}{4},
\qquad |\psi^-\rangle=\frac{|01\rangle-|10\rangle}{\sqrt2}.
$$

For a projective measurement along the unit vector $\vec n$, the assemblage is

$$
\sigma_{a|x}(p) = \frac14\Big(\mathbb{1} - p,(-1)^a,\vec n_x\cdot\vec\sigma\Big)
= p,\sigma^{\rm singlet}_{a|x} + (1-p)\frac{\mathbb{1}}{4}.
$$

Bob’s conditional Bloch vectors therefore have length $p$. Since $p$ enters the constraints linearly, we can make it an SDP variable and find the critical visibility $p_c$ in one shot:

$$
p_c = \max_{p,{\sigma_\lambda}}\ p
\quad\text{s.t.}\quad
\sum_{\lambda:,\lambda_x=a}\sigma_\lambda = p,\sigma^{\rm singlet}_{a|x}+(1-p)\frac{\mathbb{1}}{4},\ \ \sigma_\lambda\succeq0 .
$$

For $p\le p_c$ the state is unsteerable with the chosen settings, and for $p>p_c$ it is steerable.

2. The worked example

We study three things:

Critical visibility versus number of settings. We use $m=2$ (the $x,z$ axes), $m=3$ (the $x,y,z$ axes), $m=4$ (four cube diagonals), and $m=6$ (six icosahedron axes).
Steerable weight curves $\mathrm{SW}(p)$ for $m=2,3,4$.
Two settings at relative angle $\theta$, which gives a 3D surface $\mathrm{SW}(p,\theta)$. The critical curve is expected to be

$$
p_c(\theta)=\frac{1}{\sqrt{1+\sin\theta}} .
$$

Known reference values are $p_c=1/\sqrt2\approx0.7071$ for two orthogonal settings, $p_c=1/\sqrt3\approx0.5774$ for three orthogonal settings, and $p_c\to 1/2$ in the limit of infinitely many projective settings.

3. Source code

import itertools
import time

import numpy as np
import cvxpy as cp
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D  # noqa: F401

# ---------------------------------------------------------------
# 1. Basic objects: Pauli matrices, Werner state, measurements
# ---------------------------------------------------------------
I2 = np.eye(2, dtype=complex)
SX = np.array([[0, 1], [1, 0]], dtype=complex)
SY = np.array([[0, -1j], [1j, 0]], dtype=complex)
SZ = np.array([[1, 0], [0, -1]], dtype=complex)

PSI = np.array([0, 1, -1, 0], dtype=complex) / np.sqrt(2)
RHO_SINGLET = np.outer(PSI, PSI.conj())


def werner_state(p):
    """rho(p) = p |psi-><psi-| + (1-p) I/4"""
    return p * RHO_SINGLET + (1 - p) * np.eye(4, dtype=complex) / 4


def projector(n, a):
    """Projector onto outcome a (0 or 1) of the spin measurement along n."""
    n = np.asarray(n, dtype=float)
    n = n / np.linalg.norm(n)
    return 0.5 * (I2 + (-1) ** a * (n[0] * SX + n[1] * SY + n[2] * SZ))


def assemblage(rho, directions):
    """sigma_{a|x} = Tr_A[(P_{a|x} (x) I) rho]  ->  list[x][a] of 2x2 matrices"""
    out = []
    for n in directions:
        row = []
        for a in (0, 1):
            M = np.kron(projector(n, a), I2)
            full = (M @ rho).reshape(2, 2, 2, 2)
            row.append(np.einsum("ijik->jk", full))
        out.append(row)
    return out


# ---------------------------------------------------------------
# 2. SDP solver helper
# ---------------------------------------------------------------
if "CLARABEL" in cp.installed_solvers():
    SOLVE_KW = dict(solver=cp.CLARABEL)
else:
    SOLVE_KW = dict(solver=cp.SCS, eps=1e-9, max_iters=200000)


# ---------------------------------------------------------------
# 3. SDP (A): steerable weight for a given assemblage
# ---------------------------------------------------------------
def steerable_weight(sigma):
    m = len(sigma)
    strategies = list(itertools.product((0, 1), repeat=m))
    S = [cp.Variable((2, 2), hermitian=True) for _ in strategies]
    cons = [s >> 0 for s in S]
    for x in range(m):
        for a in (0, 1):
            part = sum(S[k] for k, lam in enumerate(strategies) if lam[x] == a)
            cons.append(sigma[x][a] - part >> 0)
    objective = cp.Maximize(cp.real(sum(cp.trace(s) for s in S)))
    prob = cp.Problem(objective, cons)
    prob.solve(**SOLVE_KW)
    return max(0.0, 1.0 - float(prob.value))


# ---------------------------------------------------------------
# 4. SDP (B), slow: critical visibility by bisection
# ---------------------------------------------------------------
def critical_p_bisection(directions, lo=0.3, hi=1.0, iters=22, tol=1e-6):
    for _ in range(iters):
        mid = 0.5 * (lo + hi)
        sw = steerable_weight(assemblage(werner_state(mid), directions))
        if sw > tol:
            hi = mid
        else:
            lo = mid
    return 0.5 * (lo + hi)


# ---------------------------------------------------------------
# 5. SDP (C), fast: critical visibility with ONE SDP
#    p becomes a variable:  sum_{lam: lam_x=a} sigma_lam = p*sigma^singlet_{a|x} + (1-p) I/4
# ---------------------------------------------------------------
def critical_p_direct(directions):
    m = len(directions)
    sig_s = assemblage(RHO_SINGLET, directions)
    strategies = list(itertools.product((0, 1), repeat=m))
    S = [cp.Variable((2, 2), hermitian=True) for _ in strategies]
    p = cp.Variable()
    cons = [s >> 0 for s in S] + [p >= 0, p <= 1]
    for x in range(m):
        for a in (0, 1):
            part = sum(S[k] for k, lam in enumerate(strategies) if lam[x] == a)
            cons.append(part == p * sig_s[x][a] + (1 - p) * I2 / 4)
    prob = cp.Problem(cp.Maximize(p), cons)
    prob.solve(**SOLVE_KW)
    return float(p.value)


# ---------------------------------------------------------------
# 6. Measurement settings
# ---------------------------------------------------------------
PHI = (1 + np.sqrt(5)) / 2
SETTINGS = {
    2: [(1, 0, 0), (0, 0, 1)],
    3: [(1, 0, 0), (0, 1, 0), (0, 0, 1)],
    4: [(1, 1, 1), (1, -1, -1), (-1, 1, -1), (-1, -1, 1)],
    6: [(0, 1, PHI), (0, 1, -PHI), (1, PHI, 0), (1, -PHI, 0), (PHI, 0, 1), (-PHI, 0, 1)],
}

# ---------------------------------------------------------------
# 7. Experiment 1: critical visibility vs. number of settings
# ---------------------------------------------------------------
print("=" * 64)
print("Experiment 1: critical visibility p_c of the Werner state")
print("=" * 64)
print(f"{'settings m':>10} {'strategies 2^m':>15} {'p_c (SDP)':>12} {'time [s]':>10}")
pc = {}
for m, dirs in SETTINGS.items():
    t0 = time.time()
    pc[m] = critical_p_direct(dirs)
    print(f"{m:>10d} {2**m:>15d} {pc[m]:>12.6f} {time.time() - t0:>10.3f}")
print(f"\nReference: 1/sqrt(2) = {1/np.sqrt(2):.6f},  1/sqrt(3) = {1/np.sqrt(3):.6f},  1/2 = 0.500000")

# ---------------------------------------------------------------
# 8. Speed comparison: bisection vs. one-shot SDP (m = 4)
# ---------------------------------------------------------------
print("\n" + "=" * 64)
print("Speed comparison (m = 4 settings)")
print("=" * 64)
t0 = time.time()
pc_bis = critical_p_bisection(SETTINGS[4])
t_bis = time.time() - t0
t0 = time.time()
pc_dir = critical_p_direct(SETTINGS[4])
t_dir = time.time() - t0
print(f"bisection (22 SDPs): p_c = {pc_bis:.6f}, time = {t_bis:.3f} s")
print(f"one-shot SDP       : p_c = {pc_dir:.6f}, time = {t_dir:.3f} s")
print(f"speed-up           : x{t_bis / max(t_dir, 1e-9):.1f}")

# ---------------------------------------------------------------
# 9. Experiment 2: steerable weight curves
# ---------------------------------------------------------------
p_grid = np.linspace(0.4, 1.0, 31)
curve_ms = [2, 3, 4]
sw_curves = {}
for m in curve_ms:
    sw_curves[m] = np.array(
        [steerable_weight(assemblage(werner_state(p), SETTINGS[m])) for p in p_grid]
    )

print("\n" + "=" * 64)
print("Experiment 2: steerable weight SW(p)")
print("=" * 64)
print(f"{'p':>6} " + " ".join(f"{'SW (m=' + str(m) + ')':>12}" for m in curve_ms))
for i in range(0, len(p_grid), 3):
    print(f"{p_grid[i]:>6.2f} " + " ".join(f"{sw_curves[m][i]:>12.6f}" for m in curve_ms))

# ---------------------------------------------------------------
# 10. Experiment 3: two settings with relative angle theta
# ---------------------------------------------------------------
def two_settings(theta):
    return [(0, 0, 1), (np.sin(theta), 0, np.cos(theta))]


thetas = np.linspace(0.05, np.pi / 2, 18)
ps = np.linspace(0.5, 1.0, 21)
SW_surf = np.zeros((len(thetas), len(ps)))
for i, th in enumerate(thetas):
    for j, p in enumerate(ps):
        SW_surf[i, j] = steerable_weight(assemblage(werner_state(p), two_settings(th)))

th_line = np.linspace(0.05, np.pi / 2, 40)
pc_num = np.array([critical_p_direct(two_settings(th)) for th in th_line])
pc_formula = 1 / np.sqrt(1 + np.sin(th_line))

print("\n" + "=" * 64)
print("Experiment 3: two settings at angle theta")
print("=" * 64)
print(f"{'theta [deg]':>12} {'p_c (SDP)':>12} {'1/sqrt(1+sin)':>15}")
for th in np.radians([15, 30, 45, 60, 75, 90]):
    pn = critical_p_direct(two_settings(th))
    print(f"{np.degrees(th):>12.1f} {pn:>12.6f} {1/np.sqrt(1+np.sin(th)):>15.6f}")
print(f"\nmax |SDP - formula| over the scan = {np.max(np.abs(pc_num - pc_formula)):.2e}")

# ---------------------------------------------------------------
# 11. Figure (all graphs in one figure)
# ---------------------------------------------------------------
fig = plt.figure(figsize=(16, 13))
colors = {2: "tab:red", 3: "tab:blue", 4: "tab:green"}
styles = {2: "-", 3: "-", 4: "--"}

# (a) SW vs p
ax1 = fig.add_subplot(2, 2, 1)
for m in curve_ms:
    ax1.plot(p_grid, sw_curves[m], marker="o", ms=4, ls=styles[m], color=colors[m], label=f"m = {m}")
    ax1.axvline(pc[m], color=colors[m], ls=":", lw=1)
ax1.set_xlabel("Werner parameter p")
ax1.set_ylabel("Steerable weight SW")
ax1.set_title("(a) Steerable weight vs. noise")
ax1.legend()
ax1.grid(alpha=0.3)

# (b) critical p vs m
ax2 = fig.add_subplot(2, 2, 2)
ms = list(pc.keys())
bars = ax2.bar([str(m) for m in ms], [pc[m] for m in ms], color="tab:purple", alpha=0.8)
for b, m in zip(bars, ms):
    ax2.text(b.get_x() + b.get_width() / 2, pc[m] + 0.01, f"{pc[m]:.4f}", ha="center")
ax2.axhline(1 / np.sqrt(2), color="gray", ls="--", lw=1, label="1/sqrt(2)")
ax2.axhline(0.5, color="black", ls="--", lw=1, label="1/2 (continuum limit)")
ax2.set_ylim(0.4, 0.8)
ax2.set_xlabel("Number of measurement settings m")
ax2.set_ylabel("Critical visibility p_c")
ax2.set_title("(b) More settings reveal steering at lower p")
ax2.legend()
ax2.grid(alpha=0.3, axis="y")

# (c) 3D surface
ax3 = fig.add_subplot(2, 2, 3, projection="3d")
P, T = np.meshgrid(ps, np.degrees(thetas))
surf = ax3.plot_surface(P, T, SW_surf, cmap="viridis", edgecolor="none", alpha=0.92)
ax3.plot(pc_formula, np.degrees(th_line), np.zeros_like(th_line), color="red", lw=2.5, label="critical curve")
ax3.set_xlabel("p")
ax3.set_ylabel("angle theta [deg]")
ax3.set_zlabel("SW")
ax3.set_title("(c) SW(p, theta) for two settings")
ax3.view_init(elev=28, azim=-125)
ax3.legend(loc="upper left")
fig.colorbar(surf, ax=ax3, shrink=0.6, pad=0.1)

# (d) 3D Bloch sphere
ax4 = fig.add_subplot(2, 2, 4, projection="3d")
u, v = np.mgrid[0:2 * np.pi:40j, 0:np.pi:20j]
ax4.plot_wireframe(np.cos(u) * np.sin(v), np.sin(u) * np.sin(v), np.cos(v), color="lightgray", lw=0.4)
for p_val, col in [(1.0, "tab:red"), (0.8, "tab:orange"), (pc[3], "tab:blue")]:
    sig = assemblage(werner_state(p_val), SETTINGS[3])
    pts = []
    for x in range(3):
        for a in (0, 1):
            s = sig[x][a]
            pts.append([np.real(np.trace(s @ P_)) / np.real(np.trace(s)) for P_ in (SX, SY, SZ)])
    pts = np.array(pts)
    ax4.scatter(pts[:, 0], pts[:, 1], pts[:, 2], s=70, color=col, label=f"p = {p_val:.3f}")
ax4.set_xlabel("x")
ax4.set_ylabel("y")
ax4.set_zlabel("z")
ax4.set_title("(d) Bob's conditional states (m = 3)")
ax4.legend(loc="upper left")
ax4.set_box_aspect((1, 1, 1))

plt.tight_layout()
plt.show()

4. Walking through the code

Building blocks (sections 1 and 2)

werner_state(p) builds $\rho(p)=p|\psi^-\rangle\langle\psi^-|+(1-p)\mathbb 1/4$ as a $4\times4$ matrix. projector(n, a) returns $P_{a|n}=\tfrac12\big(\mathbb 1+(-1)^a,\vec n\cdot\vec\sigma\big)$, the spin projector along an arbitrary axis, normalized internally so that non-unit vectors like $(1,1,1)$ also work.

assemblage computes $\sigma_{a|x}=\mathrm{Tr}A[(P{a|x}\otimes\mathbb 1)\rho]$. The product $(P\otimes\mathbb 1)\rho$ is reshaped into a rank-4 tensor with indices $(i,j,i’,k)$, where the first pair labels Alice and Bob on the row side and the second pair labels them on the column side. The partial trace over Alice is then the contraction "ijik->jk", which sets the two Alice indices equal and sums over them.

The solver helper prefers Clarabel, an interior-point solver that is accurate to roughly $10^{-9}$ and well suited to small SDPs. If it is missing, it falls back to SCS with a tight tolerance.

SDP (A): steerable weight (section 3)

The function enumerates all deterministic strategies with itertools.product((0, 1), repeat=m), which gives $2^m$ tuples $\lambda=(\lambda_1,\dots,\lambda_m)$. Each strategy gets a Hermitian $2\times2$ variable $\sigma_\lambda$ with the constraint s >> 0.

For every setting $x$ and outcome $a$, the expression part is $\sum_{\lambda:\lambda_x=a}\sigma_\lambda$. The constraint sigma[x][a] - part >> 0 enforces $\sigma_{a|x}-\sum_\lambda D(a|x,\lambda)\sigma_\lambda\succeq0$. The objective maximizes $\mathrm{Re},\mathrm{Tr}\sum_\lambda\sigma_\lambda=\mu$, and the function returns $\mathrm{SW}=1-\mu^\star$, clipped at zero to remove tiny negative solver noise.

SDP (B) and SDP (C): finding $p_c$ (sections 4 and 5)

critical_p_bisection is the straightforward approach. It repeatedly solves the steerable-weight SDP at a trial value of $p$, moving the upper bound down if $\mathrm{SW}>10^{-6}$ and the lower bound up otherwise. With 22 iterations this costs 22 SDPs per critical visibility.

critical_p_direct is the accelerated version. The only change is that $p$ becomes a scalar cp.Variable(), and the inequality becomes the equality

$$
\sum_{\lambda:\lambda_x=a}\sigma_\lambda=p,\sigma^{\rm singlet}_{a|x}+(1-p)\frac{\mathbb 1}{4}.
$$

Maximizing $p$ then returns $p_c$ directly. The feasible set in $p$ is an interval $[0,p_c]$, because the problem is convex, so the maximum is exactly the threshold. One SDP replaces 22, and the program stays small: the matrices are always $2\times2$, with $2^m$ of them.

Settings (section 6)

The four sets of measurement directions are:

$m=2$: the $x$ and $z$ axes.
$m=3$: the $x$, $y$ and $z$ axes. The $y$ axis makes some matrices complex, which is why the variables are declared hermitian=True.
$m=4$: four cube diagonals.
$m=6$: six axes through the vertices of an icosahedron, using the golden ratio $\varphi=(1+\sqrt5)/2$.

Experiments (sections 7 to 10)

Experiment 1 loops over the settings with the fast solver and prints a table. Section 8 compares bisection against the one-shot SDP on $m=4$, and both must agree on $p_c$. Experiment 2 sweeps $p$ from $0.4$ to $1$ and evaluates $\mathrm{SW}$ for $m=2,3,4$.

Experiment 3 fixes the first setting along $z$ and rotates the second by an angle $\theta$ in the $x$–$z$ plane. It builds an $18\times21$ grid of $\mathrm{SW}(p,\theta)$ values for the 3D surface, then computes $p_c(\theta)$ by the fast SDP at 40 angles and compares it with $1/\sqrt{1+\sin\theta}$.

Plotting (section 11)

All four graphs are placed in a single figure with a $2\times2$ layout. The m=4 curve is dashed so that it remains visible even where it lies exactly on top of the m=3 curve. In panel (d), Bob’s conditional Bloch vector is recovered from each conditional state $s$ as $r_i=\mathrm{Tr}(s,\sigma_i)/\mathrm{Tr}(s)$.

5. Console output

================================================================
Experiment 1: critical visibility p_c of the Werner state
================================================================
settings m  strategies 2^m    p_c (SDP)   time [s]
         2               4     0.707107      0.138
         3               8     0.577350      0.258
         4              16     0.577350      0.436
         6              64     0.539345      2.713

Reference: 1/sqrt(2) = 0.707107,  1/sqrt(3) = 0.577350,  1/2 = 0.500000

================================================================
Speed comparison (m = 4 settings)
================================================================
bisection (22 SDPs): p_c = 0.577351, time = 8.040 s
one-shot SDP       : p_c = 0.577350, time = 0.151 s
speed-up           : x53.2

================================================================
Experiment 2: steerable weight SW(p)
================================================================
     p     SW (m=2)     SW (m=3)     SW (m=4)
  0.40     0.000000     0.000000     0.000000
  0.46     0.000000     0.000000     0.000000
  0.52     0.000000     0.000000     0.000000
  0.58     0.000000     0.006269     0.006269
  0.64     0.000000     0.148231     0.148231
  0.70     0.000000     0.290192     0.290192
  0.76     0.180589     0.432154     0.432154
  0.82     0.385442     0.574115     0.574115
  0.88     0.590294     0.716077     0.716077
  0.94     0.795147     0.858038     0.858038
  1.00     1.000000     1.000000     1.000000

================================================================
Experiment 3: two settings at angle theta
================================================================
 theta [deg]    p_c (SDP)   1/sqrt(1+sin)
        15.0     0.891289        0.891289
        30.0     0.816497        0.816497
        45.0     0.765367        0.765367
        60.0     0.732051        0.732051
        75.0     0.713208        0.713208
        90.0     0.707107        0.707107

max |SDP - formula| over the scan = 1.85e-08

Reading the console output

Experiment 1. The SDP reproduces $p_c=1/\sqrt2\approx0.7071$ for two settings and $p_c=1/\sqrt3\approx0.5774$ for three. Four cube diagonals give the same value as three orthogonal axes, $0.5774$, so the extra setting adds nothing here. Six icosahedron axes lower the threshold to about $0.5393$, approaching the continuum limit of $1/2$. The runtime column shows the cost growing with the number of strategies $2^m$, yet even $m=6$ remains cheap.
Speed comparison. Both methods return the same $p_c$ to within the bisection tolerance, but the one-shot SDP is faster by more than an order of magnitude because it replaces 22 SDP solves with a single one.
Experiment 2. The table shows $\mathrm{SW}=0$ below the threshold and a linear rise above it up to $\mathrm{SW}=1$ at $p=1$. For two settings, for example, $\mathrm{SW}\approx0.3172$ at $p=0.8$.
Experiment 3. The SDP values of $p_c(\theta)$ match $1/\sqrt{1+\sin\theta}$ essentially to solver precision, so the closed form is confirmed numerically.

6. Graphs

Panel (a): steerable weight versus noise

Every curve is flat at zero until its own critical visibility, marked by the dotted vertical line, and then rises linearly to $\mathrm{SW}=1$ at $p=1$. In all the cases computed here the rise follows

$$
\mathrm{SW}(p)=\max!\left(0,\ \frac{p-p_c}{1-p_c}\right),
$$

so the whole curve is determined by $p_c$. The red curve ($m=2$) starts rising at $p\approx0.707$, while the blue and green curves ($m=3$ and $m=4$) start at $p\approx0.577$. The blue and green curves coincide exactly, which is why the dashed green line sits on top of the blue one: four cube diagonals reveal no more steering than three orthogonal axes for the Werner state.

Panel (b): critical visibility versus number of settings

This is the central result. Each additional well-chosen setting lets Alice detect steering in noisier states, so the threshold falls from $0.7071$ at $m=2$ to $0.5774$ at $m=3$ and $0.5393$ at $m=6$. The dashed line at $1/2$ is the continuum limit, which no finite set of settings can cross. Panel (b) also shows that the number of settings is not the only thing that matters: their geometry does too, since $m=4$ with cube diagonals gains nothing over $m=3$.

Panel (c): the 3D surface $\mathrm{SW}(p,\theta)$

The two horizontal axes are the Werner parameter $p$ and the angle $\theta$ between Alice’s two measurement axes, and the height is the steerable weight. The surface is a flat plain at $\mathrm{SW}=0$ in the unsteerable region, and it rises into a ramp that reaches $\mathrm{SW}=1$ at $p=1$ for large $\theta$. The red curve on the floor is $p_c(\theta)=1/\sqrt{1+\sin\theta}$, and it follows the foot of the ramp exactly.

At small $\theta$ the two measurements are nearly the same, so they carry almost no incompatibility and the critical visibility approaches $1$. As $\theta$ grows to $90^\circ$ the settings become maximally incompatible and the threshold falls to $1/\sqrt2$. The steerable region therefore widens as the measurements become more incompatible.

Panel (d): Bob’s conditional states on the Bloch sphere

For $m=3$, each setting and outcome produces one point on or inside the Bloch sphere. At $p=1$ (red) all six points lie on the surface at $\pm\hat x,\pm\hat y,\pm\hat z$, which means Alice steers Bob into pure states. Lowering $p$ shrinks all six points toward the center, giving orange at $p=0.8$ and blue at $p=1/\sqrt3\approx0.577$. At the critical visibility the shrunken octahedron has just become small enough that a local hidden state ensemble can reproduce it, which is where the SDP reports $\mathrm{SW}=0$.

7. Why the fast version matters

Each steering test involves $2^m$ matrix variables, so the cost grows exponentially with the number of settings. Locating a threshold by bisection multiplies that cost by the number of iterations. Turning $p$ into an SDP variable removes this factor entirely, because the threshold becomes the optimum of a single convex program. The comparison section of the console output makes the difference concrete, and the saving is largest in the regime where each SDP is already expensive, with $m=6$ and beyond.

8. Conclusion

Steerability with a finite set of measurements is a semidefinite feasibility problem. Solving it for the Werner state reproduces the textbook thresholds $1/\sqrt2$, $1/\sqrt3$, and the closed form $1/\sqrt{1+\sin\theta}$ for two settings at angle $\theta$. The same code can be applied to any two-qubit state or any set of measurements by changing werner_state and SETTINGS.

Computing the Maximal Quantum Violation of a Bell Inequality with the NPA Hierarchy

October 12, 2026

The Bell–CHSH inequality is a fence that every local hidden variable theory must respect. Quantum mechanics climbs over it, but only up to a precise height called the Tsirelson bound. In this post we compute that height numerically with a semidefinite program (SDP), using the Navascués–Pironio–Acín (NPA) hierarchy, and then explore how the bound changes when we reweight the inequality.

1. The Problem

Alice and Bob each choose between two $\pm1$-valued observables, $A_0, A_1$ and $B_0, B_1$. With correlators $E_{ij}=\langle A_i B_j\rangle$, the CHSH expression is

$$
S = E_{00} + E_{01} + E_{10} - E_{11}.
$$

Local hidden variable models obey $|S|\le 2$.
Quantum mechanics obeys $|S|\le 2\sqrt{2}$ (Tsirelson bound).
No-signalling theories can reach $|S|\le 4$.

We also study a weighted CHSH family with two parameters $\alpha,\beta\ge 0$:

$$
S(\alpha,\beta)=\alpha,(E_{00}+E_{01})+\beta,(E_{10}-E_{11}).
$$

For this family the classical bound is $2\max(\alpha,\beta)$, and Tsirelson’s vector construction gives the analytic quantum bound

$$
S_Q(\alpha,\beta)=2\sqrt{\alpha^{2}+\beta^{2}}.
$$

This gives us an exact formula to check the SDP against.

2. The NPA Hierarchy in a Nutshell

Choose a finite list of operator words ${w_i}$ built from $A_0,A_1,B_0,B_1$. These operators satisfy

$$
A_i^2=B_j^2=\mathbb{1},\qquad [A_i,B_j]=0 .
$$

Define the moment matrix

$$
\Gamma_{ij}=\langle \psi|,w_i^{\dagger}w_j,|\psi\rangle .
$$

Every quantum realisation produces a $\Gamma$ that is positive semidefinite, has unit diagonal, and whose entries are equal whenever the underlying operator products reduce to the same word. The NPA relaxation drops everything else:

$$
\max_{\Gamma}\ \sum_{ij} c_{ij},\langle A_iB_j\rangle
\quad\text{s.t.}\quad \Gamma\succeq 0,\ \ \Gamma_{\mathbb 1\mathbb 1}=1,\ \ \Gamma_{ij}=\Gamma_{kl}\ \text{if}\ w_i^\dagger w_j \equiv w_k^\dagger w_l .
$$

Each level of the hierarchy uses longer words and gives a tighter outer approximation of the quantum set:

For CHSH, already the first level is tight, so every level returns $2\sqrt2$.

3. The Complete Code

import time
import itertools
import numpy as np
import scipy.sparse as sp
import cvxpy as cp
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D

# ------------------------------------------------------------------
# Solver selection
# ------------------------------------------------------------------
SOLVER = "CLARABEL" if "CLARABEL" in cp.installed_solvers() else "SCS"

# ------------------------------------------------------------------
# Operator-word algebra: A_i^2 = B_j^2 = 1, Alice and Bob commute
# A word is a pair (alice_sequence, bob_sequence) of tuples of 0/1.
# ------------------------------------------------------------------
def reduce_word(seq):
    out = []
    for s in seq:
        if out and out[-1] == s:
            out.pop()          # X_i X_i = identity
        else:
            out.append(s)
    return tuple(out)

def party_words(max_len):
    words = [()]
    frontier = [()]
    for _ in range(max_len):
        new = []
        for w in frontier:
            for s in (0, 1):
                if not w or w[-1] != s:
                    new.append(w + (s,))
        words += new
        frontier = new
    return words

def build_words(level):
    if level == "1+AB":
        base, extra = 1, True
    else:
        base, extra = int(level), False
    A = party_words(base)
    B = party_words(base)
    words = [(a, b) for a in A for b in B if len(a) + len(b) <= base]
    if extra:
        words += [((i,), (j,)) for i in (0, 1) for j in (0, 1)]
    return words

def canon(a, b):
    k1 = (a, b)
    k2 = (a[::-1], b[::-1])    # adjoint (real relaxation)
    return min(k1, k2)

def moment_key(u, v):
    a = reduce_word(u[0][::-1] + v[0])
    b = reduce_word(u[1][::-1] + v[1])
    return canon(a, b)

def word_name(w):
    a, b = w
    s = "".join(f"A{i}" for i in a) + "".join(f"B{j}" for j in b)
    return s if s else "I"

# ------------------------------------------------------------------
# Build the NPA SDP once (coefficients are a cvxpy Parameter)
# ------------------------------------------------------------------
def make_npa(level):
    words = build_words(level)
    n = len(words)
    keys, entries = {}, []
    for i, u in enumerate(words):
        for j, v in enumerate(words):
            k = moment_key(u, v)
            if k not in keys:
                keys[k] = len(keys)
            entries.append((i * n + j, keys[k]))
    K = len(keys)
    rows = [e[0] for e in entries]
    cols = [e[1] for e in entries]
    S = sp.coo_matrix((np.ones(len(entries)), (rows, cols)), shape=(n * n, K)).tocsr()

    y = cp.Variable(K)
    Gamma = cp.reshape(cp.Constant(S) @ y, (n, n), order="F")
    constraints = [Gamma >> 0, y[keys[((), ())]] == 1]

    C = cp.Parameter((2, 2))
    objective = 0
    for i in range(2):
        for j in range(2):
            objective += C[i, j] * y[keys[((i,), (j,))]]
    prob = cp.Problem(cp.Maximize(objective), constraints)
    return dict(prob=prob, C=C, y=y, Gamma=Gamma, keys=keys, words=words, n=n, K=K)

def solve_npa(model, coeffs):
    model["C"].value = np.array(coeffs, dtype=float)
    model["prob"].solve(solver=SOLVER)
    return model["prob"].value

# ------------------------------------------------------------------
# Example 1: CHSH at several levels of the hierarchy
# ------------------------------------------------------------------
CHSH = [[1, 1], [1, -1]]
levels = [1, "1+AB", 2, 3]
level_values, level_info = [], []

print(f"Solver: {SOLVER}\n")
print(f"{'level':>6} | {'matrix size':>11} | {'#moments':>8} | {'CHSH max':>10} | {'time [s]':>8}")
print("-" * 58)
for lv in levels:
    t0 = time.time()
    model = make_npa(lv)
    val = solve_npa(model, CHSH)
    dt = time.time() - t0
    level_values.append(val)
    level_info.append(model)
    print(f"{str(lv):>6} | {model['n']:>5d}x{model['n']:<5d} | {model['K']:>8d} | {val:>10.6f} | {dt:>8.3f}")

print(f"\nTsirelson bound 2*sqrt(2) = {2 * np.sqrt(2):.6f}")

# Optimal correlators at level 1
m1 = level_info[0]
yv = m1["y"].value
print("\nOptimal correlators (level 1):")
for i in range(2):
    for j in range(2):
        print(f"  E{i}{j} = {yv[m1['keys'][((i,), (j,))]]: .6f}")
print(f"  1/sqrt(2) = {1 / np.sqrt(2):.6f}")

# ------------------------------------------------------------------
# Classical bound by brute force and an explicit qubit strategy
# ------------------------------------------------------------------
classical = max(a0 * b0 + a0 * b1 + a1 * b0 - a1 * b1
                for a0, a1, b0, b1 in itertools.product([-1, 1], repeat=4))

X = np.array([[0.0, 1.0], [1.0, 0.0]])
Z = np.diag([1.0, -1.0])
phi = np.array([1.0, 0.0, 0.0, 1.0]) / np.sqrt(2)

def obs(t):
    return np.cos(t) * Z + np.sin(t) * X

def corr(a, b):
    return float(phi @ np.kron(obs(a), obs(b)) @ phi)

a0, a1, b0, b1 = 0.0, np.pi / 2, np.pi / 4, -np.pi / 4
qubit_S = corr(a0, b0) + corr(a0, b1) + corr(a1, b0) - corr(a1, b1)
print(f"\nClassical bound (brute force)   : {classical:.6f}")
print(f"Explicit two-qubit strategy     : {qubit_S:.6f}")

# ------------------------------------------------------------------
# Example 2: weighted CHSH on a grid (fast: one compiled problem, many solves)
# ------------------------------------------------------------------
grid = np.linspace(0.0, 1.0, 11)
AA, BB = np.meshgrid(grid, grid, indexing="ij")
Z_sdp = np.zeros_like(AA)

fast = make_npa(1)
t0 = time.time()
for i, a in enumerate(grid):
    for j, b in enumerate(grid):
        Z_sdp[i, j] = solve_npa(fast, [[a, a], [b, -b]])
t_grid = time.time() - t0

Z_cl = 2.0 * np.maximum(AA, BB)
Z_an = 2.0 * np.sqrt(AA ** 2 + BB ** 2)
max_err = np.max(np.abs(Z_sdp - Z_an))
print(f"\nWeighted CHSH grid: {grid.size * grid.size} SDPs solved in {t_grid:.2f} s")
print(f"Max |SDP - 2*sqrt(a^2+b^2)| over the grid = {max_err:.2e}")

# ------------------------------------------------------------------
# Geometry of the optimal level-1 moment matrix (Tsirelson vectors)
# ------------------------------------------------------------------
G = np.array(m1["Gamma"].value)
G = 0.5 * (G + G.T)
names = [word_name(w) for w in m1["words"]]
w_eig, U = np.linalg.eigh(G)
order = np.argsort(w_eig)[::-1]
V = U[:, order[:3]] * np.sqrt(np.clip(w_eig[order[:3]], 0.0, None))

# ------------------------------------------------------------------
# One combined figure
# ------------------------------------------------------------------
fig = plt.figure(figsize=(21, 12))
fig.suptitle("NPA hierarchy for the CHSH inequality", fontsize=18, fontweight="bold")

# (1) 3D surface: SDP value vs classical bound
ax1 = fig.add_subplot(2, 3, 1, projection="3d")
ax1.plot_surface(AA, BB, Z_sdp, cmap="viridis", alpha=0.92, linewidth=0, antialiased=True)
ax1.plot_wireframe(AA, BB, Z_cl, color="crimson", linewidth=0.8, alpha=0.85)
ax1.set_xlabel(r"$\alpha$")
ax1.set_ylabel(r"$\beta$")
ax1.set_zlabel("Bell value")
ax1.set_title("Quantum maximum (SDP surface) vs classical bound (red mesh)")
ax1.view_init(elev=24, azim=-58)
ax1.legend(handles=[Line2D([0], [0], color="seagreen", lw=6, label="NPA level 1"),
                    Line2D([0], [0], color="crimson", lw=2, label="Classical")],
           loc="upper left")

# (2) 3D Tsirelson vectors
ax2 = fig.add_subplot(2, 3, 2, projection="3d")
uu, vv = np.mgrid[0:2 * np.pi:40j, 0:np.pi:20j]
ax2.plot_wireframe(np.cos(uu) * np.sin(vv), np.sin(uu) * np.sin(vv), np.cos(vv),
                   color="lightgray", alpha=0.35, linewidth=0.4)
colors = {"A0": "tab:blue", "A1": "tab:cyan", "B0": "tab:red", "B1": "tab:orange"}
for nm, col in colors.items():
    vec = V[names.index(nm)]
    ax2.quiver(0, 0, 0, vec[0], vec[1], vec[2], color=col, linewidth=2.5,
               arrow_length_ratio=0.1)
    ax2.text(vec[0] * 1.15, vec[1] * 1.15, vec[2] * 1.15, nm, color=col, fontsize=12)
ax2.set_xlim(-1.2, 1.2)
ax2.set_ylim(-1.2, 1.2)
ax2.set_zlim(-1.2, 1.2)
ax2.set_box_aspect((1, 1, 1))
ax2.set_title("Optimal unit vectors from the moment matrix")
ax2.view_init(elev=22, azim=35)

# (3) Moment matrix heatmap
ax3 = fig.add_subplot(2, 3, 3)
im = ax3.imshow(G, cmap="RdBu_r", vmin=-1, vmax=1)
ax3.set_xticks(range(len(names)))
ax3.set_yticks(range(len(names)))
ax3.set_xticklabels(names)
ax3.set_yticklabels(names)
for i in range(G.shape[0]):
    for j in range(G.shape[1]):
        ax3.text(j, i, f"{G[i, j]:.2f}", ha="center", va="center", fontsize=10)
ax3.set_title(r"Optimal moment matrix $\Gamma$ (level 1)")
fig.colorbar(im, ax=ax3, fraction=0.046, pad=0.04)

# (4) Slice beta = 1
ax4 = fig.add_subplot(2, 3, 4)
fine = np.linspace(0, 1, 200)
ax4.plot(fine, 2 * np.sqrt(fine ** 2 + 1), color="black", lw=2, label=r"$2\sqrt{\alpha^2+1}$ (analytic)")
ax4.plot(fine, 2 * np.maximum(fine, 1.0), color="crimson", lw=2, ls="--", label="Classical bound")
ax4.scatter(grid, Z_sdp[:, -1], color="tab:green", zorder=5, s=45, label="SDP (NPA level 1)")
ax4.set_xlabel(r"$\alpha$ (with $\beta=1$)")
ax4.set_ylabel("Bell value")
ax4.set_title(r"Slice of the surface at $\beta=1$")
ax4.grid(alpha=0.3)
ax4.legend()

# (5) Bar chart of bounds
ax5 = fig.add_subplot(2, 3, 5)
bar_labels = ["Classical", "NPA 1", "NPA 1+AB", "NPA 2", "NPA 3", "Qubit\nstrategy", "No-signalling\n(PR box)"]
bar_vals = [classical] + list(level_values) + [qubit_S, 4.0]
bar_cols = ["crimson"] + ["tab:green"] * 4 + ["tab:blue", "gray"]
bars = ax5.bar(bar_labels, bar_vals, color=bar_cols)
ax5.axhline(2 * np.sqrt(2), color="black", ls="--", lw=1.2, label=r"$2\sqrt{2}$")
for b_, v_ in zip(bars, bar_vals):
    ax5.text(b_.get_x() + b_.get_width() / 2, v_ + 0.04, f"{v_:.4f}", ha="center", fontsize=9)
ax5.set_ylim(0, 4.5)
ax5.set_ylabel("CHSH value")
ax5.set_title("Bounds on the CHSH value")
ax5.legend()
plt.setp(ax5.get_xticklabels(), fontsize=9)

# (6) Numerical error map
ax6 = fig.add_subplot(2, 3, 6)
err = np.abs(Z_sdp - Z_an)
im6 = ax6.imshow(np.log10(err.T + 1e-16), origin="lower", extent=[0, 1, 0, 1],
                 aspect="auto", cmap="magma")
ax6.set_xlabel(r"$\alpha$")
ax6.set_ylabel(r"$\beta$")
ax6.set_title(r"$\log_{10}|\,\mathrm{SDP}-2\sqrt{\alpha^2+\beta^2}\,|$")
fig.colorbar(im6, ax=ax6, fraction=0.046, pad=0.04)

fig.tight_layout(rect=[0, 0, 1, 0.96])
plt.savefig("npa_chsh_summary.png", dpi=150)
plt.show()

4. Code Walkthrough

4.1 Operator algebra

reduce_word implements $X_i^2=\mathbb 1$ with a stack: whenever the same letter appears twice in a row, both copies cancel. Alice’s and Bob’s letters commute, so a word is stored as a pair (alice_sequence, bob_sequence) and each half is reduced independently.

party_words(L) lists all reduced words for one party up to length $L$. With two observables these are just alternating strings: $\mathbb 1, A_0, A_1, A_0A_1, A_1A_0, A_0A_1A_0,\dots$

build_words(level) then forms every product of an Alice word and a Bob word whose total length is at most the level. The special level "1+AB" adds the four products $A_iB_j$ to the level-1 list, giving a $9\times 9$ moment matrix.

4.2 Moment keys

The entry $\Gamma_{uv}=\langle u^\dagger v\rangle$ depends only on the reduced form of $u^\dagger v$. The function moment_key computes that reduced word, and canon identifies a word with its reversed (adjoint) version, since we work with the real relaxation. Two matrix entries share one SDP variable exactly when their keys coincide. This is how the linear constraints $\Gamma_{ij}=\Gamma_{kl}$ are enforced without writing a single equality by hand.

4.3 Building the SDP

make_npa assigns one scalar variable (an entry of y) to each distinct key. A sparse $0/1$ matrix S maps y to the flattened $\Gamma$, and cp.reshape turns the result into a square matrix. Because $\Gamma$ is symmetric by construction, the reshape order is irrelevant. The constraints are:

Gamma >> 0, which is the positive-semidefiniteness condition,
y[identity] == 1, which is the normalisation $\langle\mathbb 1\rangle=1$.

The unit diagonal comes for free, because $u^\dagger u$ always reduces to the identity.

The Bell coefficients live in a cp.Parameter, so the objective $\sum_{ij}c_{ij}\langle A_iB_j\rangle$ is affine in the variables for any value of the parameter.

4.4 Speed-up for the parameter scan

A naive scan would rebuild and re-canonicalise a new problem for each of the 121 grid points. Here the problem is built once, and only C.value changes between solves. Because the parameter multiplies a variable in a disciplined-parametrized (DPP) way, cvxpy reuses the compiled problem and only re-runs the solver. The sparse map S also keeps the matrix construction cheap at higher levels.

4.5 Validation

Three independent checks surround the SDP:

A brute-force search over the 16 deterministic strategies gives the classical bound $2$.
An explicit two-qubit strategy with the state $|\Phi^+\rangle$ and measurement angles $0,\pi/2$ (Alice) and $\pm\pi/4$ (Bob) gives $4\cdot\tfrac{1}{\sqrt2}=2\sqrt2$, so the upper bound from the SDP is actually attained.
The analytic formula $2\sqrt{\alpha^2+\beta^2}$ is compared with the SDP over the whole grid.

4.6 Geometry

At level 1 the matrix $\Gamma$ is a Gram matrix of unit vectors, $\Gamma_{ij}=\vec v_i\cdot\vec v_j$. An eigendecomposition $\Gamma=V V^{T}$ recovers those vectors, and we plot the three dominant components in 3D. This is Tsirelson’s construction made visible.

5. Console Output

Solver: CLARABEL

 level | matrix size | #moments |   CHSH max | time [s]
----------------------------------------------------------
     1 |     5x5     |       11 |   2.828427 |    0.043
  1+AB |     9x9     |       17 |   2.828427 |    0.025
     2 |    13x13    |       31 |   2.828427 |    0.022
     3 |    25x25    |       61 |   2.828427 |    0.074

Tsirelson bound 2*sqrt(2) = 2.828427

Optimal correlators (level 1):
  E00 =  0.707107
  E01 =  0.707107
  E10 =  0.707107
  E11 = -0.707107
  1/sqrt(2) = 0.707107

Classical bound (brute force)   : 2.000000
Explicit two-qubit strategy     : 2.828427

Weighted CHSH grid: 121 SDPs solved in 0.62 s
Max |SDP - 2*sqrt(a^2+b^2)| over the grid = 2.64e-08

Reading the console output

The table shows the CHSH maximum at every level of the hierarchy. All four entries should equal $2\sqrt2\approx 2.828427$ to solver precision, which confirms that the first level already captures the full quantum set for this scenario.
The optimal correlators come out as $|E_{ij}|=1/\sqrt2$, with $E_{11}$ carrying the opposite sign, so that $S = 4\cdot\tfrac{1}{\sqrt2}=2\sqrt2$.
The explicit qubit strategy reproduces the same number. Since the NPA value is an upper bound and the qubit strategy is a lower bound, they sandwich the true quantum maximum.
The final line reports the largest deviation between the SDP and the closed form $2\sqrt{\alpha^2+\beta^2}$ over all 121 grid points.

6. Result Image

7. Reading the Figure

Top left: 3D surface. The coloured surface is the quantum maximum computed by the SDP over $(\alpha,\beta)\in[0,1]^2$. The red mesh is the classical bound $2\max(\alpha,\beta)$. The two coincide along the axes, where the Bell expression degenerates to a trivially satisfiable one, and they separate most strongly along the diagonal $\alpha=\beta$, where the surface peaks at $2\sqrt2$ against a classical $2$. The gap between the surface and the mesh is the room available for quantum non-locality.

Top centre: Tsirelson vectors. The four arrows are the vectors $\vec a_0,\vec a_1,\vec b_0,\vec b_1$ extracted from $\Gamma$. $\vec a_0$ and $\vec a_1$ are orthogonal, and likewise $\vec b_0$ and $\vec b_1$. Each $\vec b_j$ bisects the two Alice vectors, at $45^\circ$ to each. That is exactly the geometry that maximises

$$
\sum_{j}\bigl|c_{0j}\vec a_0+c_{1j}\vec a_1\bigr| ,
$$

and it mirrors the measurement angles of the optimal qubit strategy.

Top right: moment matrix. The block of correlators between Alice’s and Bob’s observables shows entries of magnitude $0.71\approx 1/\sqrt2$. The entries between $A_0$ and $A_1$, and between $B_0$ and $B_1$, are zero, which is the orthogonality seen in the 3D picture. The unit diagonal is the normalisation.

Bottom left: slice at $\beta=1$. Green dots (SDP) lie on the black analytic curve $2\sqrt{\alpha^2+1}$ and stay above the dashed classical line for every $\alpha<1$. At $\alpha=1$ the curve reaches $2\sqrt2$.

Bottom centre: bounds on one axis. The bars compare the classical value ($2$), the four NPA levels, the explicit qubit strategy, and the no-signalling maximum ($4$, attained by the Popescu–Rohrlich box). The NPA bars and the qubit bar have the same height, namely the Tsirelson bound, which sits strictly between the classical and no-signalling values.

Bottom right: numerical error. The map of $\log_{10}|\mathrm{SDP}-2\sqrt{\alpha^2+\beta^2}|$ shows that the SDP matches the analytic bound to solver tolerance everywhere on the grid, so the SDP and the closed form agree across the whole parameter range.

8. Takeaways

The NPA hierarchy turns the question “what is the largest Bell violation quantum mechanics allows?” into a sequence of SDPs with moment matrices that are positive semidefinite and carry the algebraic identifications between words.
For CHSH, level 1 is already exact, and the optimum $2\sqrt2$ is reached by a maximally entangled two-qubit state.
Building the SDP once with a cp.Parameter makes parameter scans cheap, and the same code can be pointed at other two-input, two-output Bell expressions simply by changing the coefficient matrix C.

Computing Quantum Entanglement Measures

October 11, 2026

Minimizing Relative Entropy and Robustness

How entangled is a quantum state? A yes/no answer (separable or entangled) is not enough for most purposes. We want a number. Many entanglement measures are defined as the minimum “distance” from the given state to the set of separable states. In this post we compute two of them numerically for two-qubit states:

the relative entropy of entanglement (a nonconvex optimization over separable states), and
the generalized robustness of entanglement (a semidefinite program).

We check both against exact formulas and visualize them over a two-parameter family of states, including 3D surfaces.

1. The Measures

Let $\mathrm{SEP}$ denote the set of separable states. A state is separable if it can be written as

$$
\sigma = \sum_k p_k, |a_k\rangle\langle a_k| \otimes |b_k\rangle\langle b_k| ,\qquad p_k\ge 0,\ \sum_k p_k = 1 .
$$

Relative entropy of entanglement

$$
E_R(\rho) = \min_{\sigma\in \mathrm{SEP}} S(\rho,|,\sigma),\qquad
S(\rho|\sigma) = \mathrm{Tr}\left[\rho\left(\log_2\rho - \log_2\sigma\right)\right].
$$

Generalized robustness of entanglement

This is the minimum amount of “noise” $\tau$ that must be mixed in to destroy the entanglement. Setting $Y=(1+s)\sigma$ turns it into a semidefinite program (SDP):

$$
R_G(\rho)=\min_{Y}\ \mathrm{Tr},Y - 1\quad \text{s.t.}\quad Y\succeq \rho,\quad Y^{T_B}\succeq 0 .
$$

For two qubits, “separable” is equivalent to “positive partial transpose” (the Peres–Horodecki criterion). This is why the constraint $Y^{T_B}\succeq 0$ is exact here.

The inequality connecting them

The log-robustness upper-bounds the relative entropy of entanglement:

$$
E_R(\rho)\ \le\ \log_2\bigl(1+R_G(\rho)\bigr).
$$

2. The Example Problem

Problem 1 (Werner states).

$$
\rho_W(p) = p,|\Phi^+\rangle\langle\Phi^+| + (1-p)\frac{I}{4},\qquad |\Phi^+\rangle=\frac{|00\rangle+|11\rangle}{\sqrt2}.
$$

With the singlet fraction $F=\langle\Phi^+|\rho_W|\Phi^+\rangle=\frac{3p+1}{4}$, the exact results for $F>1/2$ (that is, $p>1/3$) are

$$
E_R = 1 - h(F),\qquad h(x)=-x\log_2 x-(1-x)\log_2(1-x),\qquad R_G = 2F-1 = \frac{3p-1}{2},
$$

and both vanish for $p\le 1/3$. We recover these numerically.

Problem 2 (two-parameter family).

$$
\rho(p,\theta) = p,|\psi_\theta\rangle\langle\psi_\theta| + (1-p)\frac{I}{4},\qquad |\psi_\theta\rangle=\cos\theta,|00\rangle+\sin\theta,|11\rangle .
$$

The partial transpose has smallest eigenvalue $\frac{1-p}{4}-\frac{p\sin 2\theta}{2}$, so the state is separable exactly when

$$
p\le \frac{1}{1+2\sin 2\theta}.
$$

For pure states ($p=1$), the exact values are $E_R=h(\cos^2\theta)$ (the entanglement entropy) and $R_G=\sin 2\theta$.

We compute $E_R$ and $R_G$ on a $(p,\theta)$ grid and draw 3D surfaces.

3. Numerical Strategy

Relative entropy. The set $\mathrm{SEP}$ is hard to describe directly, so we parametrize it from the inside. We write

$$
\sigma(\mathbf{x}) = \sum_{k=1}^{K} p_k, |a_k b_k\rangle\langle a_k b_k|,\qquad
p_k=\frac{e^{\ell_k}}{\sum_j e^{\ell_j}},\qquad
|a\rangle=\begin{pmatrix}\cos\frac{\vartheta}{2}\ e^{i\varphi}\sin\frac{\vartheta}{2}\end{pmatrix}.
$$

This $\sigma$ is separable for every parameter value, so minimizing $S(\rho|\sigma(\mathbf x))$ with L-BFGS-B gives an upper bound on $E_R$ that converges to the true value as the optimizer converges. We use $K=8$ product states and several random restarts.

Robustness. This is a convex SDP, so we hand it to CVXPY.

4. Source Code

import time
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
from mpl_toolkits.mplot3d import Axes3D  # noqa: F401

try:
    import cvxpy as cp
except ImportError:
    import subprocess, sys
    subprocess.check_call([sys.executable, "-m", "pip", "install", "-q", "cvxpy"])
    import cvxpy as cp

np.set_printoptions(precision=5, suppress=True)
rng = np.random.default_rng(7)

# ---------------------------------------------------------------
# Basic quantum-information tools
# ---------------------------------------------------------------
def ket(v):
    return np.array(v, dtype=complex).reshape(-1, 1)

phi_plus = (ket([1, 0, 0, 0]) + ket([0, 0, 0, 1])) / np.sqrt(2)
I4 = np.eye(4, dtype=complex)

def werner(p):
    """rho = p |Phi+><Phi+| + (1-p) I/4"""
    return p * (phi_plus @ phi_plus.conj().T) + (1 - p) * I4 / 4

def family(p, theta):
    """rho = p |psi><psi| + (1-p) I/4, |psi> = cos(theta)|00> + sin(theta)|11>"""
    psi = np.cos(theta) * ket([1, 0, 0, 0]) + np.sin(theta) * ket([0, 0, 0, 1])
    return p * (psi @ psi.conj().T) + (1 - p) * I4 / 4

def partial_transpose(rho):
    """Partial transpose on the second qubit."""
    r = rho.reshape(2, 2, 2, 2)
    return r.transpose(0, 3, 2, 1).reshape(4, 4)

def fidelity_phi_plus(rho):
    return float(np.real((phi_plus.conj().T @ rho @ phi_plus)[0, 0]))

def binary_entropy(x):
    x = np.clip(x, 1e-15, 1 - 1e-15)
    return -x * np.log2(x) - (1 - x) * np.log2(1 - x)

def negativity(rho):
    ev = np.linalg.eigvalsh(partial_transpose(rho))
    return float(np.sum(np.abs(ev[ev < 0])))

def neg_entropy(rho):
    """Tr rho log2 rho"""
    w = np.linalg.eigvalsh(rho)
    w = w[w > 1e-14]
    return float(np.sum(w * np.log2(w)))

# ---------------------------------------------------------------
# (1) Relative entropy of entanglement
#     sigma = sum_k p_k |a_k b_k><a_k b_k|  (separable by construction)
# ---------------------------------------------------------------
K = 8  # number of product states

def build_sigma(x):
    logits = x[:K]
    ang = x[K:].reshape(4, K)
    w = np.exp(logits - logits.max())
    w /= w.sum()
    ta, pa, tb, pb = ang
    A = np.stack([np.cos(ta / 2), np.exp(1j * pa) * np.sin(ta / 2)], axis=1)
    B = np.stack([np.cos(tb / 2), np.exp(1j * pb) * np.sin(tb / 2)], axis=1)
    psi = (A[:, :, None] * B[:, None, :]).reshape(K, 4)
    return np.einsum("k,ki,kj->ij", w, psi, psi.conj())

def make_objective(rho):
    c = neg_entropy(rho)
    def f(x):
        sigma = build_sigma(x)
        w, V = np.linalg.eigh(sigma)
        lw = np.log2(np.clip(w, 1e-13, None))
        log_sigma = (V * lw) @ V.conj().T
        return c - float(np.real(np.trace(rho @ log_sigma)))
    return f

def relative_entropy_of_entanglement(rho, n_starts=3, x0=None, maxiter=400):
    f = make_objective(rho)
    best_val, best_x = np.inf, None
    starts = []
    if x0 is not None:
        starts.append(x0)
    for _ in range(n_starts):
        starts.append(np.concatenate([rng.normal(size=K),
                                      rng.uniform(0, 2 * np.pi, size=4 * K)]))
    for s in starts:
        res = minimize(f, s, method="L-BFGS-B", options={"maxiter": maxiter})
        if res.fun < best_val:
            best_val, best_x = res.fun, res.x
    return max(best_val, 0.0), best_x

# ---------------------------------------------------------------
# (2) Generalized robustness of entanglement (SDP)
# ---------------------------------------------------------------
def generalized_robustness(rho):
    rho_r = np.real(rho)
    Y = cp.Variable((4, 4), symmetric=True)
    cons = [Y - rho_r >> 0, cp.partial_transpose(Y, dims=[2, 2], axis=1) >> 0]
    prob = cp.Problem(cp.Minimize(cp.trace(Y) - 1), cons)
    prob.solve(solver=cp.SCS, eps=1e-8)
    return max(float(prob.value), 0.0)

# ---------------------------------------------------------------
# (3) Werner states: numerical result vs. analytic formulas
# ---------------------------------------------------------------
t_start = time.time()
p_list = np.linspace(0, 1, 13)
er_w, rg_w, neg_w, er_exact, rg_exact = [], [], [], [], []
for p in p_list:
    rho = werner(p)
    F = fidelity_phi_plus(rho)
    er, _ = relative_entropy_of_entanglement(rho, n_starts=3)
    er_w.append(er)
    rg_w.append(generalized_robustness(rho))
    neg_w.append(negativity(rho))
    er_exact.append(1 - binary_entropy(F) if F > 0.5 else 0.0)
    rg_exact.append(max(2 * F - 1, 0.0))
er_w, rg_w, neg_w = map(np.array, (er_w, rg_w, neg_w))
er_exact, rg_exact = np.array(er_exact), np.array(rg_exact)

print("=== Werner states: rho = p|Phi+><Phi+| + (1-p) I/4 ===")
print(f"{'p':>6} {'F':>7} {'E_R(num)':>10} {'E_R(exact)':>11} {'R_G(num)':>10} {'R_G(exact)':>11} {'N':>7}")
for p, a, b, c, d, n in zip(p_list, er_w, er_exact, rg_w, rg_exact, neg_w):
    F = (3 * p + 1) / 4
    print(f"{p:6.3f} {F:7.4f} {a:10.6f} {b:11.6f} {c:10.6f} {d:11.6f} {n:7.4f}")
print(f"max |E_R(num) - E_R(exact)| = {np.max(np.abs(er_w - er_exact)):.2e}")
print(f"max |R_G(num) - R_G(exact)| = {np.max(np.abs(rg_w - rg_exact)):.2e}")

# ---------------------------------------------------------------
# (4) Two-parameter family: rho(p, theta)
# ---------------------------------------------------------------
p_grid = np.linspace(0, 1, 11)
th_grid = np.linspace(0, np.pi / 4, 9)
ER = np.zeros((len(th_grid), len(p_grid)))
RG = np.zeros_like(ER)
for i, th in enumerate(th_grid):
    x_prev = None
    for j, p in enumerate(p_grid):
        rho = family(p, th)
        ER[i, j], x_prev = relative_entropy_of_entanglement(rho, n_starts=2, x0=x_prev)
        RG[i, j] = generalized_robustness(rho)

# pure-state check (p = 1): E_R = entanglement entropy = h(cos^2 theta)
er_pure_exact = binary_entropy(np.cos(th_grid) ** 2)
print("\n=== Pure states (p = 1): E_R vs. entanglement entropy ===")
for th, a, b in zip(th_grid, ER[:, -1], er_pure_exact):
    print(f"theta = {th:6.4f}   E_R(num) = {a:.6f}   entropy = {b:.6f}")
print(f"\nmax over grid of [E_R - log2(1+R_G)] = {np.max(ER - np.log2(1 + RG)):.2e}  (<= 0 up to numerical tolerance)")
print(f"Total computation time: {time.time() - t_start:.1f} s")

# ---------------------------------------------------------------
# (5) One combined figure
# ---------------------------------------------------------------
P, TH = np.meshgrid(p_grid, th_grid)
fig = plt.figure(figsize=(21, 12))

ax1 = fig.add_subplot(2, 3, 1)
pp = np.linspace(0, 1, 400)
FF = (3 * pp + 1) / 4
ax1.plot(pp, np.where(FF > 0.5, 1 - binary_entropy(FF), 0), "b-", label=r"$E_R$ exact")
ax1.plot(p_list, er_w, "ro", label=r"$E_R$ numerical")
ax1.plot(pp, np.maximum(2 * FF - 1, 0), "g-", label=r"$R_G$ exact")
ax1.plot(p_list, rg_w, "ks", mfc="none", label=r"$R_G$ numerical (SDP)")
ax1.plot(p_list, neg_w, "m^--", label=r"Negativity $N$")
ax1.axvline(1 / 3, color="gray", ls=":")
ax1.set_xlabel("p"); ax1.set_ylabel("entanglement")
ax1.set_title("(a) Werner states"); ax1.legend(); ax1.grid(alpha=0.3)

ax2 = fig.add_subplot(2, 3, 2, projection="3d")
s2 = ax2.plot_surface(P, TH, ER, cmap="viridis", edgecolor="k", linewidth=0.3)
ax2.set_xlabel("p"); ax2.set_ylabel(r"$\theta$"); ax2.set_zlabel(r"$E_R$")
ax2.set_title(r"(b) Relative entropy of entanglement $E_R(p,\theta)$")
ax2.view_init(elev=28, azim=-125)
fig.colorbar(s2, ax=ax2, shrink=0.55, pad=0.1)

ax3 = fig.add_subplot(2, 3, 3, projection="3d")
s3 = ax3.plot_surface(P, TH, RG, cmap="plasma", edgecolor="k", linewidth=0.3)
ax3.set_xlabel("p"); ax3.set_ylabel(r"$\theta$"); ax3.set_zlabel(r"$R_G$")
ax3.set_title(r"(c) Generalized robustness $R_G(p,\theta)$")
ax3.view_init(elev=28, azim=-125)
fig.colorbar(s3, ax=ax3, shrink=0.55, pad=0.1)

ax4 = fig.add_subplot(2, 3, 4)
cs = ax4.contourf(P, TH, ER, levels=20, cmap="viridis")
tt4 = np.linspace(0, np.pi / 4, 200)
ax4.plot(1 / (1 + 2 * np.sin(2 * tt4)), tt4, "r--", lw=2.5, label=r"$p=1/(1+2\sin 2\theta)$")
ax4.legend(loc="upper right")
fig.colorbar(cs, ax=ax4)
ax4.set_xlabel("p"); ax4.set_ylabel(r"$\theta$")
ax4.set_title(r"(d) $E_R$ map (dashed: separable boundary)")

ax5 = fig.add_subplot(2, 3, 5)
bound = np.log2(1 + RG).ravel()
sc = ax5.scatter(bound, ER.ravel(), c=P.ravel(), cmap="coolwarm", s=40, edgecolor="k")
mx = max(bound.max(), ER.max()) * 1.05
ax5.plot([0, mx], [0, mx], "k--", label=r"$E_R=\log_2(1+R_G)$")
fig.colorbar(sc, ax=ax5, label="p")
ax5.set_xlabel(r"$\log_2(1+R_G)$"); ax5.set_ylabel(r"$E_R$")
ax5.set_title(r"(e) Inequality $E_R \leq \log_2(1+R_G)$"); ax5.legend(); ax5.grid(alpha=0.3)

ax6 = fig.add_subplot(2, 3, 6)
tt = np.linspace(0, np.pi / 4, 200)
ax6.plot(tt, binary_entropy(np.cos(tt) ** 2), "b-", label="entanglement entropy (exact)")
ax6.plot(th_grid, ER[:, -1], "ro", label=r"$E_R$ numerical ($p=1$)")
ax6.plot(th_grid, RG[:, -1], "gs", label=r"$R_G$ numerical ($p=1$)")
ax6.plot(tt, np.sin(2 * tt), "g-", alpha=0.5, label=r"$\sin 2\theta$")
ax6.set_xlabel(r"$\theta$"); ax6.set_ylabel("entanglement")
ax6.set_title("(f) Pure states"); ax6.legend(); ax6.grid(alpha=0.3)

plt.tight_layout()
plt.show()

5. Code Walkthrough

Basic tools

werner(p) and family(p, theta) build the $4\times4$ density matrices. partial_transpose reshapes $\rho$ into a rank-4 tensor $\rho_{i j, k l}\to\rho_{i,l,,k,j}$ by swapping the two indices of the second qubit, which is the partial transpose $T_B$. negativity sums the absolute values of the negative eigenvalues of $\rho^{T_B}$. binary_entropy is clipped away from $0$ and $1$ so that the $\log$ never produces nan.

Parametrizing separable states

build_sigma maps a flat parameter vector $\mathbf x\in\mathbb R^{5K}$ to a separable density matrix:

The first $K$ entries are logits. A softmax (shifted by the maximum for numerical stability) turns them into probabilities $p_k$, so positivity and normalization hold automatically and the problem becomes unconstrained.
The remaining $4K$ entries are Bloch-sphere angles $(\vartheta_A,\varphi_A,\vartheta_B,\varphi_B)$ for each term.
The vectors $|a_k\rangle$ and $|b_k\rangle$ are built for all $k$ at once. The batched outer product A[:, :, None] * B[:, None, :] gives all $K$ product vectors $|a_kb_k\rangle$ in one operation.
A single einsum assembles $\sigma=\sum_k p_k|\psi_k\rangle\langle\psi_k|$.

Because every $\sigma$ produced this way is separable, the optimizer can never wander outside $\mathrm{SEP}$.

Objective function

Since $S(\rho|\sigma)=\mathrm{Tr},\rho\log_2\rho-\mathrm{Tr},\rho\log_2\sigma$, the first term is a constant, computed once in make_objective. Each evaluation only diagonalizes $\sigma$ once: eigh gives $\sigma=V,\mathrm{diag}(w),V^\dagger$, and $\log_2\sigma=V,\mathrm{diag}(\log_2 w),V^\dagger$. The eigenvalues are clipped at $10^{-13}$ so that a numerically zero eigenvalue cannot produce -inf. The term (V * lw) @ V.conj().T multiplies column-wise by lw, which avoids building a diagonal matrix.

Optimization and speed-ups

The landscape is nonconvex in $\mathbf x$, so relative_entropy_of_entanglement runs L-BFGS-B from several random starts and keeps the best. The code is already in its accelerated form, and four design choices keep the runtime to roughly a minute or so on a standard CPU:

The objective has no Python loops over the $K$ product states (batched arrays and einsum).
Only one $4\times4$ eigendecomposition is needed per evaluation.
On the $(p,\theta)$ grid, the optimum found for the previous $p$ is passed as a warm start (x0=x_prev) for the next one, so only 2 random restarts are needed there.
The robustness problem is a tiny convex SDP, solved in about 10–20 ms per state.

Robustness SDP

generalized_robustness is a direct translation of the SDP above. Because all of our states are real, the optimal $Y$ can be taken to be real symmetric (averaging a solution with its complex conjugate keeps it feasible and optimal), so symmetric=True suffices and keeps the problem small. The constraint Y - rho_r >> 0 encodes $Y\succeq\rho$, and cp.partial_transpose(Y, dims=[2, 2], axis=1) >> 0 encodes $Y^{T_B}\succeq0$. The objective $\mathrm{Tr},Y-1$ equals $s$.

Experiments

Part (3) sweeps 13 values of the Werner parameter $p$ and compares numerical results with the exact formulas $E_R=1-h(F)$ and $R_G=2F-1$. Part (4) evaluates both measures on an $11\times 9$ grid in $(p,\theta)$, and checks the pure-state values at $p=1$ and the inequality $E_R\le\log_2(1+R_G)$.

Plotting

All six panels live in one figure. Panels (b) and (c) are 3D surfaces created with projection="3d". The legends and axis labels use TeX syntax in raw strings.

6. Execution Results

=== Werner states: rho = p|Phi+>＜Phi+| + (1-p) I/4 ===
     p       F   E_R(num)  E_R(exact)   R_G(num)  R_G(exact)       N
 0.000  0.2500   0.000000    0.000000   0.000000    0.000000  0.0000
 0.083  0.3125   0.000000    0.000000   0.000000    0.000000  0.0000
 0.167  0.3750   0.000000    0.000000   0.000000    0.000000  0.0000
 0.250  0.4375   0.000000    0.000000   0.000000    0.000000  0.0000
 0.333  0.5000   0.000000    0.000000   0.000000    0.000000  0.0000
 0.417  0.5625   0.011301    0.011301   0.125000    0.125000  0.0625
 0.500  0.6250   0.045566    0.045566   0.250000    0.250000  0.1250
 0.583  0.6875   0.103962    0.103962   0.375000    0.375000  0.1875
 0.667  0.7500   0.188722    0.188722   0.500000    0.500000  0.2500
 0.750  0.8125   0.303788    0.303788   0.625000    0.625000  0.3125
 0.833  0.8750   0.456436    0.456436   0.750000    0.750000  0.3750
 0.917  0.9375   0.662710    0.662710   0.875000    0.875000  0.4375
 1.000  1.0000   1.000000    1.000000   1.000000    1.000000  0.5000
max |E_R(num) - E_R(exact)| = 6.45e-08
max |R_G(num) - R_G(exact)| = 3.37e-09

=== Pure states (p = 1): E_R vs. entanglement entropy ===
theta = 0.0000   E_R(num) = 0.000000   entropy = 0.000000
theta = 0.0982   E_R(num) = 0.078179   entropy = 0.078179
theta = 0.1963   E_R(num) = 0.233327   entropy = 0.233327
theta = 0.2945   E_R(num) = 0.417032   entropy = 0.417032
theta = 0.3927   E_R(num) = 0.600876   entropy = 0.600876
theta = 0.4909   E_R(num) = 0.764191   entropy = 0.764191
theta = 0.5890   E_R(num) = 0.891619   entropy = 0.891619
theta = 0.6872   E_R(num) = 0.972369   entropy = 0.972368
theta = 0.7854   E_R(num) = 1.000000   entropy = 1.000000

max over grid of [E_R - log2(1+R_G)] = 6.28e-08  (<= 0 up to numerical tolerance)
Total computation time: 157.6 s

7. Discussion of the Results

Console output

The Werner table shows that the optimization over separable states reproduces the exact result. In my run the maximum deviation was about $6\times10^{-8}$ for $E_R$ and about $3\times10^{-9}$ for $R_G$. Some concrete values:

$p$	$F$	$E_R$	$R_G$
$0.333$	$0.5000$	$0$	$0$
$0.500$	$0.6250$	$0.045566$	$0.2500$
$0.750$	$0.8125$	$0.303788$	$0.6250$
$1.000$	$1.0000$	$1.0000$	$1.0000$

Both measures vanish for $p\le 1/3$ and become positive immediately afterward, which is the Peres–Horodecki threshold $F=1/2$. The pure-state table shows that $E_R$ matches the entanglement entropy $h(\cos^2\theta)$ to about six digits at every $\theta$. Finally, $\max[E_R-\log_2(1+R_G)]$ is of order $10^{-8}$, which is zero up to numerical tolerance (it is not exactly nonpositive because the optimizer stops at a finite tolerance). The inequality holds everywhere, with equality at the maximally entangled state.

Panel (a): Werner states

Red dots ($E_R$ from the optimizer) lie on the blue curve ($1-h(F)$). Black squares ($R_G$ from the SDP) lie on the green line $2F-1$. Both are zero up to the dotted line at $p=1/3$. Note how different the shapes are: $R_G$ grows linearly in $p$ while $E_R$ grows slowly at first, like $(F-\tfrac12)^2$ near the threshold, and only catches up at $p=1$. The negativity (magenta) is exactly half of $R_G$ for these states, so it is a rescaled copy of the robustness. Different entanglement measures agree on which states are entangled but not on how much.

Panels (b) and (c): 3D surfaces

Both surfaces are flat at zero over a large region, the separable states, and rise toward the corner $p=1$, $\theta=\pi/4$, where the state is the maximally entangled Bell state. At that corner $E_R=R_G=1$. Both surfaces are monotonic in $p$ and in $\theta$, as expected: more purity and a more balanced Schmidt decomposition give more entanglement. The robustness surface (c) is noticeably more “tent-like” and rises earlier, whereas the relative-entropy surface (b) stays low over most of the domain and shoots up near $p\to1$. This is the same linear-versus-quadratic difference seen in panel (a).

Panel (d): the entanglement map

This is a top-down view of surface (b). The dashed red curve is the analytic boundary $p=1/(1+2\sin2\theta)$, which separates separable states (below and to the left) from entangled ones. The boundary passes through $p=1/3$ at $\theta=\pi/4$ (the Werner case) and reaches $p=1$ at $\theta=0$, where the pure state is a product state. The region where $E_R$ is exactly zero coincides with this analytic boundary, which confirms that the numerical optimizer correctly detects the transition.

Panel (e): the inequality $E_R\le\log_2(1+R_G)$

Each dot is one grid state, coloured by $p$. All dots lie on or below the dashed diagonal, which is the bound $E_R\le\log_2(1+R_G)$. The points touch the diagonal only at the origin (separable states) and at the Bell state (top right). Mixed states and partially entangled pure states sit strictly below it. The gap shows that the two measures encode genuinely different information about the state.

Panel (f): pure states

For pure states, the numerical $E_R$ (red dots) matches the entanglement entropy $h(\cos^2\theta)$ (blue curve), as it must, since for pure states the closest separable state is the dephased state. The numerical $R_G$ (green squares) lies on $\sin2\theta$, in agreement with the closed form $(\sum_i\sqrt{\lambda_i})^2-1$ for pure states with Schmidt coefficients $\lambda_i$. Here $\sin 2\theta\ge h(\cos^2\theta)$ for all $\theta$, again consistent with $E_R\le\log_2(1+R_G)$.

8. Summary

Distance-based entanglement measures reduce to optimization problems: a nonconvex one for the relative entropy and a convex SDP for the robustness.
Parametrizing separable states from the inside (softmax weights times Bloch-sphere product states) turns the constrained problem into an unconstrained one that L-BFGS-B handles quickly.
Vectorization, a single eigendecomposition per evaluation, and warm starts keep the full 3D study down to about a minute or so.
The numerical results match the exact Werner, pure-state and PPT-boundary formulas to $10^{-7}$ or better, and they confirm $E_R\le\log_2(1+R_G)$.

The same approach extends to other measures, such as the geometric measure, by changing the objective, and to higher-dimensional systems, where PPT is only a relaxation of separability and the inner parametrization used here becomes the more reliable route.

Optimizing Quantum Entanglement Witnesses

October 10, 2026

Finding the Best Observable to Detect a Given Entangled State

Entanglement is a resource, but how do you certify that a lab-prepared state is actually entangled? Full state tomography is expensive. An entanglement witness is a cheaper alternative: a single observable $W$ whose expectation value is non-negative on every separable state and negative on at least one entangled state.

$$
\mathrm{Tr}(W\sigma)\ge 0\quad \forall,\sigma\in\mathrm{SEP},\qquad \mathrm{Tr}(W\rho)<0 ;\Rightarrow; \rho \text{ is entangled.}
$$

Measuring one number is enough to certify entanglement. The catch is that a witness built naively from the state you intended to prepare can fail on the noisy state you actually have. In this article we optimize the witness for a given entangled state and show that the intuitive choice is not the best one.

1. The Problem Setup

We consider two qubits and a family of pure states with a Schmidt angle $\theta$ and a relative phase $\varphi$:

$$
|\psi(\theta,\varphi)\rangle=\cos\theta,|00\rangle+e^{i\varphi}\sin\theta,|11\rangle .
$$

The state we want to certify is this pure state mixed with white noise:

$$
\rho(p)=p,|\psi_0\rangle\langle\psi_0|+(1-p),\frac{\mathbb{I}_4}{4},\qquad \theta_0=\frac{\pi}{6},;; \varphi_0=\frac{\pi}{2},;; p=0.6 .
$$

Witness family

For every $(\theta,\varphi)$ define

This is a valid witness. For any separable state $\sigma$, the overlap with a pure state is bounded by the largest squared Schmidt coefficient, $\langle\psi|\sigma|\psi\rangle\le\lambda_{\max}^2$, hence $\mathrm{Tr}(W\sigma)\ge0$.

Fair comparison through normalization

Different witnesses have different scales, so we normalize each one by its trace:

$$
\tilde W=\frac{W}{\mathrm{Tr},W},\qquad \mathrm{Tr},W=4\lambda_{\max}^2-1 .
$$

The optimization problem

$$
\min_{\theta,\varphi};\mathrm{Tr}!\left[\tilde W(\theta,\varphi),\rho\right]
=\min_{\theta,\varphi};\frac{\lambda_{\max}^2(\theta)-p,|\langle\psi(\theta,\varphi)|\psi_0\rangle|^2-\frac{1-p}{4}}{4\lambda_{\max}^2(\theta)-1}.
$$

The overlap has a closed form, which lets us evaluate the whole landscape quickly:

$$
\langle\psi(\theta,\varphi)|\psi_0\rangle=\cos\theta\cos\theta_0+\sin\theta\sin\theta_0,e^{i(\varphi_0-\varphi)} .
$$

A reference answer from the PPT criterion

For two qubits, a state is entangled if and only if its partial transpose $\rho^{T_B}$ has a negative eigenvalue. For our state,

$$
\lambda_{\min}!\left(\rho^{T_B}\right)=\frac{1-p}{4}-\frac{p}{2}\sin 2\theta_0 ,
$$

which becomes negative exactly when

$$
p>\frac{1}{1+2\sin 2\theta_0}=\frac{1}{1+\sqrt3}\approx 0.366 .
$$

The witness $|v\rangle\langle v|^{T_B}$ built from the negative eigenvector $|v\rangle$ attains this value, so it is the benchmark any witness can be compared against.

2. The Complete Python Code

Everything is in a single cell. The heavy 3D landscape is computed with the vectorized closed form instead of a double loop over matrices, and a consistency check confirms that it agrees with the explicit matrix calculation.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from scipy.optimize import minimize

I4 = np.eye(4, dtype=complex)

def psi_vec(theta, phi):
    v = np.zeros(4, dtype=complex)
    v[0] = np.cos(theta)
    v[3] = np.exp(1j * phi) * np.sin(theta)
    return v

def witness(theta, phi):
    v = psi_vec(theta, phi)
    a = max(np.cos(theta) ** 2, np.sin(theta) ** 2)
    return a * I4 - np.outer(v, v.conj())

def normalized_value(theta, phi, rho):
    W = witness(theta, phi)
    return float(np.real(np.trace(W @ rho)) / np.real(np.trace(W)))

def target_state(p, theta0, phi0):
    v = psi_vec(theta0, phi0)
    return p * np.outer(v, v.conj()) + (1 - p) * I4 / 4

def partial_transpose(rho):
    r = rho.reshape(2, 2, 2, 2)
    return r.transpose(0, 3, 2, 1).reshape(4, 4)

def ppt_min_eigenvalue(rho):
    return float(np.linalg.eigvalsh(partial_transpose(rho))[0])

def optimize_witness(rho):
    best = None
    for t0 in np.linspace(0.2, np.pi / 2 - 0.2, 3):
        for f0 in np.linspace(0, 2 * np.pi, 4, endpoint=False):
            res = minimize(lambda x: normalized_value(x[0], x[1], rho), [t0, f0],
                           method="Nelder-Mead",
                           options={"xatol": 1e-8, "fatol": 1e-12, "maxiter": 2000})
            if best is None or res.fun < best.fun:
                best = res
    return best

def value_grid(TH, PH, p, theta0, phi0):
    a = np.maximum(np.cos(TH) ** 2, np.sin(TH) ** 2)
    ov = np.cos(TH) * np.cos(theta0) + np.sin(TH) * np.sin(theta0) * np.exp(1j * (phi0 - PH))
    return (a - p * np.abs(ov) ** 2 - (1 - p) / 4) / (4 * a - 1)

theta0, phi0, p = np.pi / 6, np.pi / 2, 0.6
rho = target_state(p, theta0, phi0)

res = optimize_witness(rho)
th_opt, ph_opt = res.x[0] % np.pi, res.x[1] % (2 * np.pi)
val_target = normalized_value(theta0, phi0, rho)
val_bell = normalized_value(np.pi / 4, phi0, rho)
val_opt = res.fun
val_ppt = ppt_min_eigenvalue(rho)

print("State: p = %.2f, theta0 = pi/6, phi0 = pi/2" % p)
print("Tr(W rho) with the target-state witness  : %+.6f" % val_target)
print("Tr(W rho) with the Bell-type witness     : %+.6f" % val_bell)
print("Tr(W rho) with the optimized witness     : %+.6f" % val_opt)
print("Min eigenvalue of rho^(T_B) (PPT bound)  : %+.6f" % val_ppt)
print("Optimal theta = %.6f (pi/4 = %.6f)" % (th_opt, np.pi / 4))
print("Optimal phi   = %.6f (phi0 = %.6f)" % (ph_opt, phi0))

th = np.linspace(0.001, np.pi / 2 - 0.001, 120)
ph = np.linspace(0, 2 * np.pi, 180)
TH, PH = np.meshgrid(th, ph)
Z = value_grid(TH, PH, p, theta0, phi0)
chk = normalized_value(0.7, 1.3, rho) - float(value_grid(np.array(0.7), np.array(1.3), p, theta0, phi0))
print("Consistency check (matrix vs. analytic): %.2e" % abs(chk))

ps = np.linspace(0.2, 1.0, 33)
v_target, v_bell, v_opt = [], [], []
for pp in ps:
    r = target_state(pp, theta0, phi0)
    v_target.append(normalized_value(theta0, phi0, r))
    v_bell.append(normalized_value(np.pi / 4, phi0, r))
    v_opt.append(optimize_witness(r).fun)

def crossing(x, y):
    for i in range(len(x) - 1):
        if y[i] >= 0 > y[i + 1]:
            return x[i] + (0 - y[i]) * (x[i + 1] - x[i]) / (y[i + 1] - y[i])
    return float("nan")

print("Detection threshold, target-state witness : p > %.3f" % crossing(ps, v_target))
print("Detection threshold, Bell-type witness    : p > %.3f" % crossing(ps, v_bell))
print("Detection threshold, optimized witness    : p > %.3f" % crossing(ps, v_opt))
print("Exact entanglement threshold (PPT)        : p > %.3f" % (1 / (1 + 2 * np.sin(2 * theta0))))

fig = plt.figure(figsize=(16, 12))

ax1 = fig.add_subplot(2, 2, 1, projection="3d")
surf = ax1.plot_surface(TH, PH, Z, cmap=cm.coolwarm, linewidth=0, antialiased=True, alpha=0.95)
ax1.scatter([th_opt], [ph_opt], [val_opt], color="black", s=80, label="Optimum")
ax1.scatter([theta0], [phi0], [val_target], color="lime", s=80, edgecolor="black", label="Target-state witness")
ax1.set_xlabel(r"$\theta$")
ax1.set_ylabel(r"$\varphi$")
ax1.set_zlabel(r"$\mathrm{Tr}(\tilde{W}\rho)$")
ax1.set_title("Witness landscape")
ax1.view_init(elev=28, azim=-60)
ax1.legend(loc="upper left")
fig.colorbar(surf, ax=ax1, shrink=0.6, pad=0.1)

ax2 = fig.add_subplot(2, 2, 2)
cs = ax2.contourf(TH, PH, Z, levels=30, cmap=cm.coolwarm)
ax2.contour(TH, PH, Z, levels=[0.0], colors="black", linewidths=2)
ax2.scatter([th_opt], [ph_opt], color="black", s=80, label="Optimal witness")
ax2.scatter([theta0], [phi0], color="lime", s=80, edgecolor="black", label="Target-state witness")
ax2.set_xlabel(r"$\theta$")
ax2.set_ylabel(r"$\varphi$")
ax2.set_title("Contour map (black line: Tr = 0, negative inside)")
ax2.legend(loc="upper right")
fig.colorbar(cs, ax=ax2)

ax3 = fig.add_subplot(2, 2, 3)
ax3.plot(ps, v_target, "o-", label="Target-state witness", color="tab:green")
ax3.plot(ps, v_bell, "s-", label="Bell-type witness", color="tab:orange")
ax3.plot(ps, v_opt, "^-", label="Optimized witness", color="tab:blue")
ax3.axhline(0, color="black", linewidth=1)
ax3.axvline(1 / (1 + 2 * np.sin(2 * theta0)), color="gray", linestyle="--", label="PPT threshold")
ax3.set_xlabel("Mixing parameter p")
ax3.set_ylabel(r"$\mathrm{Tr}(\tilde{W}\rho)$")
ax3.set_title("Detection power versus noise")
ax3.legend()
ax3.grid(alpha=0.3)

ax4 = fig.add_subplot(2, 2, 4)
for pp, c in zip([0.4, 0.5, 0.6, 0.8, 1.0], cm.viridis(np.linspace(0, 1, 5))):
    ax4.plot(th, value_grid(th, np.full_like(th, phi0), pp, theta0, phi0), color=c, label="p = %.1f" % pp)
ax4.axhline(0, color="black", linewidth=1)
ax4.axvline(np.pi / 4, color="gray", linestyle="--")
ax4.axvline(theta0, color="green", linestyle=":")
ax4.set_xlabel(r"$\theta$  (at $\varphi=\varphi_0$)")
ax4.set_ylabel(r"$\mathrm{Tr}(\tilde{W}\rho)$")
ax4.set_title(r"Slice along $\theta$: dashed = $\pi/4$, dotted = $\theta_0$")
ax4.legend()
ax4.grid(alpha=0.3)

plt.tight_layout()
plt.show()

3. Detailed Code Walkthrough

Building blocks

psi_vec(theta, phi) returns the four-component vector $\cos\theta,|00\rangle+e^{i\varphi}\sin\theta,|11\rangle$ in the basis ${|00\rangle,|01\rangle,|10\rangle,|11\rangle}$. Index 0 is $|00\rangle$ and index 3 is $|11\rangle$.
witness(theta, phi) implements $W=\lambda_{\max}^2\mathbb{I}-|\psi\rangle\langle\psi|$. The constant a is the larger of $\cos^2\theta$ and $\sin^2\theta$, which guarantees that no separable state gives a negative value.
normalized_value computes $\mathrm{Tr}(W\rho)/\mathrm{Tr}(W)$. Dividing by the trace puts every witness on the same scale, so the comparison is fair.
target_state builds the noisy state $\rho(p)$ as a $4\times4$ density matrix.

Partial transpose and the PPT benchmark

partial_transpose reshapes the matrix into a tensor $\rho_{ab,a’b’}$ with shape (2,2,2,2), swaps the two indices of the second qubit with transpose(0, 3, 2, 1), and reshapes it back. ppt_min_eigenvalue then returns the smallest eigenvalue of $\rho^{T_B}$. Because $\mathrm{Tr}(|v\rangle\langle v|^{T_B}\rho)=\langle v|\rho^{T_B}|v\rangle$, this eigenvalue equals the value attained by the optimal PPT-type witness, and it serves as the ground truth.

Optimization

optimize_witness minimizes $\mathrm{Tr}(\tilde W\rho)$ over $(\theta,\varphi)$ with Nelder–Mead. The landscape is periodic in $\varphi$ and has a kink at $\theta=\pi/4$, where $\lambda_{\max}^2$ switches branches. A single starting point could get stuck, so the function launches $3\times4=12$ starting points and keeps the best result. Nelder–Mead needs no gradients, which suits the non-smooth point at $\theta=\pi/4$.

Speed-up for the landscape

Evaluating the witness matrix at every point of a $120\times180$ grid would mean 21,600 matrix constructions. value_grid instead uses the analytic overlap

$$
\mathrm{Tr}(W\rho)=\lambda_{\max}^2-p,|\langle\psi|\psi_0\rangle|^2-\frac{1-p}{4},
$$

evaluated for the whole mesh by NumPy broadcasting. The line that prints the “Consistency check” compares this closed form against the explicit matrix version at one point and should show a value near zero.

Noise sweep

The loop over ps repeats the optimization for 33 noise levels and records three curves: the witness built from the target state, the Bell-type witness $(\theta=\pi/4)$, and the numerically optimized witness. The helper crossing finds the value of $p$ at which each curve crosses zero by linear interpolation. That crossing is the detection threshold of that witness.

4. Console Output

State: p = 0.60, theta0 = pi/6, phi0 = pi/2
Tr(W rho) with the target-state witness  : +0.025000
Tr(W rho) with the Bell-type witness     : -0.159808
Tr(W rho) with the optimized witness     : -0.159808
Min eigenvalue of rho^(T_B) (PPT bound)  : -0.159808
Optimal theta = 0.785398 (pi/4 = 0.785398)
Optimal phi   = 1.570796 (phi0 = 1.570796)
Consistency check (matrix vs. analytic): 0.00e+00
Detection threshold, target-state witness : p > 0.667
Detection threshold, Bell-type witness    : p > 0.366
Detection threshold, optimized witness    : p > 0.366
Exact entanglement threshold (PPT)        : p > 0.366

What the numbers mean

For the state with $p=0.6$, the target-state witness uses $\lambda_{\max}^2=\cos^2\theta_0=0.75$ and $\mathrm{Tr},W=2$. Its value is

$$
\mathrm{Tr}(\tilde W\rho)=\frac{0.75-\left[p+\frac{1-p}{4}\right]}{2}=\frac{0.75-0.70}{2}=+0.025>0 ,
$$

so it fails to detect the entanglement, even though the state is entangled. The Bell-type witness at $\theta=\pi/4$ has $\lambda_{\max}^2=1/2$ and $\mathrm{Tr},W=1$, and gives

$$
\mathrm{Tr}(\tilde W\rho)=\frac12-\left[p,\frac{1+\sin 2\theta_0}{2}+\frac{1-p}{4}\right]\approx-0.1598<0 ,
$$

which is a clear violation. The optimizer lands on $\theta^\ast=\pi/4$, $\varphi^\ast=\varphi_0$, and its value coincides with the PPT lower bound

$$
\lambda_{\min}!\left(\rho^{T_B}\right)=\frac{1-p}{4}-\frac{p}{2}\sin 2\theta_0\approx-0.1598 .
$$

The three detection thresholds printed at the end confirm the same picture. The target-state witness needs $p>2/3$, while the optimized witness reaches the exact entanglement threshold $p>1/(1+\sqrt3)\approx0.366$.

5. Graphs

All four panels are produced by one figure.

Top left: the 3D witness landscape

The surface shows $\mathrm{Tr}(\tilde W\rho)$ over the witness parameters $(\theta,\varphi)$. Blue regions are negative, meaning the witness detects entanglement. Red regions are positive, meaning it does not. The green point is the witness built from the target state, $(\theta_0,\varphi_0)$, and it lies in the positive region. The black point is the true minimum, which sits at the deepest part of the blue valley. The valley is centered on the phase $\varphi=\varphi_0$, so the phase must match the state exactly, while the Schmidt angle should be pushed towards the maximally entangled value $\pi/4$. The kink along $\theta=\pi/4$ comes from the $\max(\cdot)$ in $\lambda_{\max}^2$.

Top right: the contour map

This is the top view of the same landscape. The black line marks $\mathrm{Tr}=0$, and everything inside it is a successful witness. The green point lies outside the contour, so the intuitive witness fails. The black point lies near the center of the negative region. The elongated shape of the contour shows that detection is sensitive to the phase $\varphi$ and tolerant to moderate changes in $\theta$ around $\pi/4$.

Bottom left: detection power versus noise

Each curve shows the witness value as the noise level changes. A witness certifies entanglement wherever its curve is below zero. The green curve (target-state witness) crosses zero at $p\approx0.667$, the orange curve (Bell-type witness) at $p\approx0.366$, and the blue curve (optimized witness) also at $p\approx0.366$. The dashed vertical line marks the exact PPT threshold. The optimized witness reaches this limit, which no witness can beat for two qubits. The target-state witness wastes almost half of the detectable range.

Bottom right: slices along $\theta$

Here the phase is fixed at $\varphi_0$ and $\theta$ is scanned for several noise strengths. Every curve has its minimum at the dashed line $\theta=\pi/4$, not at the dotted line $\theta_0$. As $p$ decreases, the curves rise and the dip below zero shrinks. Near $p=0.4$ the minimum is only just negative, and it reaches zero at the threshold $p\approx0.366$.

6. Conclusion

A witness tailored to the state you intended to prepare is not necessarily the best witness for the state you actually hold. For the weakly entangled state $\theta_0=\pi/6$ mixed with white noise, the target-state witness detects entanglement only for $p>2/3$. The optimized witness, with $\theta=\pi/4$ and the phase matched to $\varphi_0$, detects it all the way down to $p\approx0.366$, which is the exact boundary given by the PPT criterion.

The same strategy extends to larger systems. One parametrizes a family of valid witnesses, evaluates $\mathrm{Tr}(\tilde W\rho)$ for the state at hand, and minimizes over the family. In higher dimensions, where the PPT criterion is no longer sufficient, this optimization over witnesses becomes the practical tool for certifying entanglement.

Testing Separability with PPT and Symmetric Extension Hierarchies in Python

October 9, 2026

Deciding whether a quantum state is separable or entangled is a hard problem in general. Two practical relaxations get us surprisingly far: the PPT criterion and the symmetric extension hierarchy. This post tests both on a concrete two-qubit example and visualizes how the relaxations behave.

The Problem

A bipartite state $\rho_{AB}$ is separable if it can be written as

$$
\rho_{AB}=\sum_i q_i,\rho_A^{(i)}\otimes\rho_B^{(i)},\qquad q_i\ge 0,\ \sum_i q_i=1 .
$$

Checking this directly is NP-hard, so we use necessary conditions that can be tested with semidefinite programming (SDP).

PPT relaxation. Every separable state has a positive semidefinite partial transpose:

$$
\rho_{AB}\ \text{separable}\ \Rightarrow\ \rho_{AB}^{T_B}\succeq 0 .
$$

Symmetric extension relaxation. A state $\rho_{AB}$ is $k$-extendible if there is a state $\sigma_{AB_1\cdots B_k}$ such that

Separable states are $k$-extendible for every $k$, and conversely, a state that is $k$-extendible for all $k$ is separable. Adding PPT conditions on $\sigma$ gives the stronger PPT symmetric extension test of Doherty, Parrilo and Spedalieri:

$$
\sigma^{T_{B_{k-m+1}\cdots B_k}}\succeq 0,\qquad m=1,\dots,k .
$$

The Example

We use a family of noisy two-qubit entangled states:

$$
\rho(p,\theta)=p,|\Phi_\theta\rangle\langle\Phi_\theta|+(1-p),\frac{I}{4},\qquad
|\Phi_\theta\rangle=\cos\theta,|00\rangle+\sin\theta,|11\rangle .
$$

For $\theta=\pi/4$ this is a Werner-type state built from a Bell state. The partial transpose has the smallest eigenvalue

$$
\lambda_{\min}!\left(\rho^{T_B}\right)=\frac{1-p}{4}-\frac{p,\sin 2\theta}{2},
$$

so for two qubits, where PPT is both necessary and sufficient, the exact separability boundary is

$$
p^\ast(\theta)=\frac{1}{1+2\sin 2\theta},\qquad p^\ast(\pi/4)=\frac13 .
$$

For each $\theta$ and each hierarchy level $k$, we solve the SDP

$$
\max\ p\quad\text{s.t.}\quad \sigma\succeq0,\ \operatorname{Tr}\sigma=1,\ \sigma\ \text{symmetric on }B,\ \operatorname{Tr}_{B_2\cdots B_k}\sigma=p,|\Phi_\theta\rangle\langle\Phi_\theta|+(1-p)\tfrac{I}{4}.
$$

This gives the largest noise-free weight $p$ that still passes the test. The true boundary is $p^\ast$, so a relaxation is tight when it reproduces it.

Source Code

import time
import warnings
warnings.filterwarnings("ignore")
import numpy as np
import cvxpy as cp
import matplotlib.pyplot as plt

# ------------------------------------------------------------
# 1. States and the PPT criterion
# ------------------------------------------------------------
I4 = np.eye(4) / 4.0

def phi(theta):
    v = np.zeros(4)
    v[0], v[3] = np.cos(theta), np.sin(theta)
    return v

def rho_state(p, theta):
    v = phi(theta)
    return p * np.outer(v, v) + (1.0 - p) * I4

def partial_transpose_B(rho):
    return rho.reshape(2, 2, 2, 2).transpose(0, 3, 2, 1).reshape(4, 4)

def ppt_min_eig(p, theta):
    return np.linalg.eigvalsh(partial_transpose_B(rho_state(p, theta)))[0]

def ppt_threshold(theta):
    return 1.0 / (1.0 + 2.0 * np.sin(2.0 * theta))

# ------------------------------------------------------------
# 2. Symmetric extension (and PPT symmetric extension) as an SDP
# ------------------------------------------------------------
def swap_matrix(n, i, j):
    d = 2 ** n
    P = np.zeros((d, d))
    for idx in range(d):
        bits = [(idx >> (n - 1 - q)) & 1 for q in range(n)]
        bits[i], bits[j] = bits[j], bits[i]
        new = 0
        for b in bits:
            new = (new << 1) | b
        P[new, idx] = 1.0
    return P

def trace_out_last(X, n):
    d = 2 ** (n - 1)
    out = 0
    for i in range(2):
        e = np.zeros((2, 1)); e[i, 0] = 1.0
        K = np.kron(np.eye(d), e.T)
        out = out + K @ X @ K.T
    return out

def transpose_qubit(X, q, n):
    out = 0
    for i in range(2):
        for j in range(2):
            E = np.zeros((2, 2)); E[i, j] = 1.0
            K = np.kron(np.kron(np.eye(2 ** q), E), np.eye(2 ** (n - q - 1)))
            out = out + K @ X @ K
    return out

def solve(prob):
    try:
        prob.solve(solver=cp.CLARABEL)
    except Exception:
        prob.solve(solver=cp.SCS, eps=1e-7, max_iters=100000)

def max_visibility(theta, k, ppt=False):
    n = k + 1
    d = 2 ** n
    v = phi(theta)
    target_pure = np.outer(v, v)
    sigma = cp.Variable((d, d), symmetric=True)
    p = cp.Variable()
    cons = [sigma >> 0, cp.trace(sigma) == 1, p >= 0, p <= 1]
    for j in range(1, k):
        P = swap_matrix(n, j, j + 1)
        cons.append(P @ sigma @ P.T == sigma)
    marg = sigma
    for m in range(n, 2, -1):
        marg = trace_out_last(marg, m)
    cons.append(marg == p * target_pure + (1 - p) * I4)
    if ppt:
        for m in range(1, k + 1):
            Y = sigma
            for q in range(n - m, n):
                Y = transpose_qubit(Y, q, n)
            Y = (Y + Y.T) / 2
            cons.append(Y >> 0)
    prob = cp.Problem(cp.Maximize(p), cons)
    solve(prob)
    return float(p.value)

# ------------------------------------------------------------
# 3. Example A: Werner-type state (theta = pi/4)
# ------------------------------------------------------------
print("=== Example A: Werner-type state, rho(p) = p|Phi+><Phi+| + (1-p) I/4 ===")
for p in [0.2, 1/3, 0.5, 0.8]:
    ev = np.linalg.eigvalsh(partial_transpose_B(rho_state(p, np.pi / 4)))
    verdict = "PPT (separable)" if ev[0] >= -1e-12 else "NPT (entangled)"
    print(f"p = {p:.4f}  min eig of rho^TB = {ev[0]:+.5f}  -> {verdict}")

ks = [1, 2, 3, 4]
t0 = time.time()
werner_sym = [max_visibility(np.pi / 4, k) for k in ks]
werner_dps = [max_visibility(np.pi / 4, k, ppt=True) for k in ks]
print(f"\nElapsed: {time.time() - t0:.1f} s")
print("\n k | sym. extension p_max | (k+2)/(3k) | PPT sym. ext. p_max")
for k, a, b in zip(ks, werner_sym, werner_dps):
    print(f" {k} | {a:20.5f} | {(k + 2) / (3 * k):10.5f} | {b:19.5f}")
print(f"\nPPT threshold (exact separability limit): {ppt_threshold(np.pi / 4):.5f}")

# ------------------------------------------------------------
# 4. Example B: non-symmetric family over theta
# ------------------------------------------------------------
thetas = np.linspace(np.pi / 16, np.pi / 4, 8)
kk = [2, 3]
sym_th = {k: [max_visibility(t, k) for t in thetas] for k in kk}
dps_th = {k: [max_visibility(t, k, ppt=True) for t in thetas] for k in kk}
ppt_th = np.array([ppt_threshold(t) for t in thetas])

print("\n=== Example B: rho(p, theta) = p|Phi(theta)><Phi(theta)| + (1-p) I/4 ===")
print(" theta/pi | PPT limit | sym k=2 | sym k=3 | PPT-sym k=2 | PPT-sym k=3")
for i, t in enumerate(thetas):
    print(f" {t/np.pi:8.4f} | {ppt_th[i]:9.5f} | {sym_th[2][i]:7.5f} | {sym_th[3][i]:7.5f}"
          f" | {dps_th[2][i]:11.5f} | {dps_th[3][i]:11.5f}")

# ------------------------------------------------------------
# 5. Visualization (single figure)
# ------------------------------------------------------------
fig = plt.figure(figsize=(18, 13))

ax1 = fig.add_subplot(2, 2, 1, projection="3d")
th_g = np.linspace(0.01, np.pi / 4, 60)
p_g = np.linspace(0, 1, 60)
TH, PG = np.meshgrid(th_g, p_g)
Z = np.vectorize(ppt_min_eig)(PG, TH)
ax1.plot_surface(TH / np.pi, PG, Z, cmap="coolwarm", alpha=0.9, linewidth=0)
ax1.plot_surface(TH / np.pi, PG, np.zeros_like(Z), color="gray", alpha=0.25)
tc = np.linspace(0.01, np.pi / 4, 100)
ax1.plot(tc / np.pi, ppt_threshold(tc), np.zeros_like(tc), "k-", lw=3, label="separability boundary")
ax1.set_xlabel(r"$\theta/\pi$"); ax1.set_ylabel(r"$p$"); ax1.set_zlabel(r"$\lambda_{\min}(\rho^{T_B})$")
ax1.set_title("PPT test: minimum eigenvalue of the partial transpose")
ax1.view_init(elev=25, azim=-130)
ax1.legend()

ax2 = fig.add_subplot(2, 2, 2)
ax2.plot(thetas / np.pi, ppt_th, "k-o", lw=2.5, label="PPT limit (exact)")
colors = {2: "tab:blue", 3: "tab:red"}
for k in kk:
    ax2.plot(thetas / np.pi, sym_th[k], "--s", color=colors[k], label=f"{k}-symmetric extension")
    ax2.plot(thetas / np.pi, dps_th[k], ":^", color=colors[k], label=f"PPT {k}-symmetric extension")
ax2.fill_between(thetas / np.pi, 0, ppt_th, color="green", alpha=0.12)
ax2.set_xlabel(r"$\theta/\pi$"); ax2.set_ylabel(r"largest $p$ passing the test")
ax2.set_title("Relaxation thresholds versus the exact boundary")
ax2.grid(alpha=0.3); ax2.legend()

ax3 = fig.add_subplot(2, 2, 3)
ax3.plot(ks, werner_sym, "s--", color="tab:blue", lw=2, label="symmetric extension")
ax3.plot(ks, werner_dps, "^:", color="tab:red", lw=2, label="PPT symmetric extension")
kfine = np.linspace(1, 4, 100)
ax3.plot(kfine, (kfine + 2) / (3 * kfine), "-", color="tab:gray", label=r"$(k+2)/(3k)$")
ax3.axhline(1 / 3, color="k", lw=2, label=r"separability limit $1/3$")
ax3.set_xticks(ks)
ax3.set_xlabel(r"number of copies $k$"); ax3.set_ylabel(r"largest $p$ passing the test")
ax3.set_title("Werner-type state: hierarchy converges to 1/3")
ax3.grid(alpha=0.3); ax3.legend()

ax4 = fig.add_subplot(2, 2, 4, projection="3d")
xs, ys, zs = [], [], []
for ik, k in enumerate(kk):
    for it in range(len(thetas)):
        xs.append(thetas[it] / np.pi); ys.append(k); zs.append(sym_th[k][it] - ppt_th[it])
ax4.bar3d(np.array(xs) - 0.004, np.array(ys) - 0.15, np.zeros(len(zs)), 0.008, 0.3, zs,
          color=plt.cm.viridis(np.array(zs) / max(zs)), alpha=0.9)
ax4.set_xlabel(r"$\theta/\pi$"); ax4.set_ylabel(r"$k$"); ax4.set_zlabel("gap to exact limit")
ax4.set_yticks(kk)
ax4.set_title("Gap of the symmetric-extension test")
ax4.view_init(elev=25, azim=-120)

plt.tight_layout()
plt.savefig("separability_results.png", dpi=150)
plt.show()

Code Walkthrough

1. States and the PPT criterion

phi(theta) builds $|\Phi_\theta\rangle=\cos\theta|00\rangle+\sin\theta|11\rangle$ as a length-4 vector, and rho_state(p, theta) mixes its projector with the maximally mixed state $I/4$.

partial_transpose_B reshapes the $4\times4$ matrix into a tensor with indices $(a,b,a’,b’)$ and swaps the two $B$ indices. That is exactly the partial transpose on $B$:

ppt_min_eig returns the smallest eigenvalue, and ppt_threshold implements the closed form $p^\ast(\theta)=1/(1+2\sin2\theta)$ derived above. This is the baseline against which both hierarchy tests are compared.

2. Symmetric extension as an SDP

The extension lives on $n=k+1$ qubits: qubit 0 is $A$ and qubits $1,\dots,k$ are the copies of $B$.

swap_matrix(n, i, j) returns the permutation matrix that exchanges qubits $i$ and $j$. Invariance under all permutations of $B$ is enforced by requiring invariance under the adjacent swaps $(j,j+1)$, which generate the full symmetric group.
trace_out_last removes the last qubit using $\operatorname{Tr}_{\text{last}}X=\sum_i (I\otimes\langle i|),X,(I\otimes|i\rangle)$. Applying it repeatedly produces the marginal on $AB_1$.
transpose_qubit implements the partial transpose of one qubit as a sum of four terms, using
$$
X^{T_q}=\sum_{i,j}E_{ij},X,E_{ij},\qquad E_{ij}=|i\rangle\langle j| ,
$$
so it only uses affine matrix products. This works in every CVXPY version and does not rely on a built-in partial-transpose function.
max_visibility assembles the SDP. The variable is a real symmetric matrix. Since $\rho(p,\theta)$ is real, this loses no generality: if a complex Hermitian extension exists, its real part is again a valid extension. This halves the number of variables and keeps the solver fast.
The constraint marg == p * target_pure + (1 - p) * I4 is linear in $(\sigma,p)$, so maximizing $p$ is a single SDP and needs no bisection.
With ppt=True, the code adds $\sigma^{T_{B_{k-m+1}\cdots B_k}}\succeq0$ for $m=1,\dots,k$, which turns the test into the PPT symmetric extension hierarchy.
solve uses the Clarabel interior-point solver and falls back to SCS if Clarabel is unavailable.

3. Speed

The SDP matrix has size $2^{k+1}$, so $k\le4$ means at most $32\times32$. The real-symmetric variable, the linear objective (no bisection) and the interior-point solver keep the whole script, including all grids and the figure, to a few seconds. No further acceleration is needed.

4. Example A: the Werner-type state

At $\theta=\pi/4$ the script first prints the minimum eigenvalue of $\rho^{T_B}$ for four values of $p$. It then computes the extension thresholds for $k=1,\dots,4$ and compares them with $(k+2)/(3k)$.

5. Example B: sweeping the entanglement angle

Eight values of $\theta\in[\pi/16,\pi/4]$ are tested at $k=2,3$, both with and without the PPT constraints, and compared with the exact boundary $p^\ast(\theta)$.

6. Visualization

All four panels are drawn in one figure and saved as separability_results.png.

Execution Results

=== Example A: Werner-type state, rho(p) = p|Phi+>＜Phi+| + (1-p) I/4 ===
p = 0.2000  min eig of rho^TB = +0.10000  -> PPT (separable)
p = 0.3333  min eig of rho^TB = +0.00000  -> PPT (separable)
p = 0.5000  min eig of rho^TB = -0.12500  -> NPT (entangled)
p = 0.8000  min eig of rho^TB = -0.35000  -> NPT (entangled)

Elapsed: 4.2 s

 k | sym. extension p_max | (k+2)/(3k) | PPT sym. ext. p_max
 1 |              1.00000 |    1.00000 |             0.33333
 2 |              0.66667 |    0.66667 |             0.33333
 3 |              0.55556 |    0.55556 |             0.33333
 4 |              0.50000 |    0.50000 |             0.33333

PPT threshold (exact separability limit): 0.33333

=== Example B: rho(p, theta) = p|Phi(theta)>＜Phi(theta)| + (1-p) I/4 ===
 theta/pi | PPT limit | sym k=2 | sym k=3 | PPT-sym k=2 | PPT-sym k=3
   0.0625 |   0.56645 | 0.90423 | 0.84031 |     0.56645 |     0.56645
   0.0893 |   0.48448 | 0.84405 | 0.75602 |     0.48448 |     0.48448
   0.1161 |   0.42869 | 0.78960 | 0.68803 |     0.42869 |     0.42869
   0.1429 |   0.39007 | 0.74457 | 0.63663 |     0.39007 |     0.39007
   0.1696 |   0.36358 | 0.70998 | 0.59954 |     0.36358 |     0.36358
   0.1964 |   0.34629 | 0.68572 | 0.57460 |     0.34629 |     0.34629
   0.2232 |   0.33650 | 0.67140 | 0.56024 |     0.33650 |     0.33650
   0.2500 |   0.33333 | 0.66667 | 0.55556 |     0.33333 |     0.33333

Reading the Results

Panel 1 (3D): the PPT landscape

The surface shows $\lambda_{\min}(\rho^{T_B})$ over $(\theta,p)$. The gray plane is zero, and the black curve where the surface crosses it is the separability boundary $p^\ast(\theta)$. Below the curve the partial transpose is positive and the state is separable; above it the eigenvalue goes negative and the state is entangled. The surface dips deepest near $\theta=\pi/4$ and $p=1$, the maximally entangled corner, where $\lambda_{\min}=-1/2$.

Panel 2: thresholds against the exact boundary

The black curve is the exact boundary. The plain symmetric extension test ($k=2,3$) sits well above it. States between the black curve and the dashed curves are entangled but still pass the test, so the test fails to detect them. Raising $k$ from 2 to 3 pulls the dashed curve down toward the boundary, which is the hierarchy at work. The PPT symmetric extension curves lie exactly on the black curve, so for two qubits that stronger test is already tight at $k=2$.

Panel 3: the Werner-type state as $k$ grows

For $\theta=\pi/4$, the plain symmetric extension thresholds agree with the formula

$$
p_k=\frac{k+2}{3k}\ :\qquad 1,\ \frac23,\ \frac59,\ \frac12\quad (k=1,2,3,4),
$$

and tend to the separability limit $\lim_{k\to\infty}p_k=\tfrac13$ from above. The convergence is slow: $k=4$ still accepts every state with $p\le 1/2$. The PPT-augmented version stays flat at $1/3$ for every $k$, since adding the PPT conditions already reproduces the exact answer for two qubits.

Panel 4 (3D): the gap

Each bar shows $p_k(\theta)-p^\ast(\theta)$ for the plain symmetric extension test. Bars at $k=3$ are consistently lower than bars at $k=2$, so the extra copy always removes part of the gap. The gap is smallest at small $\theta$ (weakly entangled states) and grows toward $\theta=\pi/4$ for both levels.

Takeaways

For two qubits, PPT is an exact separability test, and the PPT symmetric extension hierarchy reproduces it already at the first levels.
The plain symmetric extension hierarchy converges only gradually, with the Werner-type thresholds following $(k+2)/(3k)\to1/3$.
The same SDP formulation carries over unchanged to higher dimensions, where PPT is no longer sufficient and the hierarchy becomes the practical tool for certifying entanglement.

Quantum Optimal Decoding

October 8, 2026

Computing the Best Recovery Map After Noise with Python

When a quantum state passes through a noisy channel, the information is not gone, but it is smeared across a larger Hilbert space. The question is how to undo the damage as faithfully as physics allows. The map that does this best is called the optimal recovery map (or optimal decoder). In this article we compute it for a concrete case: a logical qubit stored in four physical qubits that suffer amplitude damping. We solve for the best recovery by an iteration that is equivalent to a semidefinite program (SDP), compare it with two standard decoders, and visualize the result in 2D and 3D.

1. The problem

We store one logical qubit in four physical qubits with the Leung–Nielsen–Chuang–Yamamoto code:

$$
|0_L\rangle=\frac{|0000\rangle+|1111\rangle}{\sqrt2},\qquad
|1_L\rangle=\frac{|0011\rangle+|1100\rangle}{\sqrt2}.
$$

The encoding is the isometry $V=|0_L\rangle\langle 0|+|1_L\rangle\langle 1|:\mathbb C^2\to\mathbb C^{16}$. Each physical qubit independently undergoes amplitude damping with parameter $\gamma$:

$$
K_0=\begin{pmatrix}1&0\0&\sqrt{1-\gamma}\end{pmatrix},\qquad
K_1=\begin{pmatrix}0&\sqrt{\gamma}\0&0\end{pmatrix},
$$

so the four-qubit noise channel has $16$ Kraus operators $K_{k}=K_{k_1}\otimes K_{k_2}\otimes K_{k_3}\otimes K_{k_4}$. Combining encoding and noise gives the effective noisy encoding with Kraus operators

$$
A_k=K_kV\in\mathbb C^{16\times 2}.
$$

Goal: find a recovery channel $\mathcal R(X)=\sum_j R_jXR_j^\dagger$ with $R_j\in\mathbb C^{2\times16}$ and $\sum_jR_j^\dagger R_j=I_{16}$ that maximizes the entanglement fidelity

$$
F_e(\mathcal R)=\frac1{d^2}\sum_{j,k}\Big|\mathrm{tr}\big(R_jA_k\big)\Big|^2,\qquad d=2 .
$$

The average fidelity over all pure input states follows from $F_e$:

$$
\bar F=\frac{dF_e+1}{d+1}.
$$

For an uncoded qubit, the best we can do is the identity map, and the entanglement fidelity is

$$
F_e^{\rm raw}=\frac{\big(1+\sqrt{1-\gamma}\big)^2}{4}\approx 1-\frac{\gamma}{2}.
$$

2. The optimization problem is an SDP

Write the Choi matrix of the recovery as $J=\sum_j|R_j\rangle!\rangle\langle!\langle R_j|$, where $|X\rangle!\rangle=\sum_{a,b}X_{ab},|b\rangle|a\rangle$ is the column-stacking vectorization. Define

$$
\Omega=\sum_k|\bar A_k\rangle!\rangle\langle!\langle\bar A_k|,\qquad
|\bar A_k\rangle!\rangle=\sum_{a,b}\overline{(A_k)_{ba}},|b\rangle|a\rangle .
$$

Then $\mathrm{tr}(R_jA_k)=\langle!\langle\bar A_k|R_j\rangle!\rangle$, and the fidelity becomes linear in $J$:

Here is the partial trace over the two-dimensional output factor, and the constraint is exactly trace preservation, $\sum_jR_j^\dagger R_j=I_{16}$.

3. A solver-free algorithm: polar-decomposition iteration

Stack all Kraus operators into one tall matrix $W=[R_1;R_2;\dots;R_m]\in\mathbb C^{2m\times16}$. Trace preservation says $W^\dagger W=I_{16}$, so $W$ is an isometry. The fidelity is a convex quadratic function of $W$:

$$
T_{jk}=\mathrm{tr}(R_jA_k),\qquad
f(W)=\frac1{d^2}\sum_{jk}|T_{jk}|^2 .
$$

Its gradient (with respect to $\bar W$, up to the factor $1/d^2$) has blocks

$$
G_j=\sum_kT_{jk},A_k^\dagger .
$$

The update maximizes the linearization of $f$ over all isometries. With the thin singular value decomposition $G=U\Sigma V_h$, the maximizer of $\mathrm{Re},\mathrm{tr}(G^\dagger W)$ over isometries is

$$
W_{\rm new}=U,V_h .
$$

Because $f$ is convex, $f(W_{\rm new})\ge f(W)+\frac{2}{d^2}\mathrm{Re},\mathrm{tr}\big[G^\dagger(W_{\rm new}-W)\big]\ge f(W)$, so the fidelity increases monotonically. Choosing $m=32$ Kraus operators is fully general, because the Choi matrix is $32\times32$.

4. Two reference decoders

Naive decoding. Project onto the code space and undo the encoding with $V^\dagger$. Anything outside the code space is replaced by the maximally mixed qubit so that the map is trace preserving.

Petz recovery map. With the reference state $\sigma=VV^\dagger/2$ and the noise channel $\mathcal N$,

$$
\mathcal R_{P}(X)=\sigma^{1/2},\mathcal N^\dagger!\Big(\mathcal N(\sigma)^{-1/2},X,\mathcal N(\sigma)^{-1/2}\Big),\sigma^{1/2},
$$

followed by the decoding $V^\dagger$. In Kraus form, $R_k=V^\dagger\sigma^{1/2}K_k^\dagger,\mathcal N(\sigma)^{-1/2}$. The Petz map is a famous near-optimal recovery map, which makes it a strong benchmark.

5. Complete source code

import time
from functools import reduce
from itertools import product

import numpy as np
import matplotlib.pyplot as plt

np.set_printoptions(precision=6, suppress=True)
plt.style.use("dark_background")

t_start = time.time()
rng = np.random.default_rng(2026)

# ------------------------------------------------------------
# 1. Code space and noise model
# ------------------------------------------------------------
N_QUBITS = 4
DIM = 2 ** N_QUBITS      # physical Hilbert-space dimension (16)
D_LOG = 2                # logical qubit dimension


def ket(bits):
    v = np.zeros(DIM, dtype=complex)
    v[int(bits, 2)] = 1.0
    return v


zero_L = (ket("0000") + ket("1111")) / np.sqrt(2)
one_L = (ket("0011") + ket("1100")) / np.sqrt(2)
V = np.stack([zero_L, one_L], axis=1)          # encoding isometry, shape (16, 2)


def amp_damp_kraus(g):
    k0 = np.array([[1.0, 0.0], [0.0, np.sqrt(1.0 - g)]], dtype=complex)
    k1 = np.array([[0.0, np.sqrt(g)], [0.0, 0.0]], dtype=complex)
    return [k0, k1]


def multi_qubit_kraus(g, n):
    single = amp_damp_kraus(g)
    ops = [reduce(np.kron, [single[c] for c in combo]) for combo in product(range(2), repeat=n)]
    return np.array(ops)                        # shape (2^n, 2^n, 2^n)


# ------------------------------------------------------------
# 2. Entanglement fidelity (linear in the Choi matrix of the recovery)
# ------------------------------------------------------------
def entanglement_fidelity(R, A):
    # R: (m, 2, n_out) recovery Kraus operators, A: (K, n_out, 2) effective noisy-encoding Kraus operators
    T = np.einsum("jab,kba->jk", R, A)
    return float(np.sum(np.abs(T) ** 2) / D_LOG ** 2)


# ------------------------------------------------------------
# 3. Optimal recovery by polar-decomposition (Reimpell-Werner type) iteration
# ------------------------------------------------------------
def random_isometry(rows, cols, generator):
    X = generator.normal(size=(rows, cols)) + 1j * generator.normal(size=(rows, cols))
    Q, _ = np.linalg.qr(X)
    return Q


def optimal_recovery(A, m=32, n_iter=400, n_starts=3, generator=rng, tol=1e-14):
    n_out = A.shape[1]
    A_dag = np.conj(np.transpose(A, (0, 2, 1)))
    best_F, best_R, histories = -1.0, None, []
    for _ in range(n_starts):
        W = random_isometry(m * D_LOG, n_out, generator)
        hist = []
        for _ in range(n_iter):
            R = W.reshape(m, D_LOG, n_out)
            T = np.einsum("jab,kba->jk", R, A)
            hist.append(float(np.sum(np.abs(T) ** 2) / D_LOG ** 2))
            G = np.einsum("jk,kab->jab", T, A_dag).reshape(m * D_LOG, n_out)
            U, _, Vh = np.linalg.svd(G, full_matrices=False)
            W = U @ Vh
            if len(hist) > 1 and abs(hist[-1] - hist[-2]) < tol:
                break
        R = W.reshape(m, D_LOG, n_out)
        F = entanglement_fidelity(R, A)
        hist.append(F)
        histories.append(hist)
        if F > best_F:
            best_F, best_R = F, R
    return best_F, best_R, histories


# ------------------------------------------------------------
# 4. Baseline recoveries
# ------------------------------------------------------------
def complete_tp(R):
    n_out = R.shape[2]
    S = np.einsum("jab,jac->bc", np.conj(R), R)
    D = np.eye(n_out) - S
    D = (D + D.conj().T) / 2
    w, U = np.linalg.eigh(D)
    extra = []
    for lam, phi in zip(w, U.T):
        if lam > 1e-12:
            for a in range(D_LOG):
                M = np.zeros((D_LOG, n_out), dtype=complex)
                M[a, :] = np.sqrt(lam / 2.0) * np.conj(phi)
                extra.append(M)
    if extra:
        return np.concatenate([R, np.array(extra)], axis=0)
    return R


def naive_recovery():
    R0 = np.conj(V.T).reshape(1, D_LOG, DIM)    # project onto the code space and decode
    return complete_tp(R0)


def psd_power(M, power):
    w, U = np.linalg.eigh((M + M.conj().T) / 2)
    w = np.clip(w, 0.0, None)
    out = np.zeros_like(w)
    mask = w > 1e-12
    out[mask] = w[mask] ** power
    return (U * out) @ U.conj().T


def petz_recovery(K):
    sigma = V @ V.conj().T / D_LOG                       # encoded maximally mixed state
    N_sigma = np.einsum("kab,bc,kdc->ad", K, sigma, np.conj(K))
    s_half = psd_power(sigma, 0.5)
    n_inv_half = psd_power(N_sigma, -0.5)
    R = np.array([np.conj(V.T) @ s_half @ np.conj(k.T) @ n_inv_half for k in K])
    return complete_tp(R)


# ------------------------------------------------------------
# 5. Bloch-sphere representation of a qubit channel
# ------------------------------------------------------------
PAULI = [
    np.array([[0, 1], [1, 0]], dtype=complex),
    np.array([[0, -1j], [1j, 0]], dtype=complex),
    np.array([[1, 0], [0, -1]], dtype=complex),
]


def bloch_affine(E):
    S = np.einsum("mab,mcb->ac", E, np.conj(E))
    t = np.array([np.real(np.trace(P @ S)) / 2.0 for P in PAULI])
    M = np.zeros((3, 3))
    for j in range(3):
        Xj = np.einsum("mab,bc,mdc->ad", E, PAULI[j], np.conj(E))
        for i in range(3):
            M[i, j] = np.real(np.trace(PAULI[i] @ Xj)) / 2.0
    return M, t


def effective_kraus(R, A):
    return np.einsum("jab,kbc->jkac", R, A).reshape(-1, D_LOG, D_LOG)


# ------------------------------------------------------------
# 6. Sweep over the damping parameter
# ------------------------------------------------------------
gammas = np.geomspace(0.003, 0.30, 28)
F_raw, F_naive, F_petz, F_opt = [], [], [], []
affine_raw, affine_opt = [], []
tp_error = 0.0

R_naive = naive_recovery()

for g in gammas:
    K_phys = multi_qubit_kraus(g, N_QUBITS)
    A = np.matmul(K_phys, V)                              # (16, 16, 2)

    K_single = np.array(amp_damp_kraus(g))                # uncoded qubit
    R_id = np.eye(D_LOG, dtype=complex).reshape(1, D_LOG, D_LOG)
    F_raw.append(entanglement_fidelity(R_id, K_single))
    affine_raw.append(bloch_affine(K_single))

    F_naive.append(entanglement_fidelity(R_naive, A))

    R_petz = petz_recovery(K_phys)
    F_petz.append(entanglement_fidelity(R_petz, A))

    F_best, R_best, _ = optimal_recovery(A)
    F_opt.append(F_best)
    affine_opt.append(bloch_affine(effective_kraus(R_best, A)))
    S = np.einsum("jab,jac->bc", np.conj(R_best), R_best)
    tp_error = max(tp_error, float(np.max(np.abs(S - np.eye(DIM)))))

F_raw, F_naive, F_petz, F_opt = map(np.array, (F_raw, F_naive, F_petz, F_opt))

# ------------------------------------------------------------
# 7. Console report
# ------------------------------------------------------------
print("Isometry check  max|V^dag V - I| =", float(np.max(np.abs(V.conj().T @ V - np.eye(D_LOG)))))
K_chk = multi_qubit_kraus(0.1, N_QUBITS)
print("Noise TP check  max|sum K^dag K - I| =",
      float(np.max(np.abs(np.einsum("kab,kac->bc", np.conj(K_chk), K_chk) - np.eye(DIM)))))
print("Recovery TP check (worst over sweep)  max|sum R^dag R - I| =", tp_error)
print()
print("Entanglement fidelity F_e")
print(f"{'gamma':>8} {'no code':>10} {'naive':>10} {'Petz':>10} {'optimal':>10}")
for i in range(0, len(gammas), 3):
    print(f"{gammas[i]:8.4f} {F_raw[i]:10.6f} {F_naive[i]:10.6f} {F_petz[i]:10.6f} {F_opt[i]:10.6f}")
print()

n_fit = 8
for name, arr in [("no code", F_raw), ("naive", F_naive), ("Petz", F_petz), ("optimal", F_opt)]:
    slope = np.polyfit(np.log(gammas[:n_fit]), np.log(np.maximum(1.0 - arr[:n_fit], 1e-16)), 1)[0]
    print(f"Scaling exponent of 1 - F_e at small gamma ({name:>8}): {slope:.3f}")
print()

g_demo = 0.10
A_demo = np.matmul(multi_qubit_kraus(g_demo, N_QUBITS), V)
F_demo, R_demo, hists = optimal_recovery(A_demo, n_starts=3)
print(f"gamma = {g_demo}: optimal F_e = {F_demo:.8f}, average state fidelity = {(D_LOG * F_demo + 1) / (D_LOG + 1):.8f}")
print("Final fidelity of each random start:", [round(h[-1], 8) for h in hists])
print("Iterations used by each start:", [len(h) - 1 for h in hists])
print(f"Elapsed time: {time.time() - t_start:.1f} s")

# ------------------------------------------------------------
# 8. Visualisation (single figure)
# ------------------------------------------------------------
fig = plt.figure(figsize=(22, 13))
gs = fig.add_gridspec(2, 3)

ax1 = fig.add_subplot(gs[0, 0])
ax1.plot(gammas, F_raw, "o-", color="#ff6b6b", label="no code")
ax1.plot(gammas, F_naive, "s-", color="#feca57", label="4-qubit code, naive decoding")
ax1.plot(gammas, F_petz, "^-", color="#48dbfb", label="4-qubit code, Petz recovery")
ax1.plot(gammas, F_opt, "d-", color="#1dd1a1", label="4-qubit code, optimal recovery")
ax1.set_xlabel("damping parameter $\\gamma$")
ax1.set_ylabel("entanglement fidelity $F_e$")
ax1.set_title("Fidelity versus noise strength")
ax1.grid(alpha=0.25)
ax1.legend(loc="lower left")

ax2 = fig.add_subplot(gs[0, 1])
ax2.loglog(gammas, np.maximum(1 - F_raw, 1e-16), "o-", color="#ff6b6b", label="no code")
ax2.loglog(gammas, np.maximum(1 - F_naive, 1e-16), "s-", color="#feca57", label="naive")
ax2.loglog(gammas, np.maximum(1 - F_petz, 1e-16), "^-", color="#48dbfb", label="Petz")
ax2.loglog(gammas, np.maximum(1 - F_opt, 1e-16), "d-", color="#1dd1a1", label="optimal")
ax2.loglog(gammas, gammas * (1 - F_raw[0]) / gammas[0], "--", color="white", alpha=0.5, label="$\\propto \\gamma$")
ax2.loglog(gammas, gammas ** 2 * (1 - F_opt[0]) / gammas[0] ** 2, ":", color="white", alpha=0.7, label="$\\propto \\gamma^2$")
ax2.set_xlabel("damping parameter $\\gamma$")
ax2.set_ylabel("infidelity $1-F_e$")
ax2.set_title("Error scaling (log-log)")
ax2.grid(alpha=0.25, which="both")
ax2.legend(loc="lower right")

ax3 = fig.add_subplot(gs[0, 2])
F_star = max(h[-1] for h in hists)
for idx, h in enumerate(hists):
    gap = np.maximum(F_star - np.array(h[:-1]), 1e-16)
    ax3.semilogy(np.arange(len(gap)), gap, lw=2, label=f"random start {idx + 1}")
ax3.set_xlim(0, 120)
ax3.set_xlabel("iteration")
ax3.set_ylabel("$F_e^{\\star} - F_e^{(t)}$")
ax3.set_title(f"Convergence of the iteration ($\\gamma={g_demo}$)")
ax3.grid(alpha=0.25, which="both")
ax3.legend()

i_show = len(gammas) - 1
u = np.linspace(0, 2 * np.pi, 60)
v = np.linspace(0, np.pi, 40)
sx = np.outer(np.cos(u), np.sin(v))
sy = np.outer(np.sin(u), np.sin(v))
sz = np.outer(np.ones_like(u), np.cos(v))
sphere = np.stack([sx, sy, sz])


def draw_bloch_image(ax, M, t, color, title):
    ax.plot_wireframe(sx, sy, sz, color="white", alpha=0.12, rstride=4, cstride=4)
    img = np.einsum("ij,jxy->ixy", M, sphere) + t[:, None, None]
    ax.plot_surface(img[0], img[1], img[2], color=color, alpha=0.65, linewidth=0)
    ax.set_xlim(-1, 1)
    ax.set_ylim(-1, 1)
    ax.set_zlim(-1, 1)
    ax.set_box_aspect((1, 1, 1))
    ax.set_xlabel("x")
    ax.set_ylabel("y")
    ax.set_zlabel("z")
    ax.set_title(title)
    ax.view_init(elev=20, azim=35)


ax4 = fig.add_subplot(gs[1, 0], projection="3d")
draw_bloch_image(ax4, affine_raw[i_show][0], affine_raw[i_show][1], "#ff6b6b",
                 f"Bloch sphere image: no code ($\\gamma={gammas[i_show]:.2f}$)")

ax5 = fig.add_subplot(gs[1, 1], projection="3d")
draw_bloch_image(ax5, affine_opt[i_show][0], affine_opt[i_show][1], "#1dd1a1",
                 f"Bloch sphere image: optimal recovery ($\\gamma={gammas[i_show]:.2f}$)")

thetas = np.linspace(0, np.pi, 41)
Z_raw = np.zeros((len(thetas), len(gammas)))
Z_opt = np.zeros((len(thetas), len(gammas)))
for ig in range(len(gammas)):
    for it, th in enumerate(thetas):
        n_in = np.array([np.sin(th), 0.0, np.cos(th)])
        for Z, (M, t) in ((Z_raw, affine_raw[ig]), (Z_opt, affine_opt[ig])):
            n_out = M @ n_in + t
            Z[it, ig] = 0.5 * (1.0 + float(n_in @ n_out))

GG, TT = np.meshgrid(gammas, thetas)
ax6 = fig.add_subplot(gs[1, 2], projection="3d")
ax6.plot_surface(GG, TT, Z_raw, cmap="autumn", alpha=0.75, linewidth=0)
ax6.plot_surface(GG, TT, Z_opt, cmap="viridis", alpha=0.95, linewidth=0)
ax6.set_xlabel("$\\gamma$")
ax6.set_ylabel("input polar angle $\\theta$")
ax6.set_zlabel("state fidelity")
ax6.set_title("State fidelity surface: no code (warm) vs optimal recovery (cool)")
ax6.view_init(elev=24, azim=-125)

plt.tight_layout()
plt.show()

6. Code walkthrough

Section 1: code space and noise

ket builds a computational basis vector of the 16-dimensional space from a bit string. Stacking zero_L and one_L as columns gives the isometry V of shape $16\times2$. amp_damp_kraus returns the two single-qubit Kraus operators, and multi_qubit_kraus forms all $2^4=16$ tensor products with reduce(np.kron, ...), which is the Kraus set of the four-qubit noise channel. Later, np.matmul(K_phys, V) broadcasts over the Kraus index and produces the effective operators $A_k=K_kV$ in a single call.

Section 2: the fidelity function

entanglement_fidelity is the formula $F_e=\frac1{d^2}\sum_{jk}|\mathrm{tr}(R_jA_k)|^2$. The einsum string "jab,kba->jk" evaluates for every pair $(j,k)$ at once. The same function scores every strategy, which keeps the comparison fair: the uncoded qubit uses $R=I$ with $A_k=K_k$, and the coded cases use $16\times2$ operators $A_k$.

Section 3: the optimizer

This is the heart of the article. random_isometry produces a random starting $W$ by taking the QR factorization of a complex Gaussian matrix, which gives orthonormal columns. In each iteration of optimal_recovery:

W.reshape(m, 2, n_out) unstacks $W$ into the $m$ Kraus operators $R_j$.
T is the matrix $T_{jk}=\mathrm{tr}(R_jA_k)$, and its squared magnitude sum is the current fidelity, which is stored in hist.
G = einsum("jk,kab->jab", T, A_dag) computes $G_j=\sum_kT_{jk}A_k^\dagger$, and reshape stacks the blocks back into a $64\times16$ matrix.
The thin SVD G = U Σ Vh gives the closest isometry W = U @ Vh, which is the polar factor and the exact maximizer of the linearized objective.

The loop stops early if the fidelity changes by less than tol. Three random starts are run, and the best final fidelity is kept. If all starts agree, that is strong numerical evidence that the iteration has found the global optimum of the SDP.

Speed. Everything is vectorized with einsum, and the only decomposition per iteration is the SVD of a $64\times16$ matrix. A complete $\gamma$ sweep with 28 values and 3 random starts each finishes within a few seconds on a standard CPU, so no further acceleration is needed.

Section 4: baselines

complete_tp takes any set of Kraus operators whose $\sum R^\dagger R$ is dominated by the identity and appends extra operators that send the missing part of the space to the maximally mixed qubit. It does this with the eigendecomposition of $D=I-\sum R^\dagger R$: for every eigenvector $|\phi\rangle$ with eigenvalue $\lambda$, it adds $\sqrt{\lambda/2},|a\rangle\langle\phi|$ for $a=0,1$, which restores trace preservation exactly. naive_recovery is just $V^\dagger$ plus this completion. psd_power raises a positive semidefinite matrix to a fractional power on its support, which is used for $\sigma^{1/2}$ and the pseudo-inverse square root $\mathcal N(\sigma)^{-1/2}$. petz_recovery builds $R_k=V^\dagger\sigma^{1/2}K_k^\dagger\mathcal N(\sigma)^{-1/2}$ for every Kraus operator of the noise.

Section 5: Bloch-sphere tools

Any qubit channel acts affinely on Bloch vectors, $\vec n\mapsto M\vec n+\vec t$. bloch_affine extracts $M$ and $\vec t$ from a Kraus set using the Pauli matrices: $\vec t_i=\tfrac12\mathrm{tr}(\sigma_i\sum_mE_mE_m^\dagger)$ and $M_{ij}=\tfrac12\mathrm{tr}(\sigma_i\sum_mE_m\sigma_jE_m^\dagger)$. effective_kraus composes recovery and noisy encoding, $E_{jk}=R_jA_k$, so that the whole “encode, damp, recover” pipeline becomes one qubit-to-qubit channel that can be drawn.

Section 6: the $\gamma$ sweep

For each damping value the loop builds the noise Kraus operators, scores the uncoded qubit, the naive decoder, the Petz map and the optimal map, and stores the Bloch-affine parameters for the 3D plots. tp_error records the worst violation of $\sum_jR_j^\dagger R_j=I$ over the sweep, which certifies that the computed optimal map is a legitimate channel.

Section 7: console report

The report prints sanity checks (isometry, trace preservation of noise and recovery), a fidelity table, the fitted scaling exponents of the infidelity, and a detailed run at $\gamma=0.1$ with the final fidelity of each random start. The exponents come from a straight-line fit of $\log(1-F_e)$ against $\log\gamma$ over the eight smallest damping values.

Section 8: the figure

All six panels live in a single figure built with add_gridspec(2, 3). The first row has three 2D panels: fidelity versus $\gamma$, the log-log infidelity, and the convergence of the iteration. The second row has three 3D panels: the image of the Bloch sphere without coding, the image under the optimal recovery, and a fidelity surface over $(\gamma,\theta)$. For a pure input with Bloch vector $\vec n_{\rm in}=(\sin\theta,0,\cos\theta)$ and output vector $\vec n_{\rm out}=M\vec n_{\rm in}+\vec t$, the state fidelity is

$$
F(\theta,\gamma)=\frac{1+\vec n_{\rm in}\cdot\vec n_{\rm out}}{2},
$$

which is computed directly without building any density matrix.

7. Execution results

Isometry check  max|V^dag V - I| = 2.220446049250313e-16
Noise TP check  max|sum K^dag K - I| = 4.440892098500626e-16
Recovery TP check (worst over sweep)  max|sum R^dag R - I| = 2.4424906541753444e-15

Entanglement fidelity F_e
   gamma    no code      naive       Petz    optimal
  0.0030   0.998499   0.995509   0.999984   0.999989
  0.0050   0.997496   0.992519   0.999956   0.999969
  0.0083   0.995822   0.987548   0.999878   0.999913
  0.0139   0.993025   0.979306   0.999660   0.999758
  0.0232   0.988352   0.965695   0.999053   0.999326
  0.0387   0.980531   0.943367   0.997359   0.998124
  0.0646   0.967414   0.907160   0.992639   0.994781
  0.1078   0.945324   0.849589   0.979520   0.985493
  0.1798   0.907851   0.761122   0.943525   0.959785
  0.3000   0.843330   0.633250   0.849138   0.889663

Scaling exponent of 1 - F_e at small gamma ( no code): 1.001
Scaling exponent of 1 - F_e at small gamma (   naive): 0.996
Scaling exponent of 1 - F_e at small gamma (    Petz): 2.001
Scaling exponent of 1 - F_e at small gamma ( optimal): 2.000

gamma = 0.1: optimal F_e = 0.98751670, average state fidelity = 0.99167780
Final fidelity of each random start: [0.9875167, 0.9875167, 0.9875167]
Iterations used by each start: [400, 400, 400]
Elapsed time: 10.7 s

8. How to read the results

Panel 1: fidelity versus noise strength

Four curves start near $F_e=1$ at small $\gamma$ and fall as the noise grows. The optimal recovery sits on top, the Petz map follows closely behind, and the uncoded qubit comes next. The naive decoder is the lowest curve, and it is even worse than not encoding at all. The reason is that the four-qubit code is designed so that a single damping jump can be detected and corrected, but the “no-jump” evolution, in which nothing visibly happens, still distorts the code words at first order in $\gamma$. A decoder that merely projects and decodes ignores that distortion. The optimal map compensates for it.

Panel 2: error scaling on a log-log plot

This panel contains the main message. On a log-log plot, a power law $1-F_e\propto\gamma^p$ is a straight line of slope $p$. The uncoded qubit and the naive decoder both have slope $p\approx1$, matching $1-F_e^{\rm raw}\approx\gamma/2$. The Petz map and the optimal map have slope $p\approx2$: the code together with a good recovery cancels all first-order errors, and only second-order processes (two damping events, or a jump combined with a distortion) survive. The optimal curve lies consistently below the Petz curve by a roughly constant factor, so Petz is excellent but not exactly optimal. At small noise, the benefit of coding is therefore not a few percent but orders of magnitude.

Panel 3: convergence of the iteration

The vertical axis shows the gap between each iterate’s fidelity and the converged value. Starting from random isometries, the gap collapses by many orders of magnitude within the first handful of iterations, then flattens at a level of about $10^{-9}$ or below. The three starts reach the same fidelity to roughly eight significant digits, which is what we expect when the iteration lands on the global optimum of the underlying SDP. As a cross-check, a generic interior-point SDP solver gives the same optimum to about eight digits. The small plateau differences between starts are caused by directions along which the objective is almost flat.

Panels 4 and 5: the Bloch sphere before and after

The faint wireframe is the original Bloch sphere. A channel maps it to an ellipsoid, and a perfect channel would leave it untouched. Without coding, amplitude damping squeezes the sphere toward the north pole $|0\rangle$: the equatorial radius shrinks to $\sqrt{1-\gamma}$, the polar axis shrinks to $1-\gamma$, and the whole shape shifts upward by $\gamma$. This is the physical signature of energy relaxation, which loses coherence and also biases the state toward the ground state. With the code and the optimal recovery at the same $\gamma=0.30$, the image is almost a centered sphere. A small residual shrinkage remains because even the optimal recovery cannot fully undo second-order processes.

Panel 6: the fidelity surface

The surface shows the fidelity of a single pure input as a function of the damping strength and the input polar angle $\theta$, where $\theta=0$ is $|0\rangle$ and $\theta=\pi$ is $|1\rangle$. The warm surface (no code) is perfectly flat at the value 1 along $\theta=0$, because $|0\rangle$ is the fixed point of amplitude damping. It collapses at $\theta=\pi$, because $|1\rangle$ is the state that decays. The cool surface (optimal recovery) is much more uniform in $\theta$. Because the recovery is optimized for the average entanglement fidelity, it gives up a small amount at the naturally protected state $|0\rangle$ in exchange for large gains at the vulnerable states. This uniformity is exactly what we want from a quantum memory, which must preserve every state equally well, not just the convenient ones.

9. Summary

The optimal recovery map is the solution of a semidefinite program: maximize a linear function of the Choi matrix subject to positivity and trace preservation.
Exploiting the structure of the entanglement fidelity, a simple polar-decomposition iteration solves the problem monotonically using only NumPy.
For the four-qubit amplitude-damping code, the optimal recovery has infidelity scaling as $\gamma^2$, the Petz map has the same scaling with a slightly larger prefactor, and naive project-and-decode decoding scales only as $\gamma$ and is worse than doing nothing.
The recovery map matters as much as the code itself: the same encoded state can be restored with first-order or second-order errors depending solely on how we decode.

The Diamond Norm Computed with a Semidefinite Program

October 7, 2026

Telling Two Quantum Channels Apart

How well can you tell two quantum channels apart? Suppose a black box applies either channel $\Phi_0$ or channel $\Phi_1$, each with probability $1/2$. You may feed it any input, including half of an entangled state, and measure the output however you like. The best possible success probability is

$$
P_{\text{succ}}^{\max}=\frac12+\frac14,\big|\Phi_0-\Phi_1\big|_{\diamond},
$$

where $|\cdot|_\diamond$ is the diamond norm. This post defines the norm, turns it into a semidefinite program (SDP), and solves it in Python for concrete channels.

1. The diamond norm

For a Hermiticity-preserving map $\Delta:\mathcal L(\mathcal H_{\text{in}})\to\mathcal L(\mathcal H_{\text{out}})$, the diamond norm is

where the maximum runs over all states on an input system $A$ plus a reference (ancilla) system $R$. The trace norm is $|X|_1=\mathrm{Tr}\sqrt{X^\dagger X}$.

The ancilla matters. A channel may look identical to another on every single input state and still be distinguishable once the input is entangled with a reference.

2. The Choi matrix

For a channel $\Phi$ with $d_{\text{in}}=d_{\text{out}}=d$, define

$$
J(\Phi)=\sum_{i,j=1}^{d}|i\rangle\langle j|\otimes\Phi\big(|i\rangle\langle j|\big).
$$

Then $J(\Phi)$ is positive semidefinite iff $\Phi$ is completely positive. For a difference of two channels, $J(\Delta)=J(\Phi_0)-J(\Phi_1)$ is a Hermitian matrix. It is generally indefinite.

3. The SDP (Watrous’ formulation)

The diamond norm of $\Delta$ is the optimal value of

The variables $Y_0,Y_1$ are Hermitian matrices of size $d_{\text{in}}d_{\text{out}}$, and $t_0,t_1$ are real scalars. The two scalar constraints are exactly the statements $|\mathrm{Tr}_{\text{out}}Y_k|_\infty\le t_k$, because the $Y_k$ are positive semidefinite. Everything is linear in the unknowns, so this is a standard SDP.

4. The example problems

Example 1: identity versus depolarizing. The depolarizing channel is $\mathcal D_p(\rho)=(1-p)\rho+p,\mathrm{Tr}(\rho),\mathbb 1/2$. Then $\mathrm{id}-\mathcal D_p=p,(\mathrm{id}-\mathcal D_1)$. The optimal input is the maximally entangled state $|\Phi^+\rangle$, and

$$
|\mathrm{id}-\mathcal D_p|_\diamond=p,\Big|,|\Phi^+\rangle\langle\Phi^+|-\tfrac{\mathbb 1}{4}\Big|_1=\frac{3p}{2}.
$$

Example 2: identity versus a $Z$-rotation. Let $\mathcal R_\theta(\rho)=U_\theta\rho U_\theta^\dagger$ with $U_\theta=\mathrm{diag}(e^{-i\theta/2},e^{i\theta/2})$. The difference of two unitary channels has the closed form

(the minimum is over density matrices). For the $Z$-rotation this evaluates to

$$
|\mathrm{id}-\mathcal R_\theta|_\diamond=2,\big|\sin(\theta/2)\big| .
$$

Example 3: depolarizing versus amplitude damping. The amplitude damping channel has Kraus operators

$$
K_0=\begin{pmatrix}1&0\0&\sqrt{1-\gamma}\end{pmatrix},\qquad
K_1=\begin{pmatrix}0&\sqrt{\gamma}\0&0\end{pmatrix}.
$$

No simple formula exists for $|\mathcal D_p-\mathcal A_\gamma|_\diamond$. We compute it on a grid of $(p,\gamma)$ and read off the optimal success probability $\tfrac12+\tfrac14|\mathcal D_p-\mathcal A_\gamma|_\diamond$.

Entanglement advantage. For comparison we also compute the best distinguishability achievable without an ancilla,

$$
\max_{|\psi\rangle}\big|\Phi_0(\psi)-\Phi_1(\psi)\big|_1,
$$

by sampling the Bloch sphere. The gap to the diamond norm is the benefit of entangled inputs.

5. Source code

import time
import warnings
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec

try:
    import cvxpy as cp
except ImportError:
    import subprocess, sys
    subprocess.check_call([sys.executable, "-m", "pip", "install", "-q", "cvxpy"])
    import cvxpy as cp

# ------------------------------------------------------------
# 1. Channels (all single-qubit: d_in = d_out = 2)
# ------------------------------------------------------------
D = 2
I2 = np.eye(2, dtype=complex)
Z = np.array([[1, 0], [0, -1]], dtype=complex)


def depolarizing(p):
    return lambda X: (1 - p) * X + p * np.trace(X) * I2 / 2


def dephasing(p):
    return lambda X: (1 - p) * X + p * Z @ X @ Z


def amplitude_damping(g):
    K0 = np.array([[1, 0], [0, np.sqrt(1 - g)]], dtype=complex)
    K1 = np.array([[0, np.sqrt(g)], [0, 0]], dtype=complex)
    return lambda X: K0 @ X @ K0.conj().T + K1 @ X @ K1.conj().T


def z_rotation(theta):
    U = np.diag([np.exp(-1j * theta / 2), np.exp(1j * theta / 2)])
    return lambda X: U @ X @ U.conj().T


def identity_channel():
    return lambda X: X


def choi(channel, d=D):
    """J = sum_ij |i><j| (input) kron Phi(|i><j|) (output)."""
    J = np.zeros((d * d, d * d), dtype=complex)
    for i in range(d):
        for j in range(d):
            E = np.zeros((d, d), dtype=complex)
            E[i, j] = 1.0
            J += np.kron(E, channel(E))
    return J


# ------------------------------------------------------------
# 2. Diamond norm via SDP (Watrous' dual formulation)
# ------------------------------------------------------------
def pick_solver():
    solvers = cp.installed_solvers()
    return cp.CLARABEL if "CLARABEL" in solvers else cp.SCS


SOLVER = pick_solver()


def diamond_norm(J, d_in=D, d_out=D):
    n = d_in * d_out
    Y0 = cp.Variable((n, n), hermitian=True)
    Y1 = cp.Variable((n, n), hermitian=True)
    t0 = cp.Variable()
    t1 = cp.Variable()

    # Partial trace over the output system, written as a sum of constant sandwiches
    A = [np.kron(np.eye(d_in), np.eye(d_out)[k].reshape(1, -1)) for k in range(d_out)]

    def ptrace_out(Y):
        return sum(Ak @ Y @ Ak.T for Ak in A)

    block = cp.bmat([[Y0, -J], [-J.conj().T, Y1]])
    constraints = [
        block >> 0,
        ptrace_out(Y0) << t0 * np.eye(d_in),
        ptrace_out(Y1) << t1 * np.eye(d_in),
    ]
    prob = cp.Problem(cp.Minimize((t0 + t1) / 2), constraints)
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        prob.solve(solver=SOLVER)
    return float(prob.value)


def diamond_distance(ch0, ch1):
    J = choi(ch0) - choi(ch1)
    return diamond_norm(J)


def helstrom_success(dn):
    return 0.5 + 0.25 * dn


# ------------------------------------------------------------
# 3. Trace norm without ancilla (single input state) for comparison
# ------------------------------------------------------------
def trace_norm(M):
    return float(np.sum(np.abs(np.linalg.eigvalsh((M + M.conj().T) / 2))))


def bloch_state(theta, phi):
    psi = np.array([np.cos(theta / 2), np.exp(1j * phi) * np.sin(theta / 2)])
    return np.outer(psi, psi.conj())


def best_single_input(ch0, ch1, n_theta=61, n_phi=121):
    best = 0.0
    for th in np.linspace(0, np.pi, n_theta):
        for ph in np.linspace(0, 2 * np.pi, n_phi):
            rho = bloch_state(th, ph)
            best = max(best, trace_norm(ch0(rho) - ch1(rho)))
    return best


# ------------------------------------------------------------
# 4. Experiments
# ------------------------------------------------------------
t_start = time.time()

# (A) Depolarizing(p) vs identity : analytic value 3p/2
ps = np.linspace(0, 1, 11)
dn_depol = np.array([diamond_distance(identity_channel(), depolarizing(p)) for p in ps])
an_depol = 1.5 * ps

# (B) Z-rotation(theta) vs identity : analytic value 2|sin(theta/2)|
thetas = np.linspace(0, 2 * np.pi, 25)
dn_rot = np.array([diamond_distance(identity_channel(), z_rotation(t)) for t in thetas])
an_rot = 2 * np.abs(np.sin(thetas / 2))

# (C) Depolarizing(p) vs amplitude damping(gamma) on a grid
n_grid = 13
p_grid = np.linspace(0, 1, n_grid)
g_grid = np.linspace(0, 1, n_grid)
surface = np.zeros((n_grid, n_grid))
for a, p in enumerate(p_grid):
    for b, g in enumerate(g_grid):
        surface[a, b] = diamond_distance(depolarizing(p), amplitude_damping(g))
success = helstrom_success(surface)

# (D) Single-input (no ancilla) vs diamond norm on the Bloch sphere for a fixed pair
p_fix, g_fix = 0.6, 0.3
ch_a, ch_b = depolarizing(p_fix), amplitude_damping(g_fix)
dn_fix = diamond_distance(ch_a, ch_b)
n_th, n_ph = 40, 80
TH, PH = np.meshgrid(np.linspace(0, np.pi, n_th), np.linspace(0, 2 * np.pi, n_ph), indexing="ij")
BLOCH = np.zeros_like(TH)
for i in range(n_th):
    for j in range(n_ph):
        rho = bloch_state(TH[i, j], PH[i, j])
        BLOCH[i, j] = trace_norm(ch_a(rho) - ch_b(rho))
best_fix = BLOCH.max()

# (E) Entanglement advantage for three channel pairs
cases = [
    ("id vs full\ndepolarizing", identity_channel(), depolarizing(1.0)),
    ("id vs\nZ-dephasing", identity_channel(), dephasing(1.0)),
    ("depol(0.6) vs\nampl.damp(0.3)", ch_a, ch_b),
]
single_vals, diamond_vals = [], []
for name, c0, c1 in cases:
    single_vals.append(best_single_input(c0, c1))
    diamond_vals.append(diamond_distance(c0, c1))

elapsed = time.time() - t_start

# ------------------------------------------------------------
# 5. Console report
# ------------------------------------------------------------
print("Solver:", SOLVER)
print("=" * 64)
print("(A) identity vs depolarizing(p)   [analytic: 3p/2]")
print(f"{'p':>6} {'SDP':>10} {'analytic':>10} {'abs.err':>10}")
for p, s, a in zip(ps, dn_depol, an_depol):
    print(f"{p:6.2f} {s:10.6f} {a:10.6f} {abs(s - a):10.2e}")
print(f"max error = {np.max(np.abs(dn_depol - an_depol)):.2e}")
print("=" * 64)
print("(B) identity vs Z-rotation(theta)   [analytic: 2|sin(theta/2)|]")
print(f"{'theta/pi':>9} {'SDP':>10} {'analytic':>10} {'abs.err':>10}")
for t, s, a in zip(thetas[::4], dn_rot[::4], an_rot[::4]):
    print(f"{t / np.pi:9.3f} {s:10.6f} {a:10.6f} {abs(s - a):10.2e}")
print(f"max error = {np.max(np.abs(dn_rot - an_rot)):.2e}")
print("=" * 64)
print(f"(C) depolarizing vs amplitude damping on a {n_grid}x{n_grid} grid")
imax = np.unravel_index(np.argmax(surface), surface.shape)
print(f"max diamond distance = {surface.max():.6f} at p = {p_grid[imax[0]]:.3f}, gamma = {g_grid[imax[1]]:.3f}")
print(f"Helstrom success probability range: {success.min():.4f} .. {success.max():.4f}")
print("=" * 64)
print(f"(D) fixed pair p = {p_fix}, gamma = {g_fix}")
print(f"diamond norm (SDP, entangled input)  = {dn_fix:.6f}")
print(f"best single-input trace norm (sphere) = {best_fix:.6f}")
print("=" * 64)
print("(E) single input vs ancilla-assisted (diamond norm)")
for (name, _, _), s, d in zip(cases, single_vals, diamond_vals):
    print(f"{name.replace(chr(10), ' '):<34} single = {s:.6f}   diamond = {d:.6f}")
print("=" * 64)
print(f"Total elapsed time: {elapsed:.1f} s")

# ------------------------------------------------------------
# 6. One combined figure
# ------------------------------------------------------------
plt.style.use("dark_background")
fig = plt.figure(figsize=(21, 12))
gs = GridSpec(2, 3, figure=fig, hspace=0.32, wspace=0.28)

ax1 = fig.add_subplot(gs[0, 0])
ax1.plot(ps, an_depol, "-", color="#00e5ff", lw=2, label="analytic $3p/2$")
ax1.plot(ps, dn_depol, "o", color="#ff4081", ms=8, label="SDP")
ax1.set_title("(A) id vs depolarizing channel")
ax1.set_xlabel("depolarizing probability $p$")
ax1.set_ylabel(r"$\|\mathrm{id}-\mathcal{D}_p\|_\diamond$")
ax1.grid(alpha=0.25)
ax1.legend()

ax2 = fig.add_subplot(gs[0, 1])
ax2.plot(thetas / np.pi, an_rot, "-", color="#00e5ff", lw=2, label=r"analytic")
ax2.plot(thetas / np.pi, dn_rot, "o", color="#ffea00", ms=7, label="SDP")
ax2.set_title("(B) id vs $Z$-rotation")
ax2.set_xlabel(r"rotation angle $\theta/\pi$")
ax2.set_ylabel(r"$\|\mathrm{id}-\mathcal{R}_\theta\|_\diamond$")
ax2.grid(alpha=0.25)
ax2.legend()

PG, GG = np.meshgrid(p_grid, g_grid, indexing="ij")
ax3 = fig.add_subplot(gs[0, 2], projection="3d")
surf = ax3.plot_surface(PG, GG, surface, cmap="plasma", edgecolor="none", alpha=0.95)
ax3.set_title("(C) depolarizing vs amplitude damping")
ax3.set_xlabel("$p$")
ax3.set_ylabel(r"$\gamma$")
ax3.set_zlabel(r"$\|\cdot\|_\diamond$")
ax3.view_init(elev=28, azim=-125)
fig.colorbar(surf, ax=ax3, shrink=0.6, pad=0.1)

ax4 = fig.add_subplot(gs[1, 0], projection="3d")
Xs = np.sin(TH) * np.cos(PH)
Ys = np.sin(TH) * np.sin(PH)
Zs = np.cos(TH)
norm_c = plt.Normalize(vmin=BLOCH.min(), vmax=BLOCH.max())
ax4.plot_surface(Xs, Ys, Zs, facecolors=plt.cm.viridis(norm_c(BLOCH)), rstride=1, cstride=1,
                 linewidth=0, antialiased=False, shade=False)
m = plt.cm.ScalarMappable(norm=norm_c, cmap="viridis")
m.set_array(BLOCH)
fig.colorbar(m, ax=ax4, shrink=0.6, pad=0.1, label="trace norm of output difference")
ax4.set_title(f"(D) single-input distinguishability on Bloch sphere\n"
              f"max = {best_fix:.3f}  vs  diamond norm = {dn_fix:.3f}")
ax4.set_xlabel("x")
ax4.set_ylabel("y")
ax4.set_zlabel("z")
ax4.set_box_aspect((1, 1, 1))

ax5 = fig.add_subplot(gs[1, 1])
cs = ax5.contourf(PG, GG, success, levels=20, cmap="magma")
ax5.contour(PG, GG, success, levels=10, colors="white", linewidths=0.5, alpha=0.5)
fig.colorbar(cs, ax=ax5)
ax5.plot([p_fix], [g_fix], "c*", ms=16)
ax5.set_title(r"(E) optimal success probability $\frac{1}{2}+\frac{1}{4}\|\Phi_0-\Phi_1\|_\diamond$")
ax5.set_xlabel("$p$")
ax5.set_ylabel(r"$\gamma$")

ax6 = fig.add_subplot(gs[1, 2])
x = np.arange(len(cases))
w = 0.36
ax6.bar(x - w / 2, single_vals, w, color="#26c6da", label="best single input (no ancilla)")
ax6.bar(x + w / 2, diamond_vals, w, color="#ff4081", label="diamond norm (ancilla-assisted)")
ax6.set_xticks(x)
ax6.set_xticklabels([c[0] for c in cases])
ax6.set_ylabel("distinguishability")
ax6.set_title("(F) entanglement advantage")
ax6.grid(alpha=0.25, axis="y")
ax6.legend(loc="upper left")
ax6.set_ylim(0, 2.6)

fig.suptitle("Quantum Channel Discrimination and the Diamond Norm (SDP)", fontsize=20, y=0.97)
plt.show()

6. Code walkthrough

6.1 Channels and Choi matrices

Each channel is a Python function acting on a $2\times2$ matrix. Because every channel is linear, the same function works on density matrices and on the matrix units $|i\rangle\langle j|$ that appear in the Choi construction. choi() builds

$$
J(\Phi)=\sum_{i,j}|i\rangle\langle j|\otimes\Phi(|i\rangle\langle j|)
$$

with np.kron(E, channel(E)). The input system is the first tensor factor and the output system is the second. The SDP uses the same ordering, which is why a partial trace over the second factor is the right operation.

6.2 The SDP in `diamond_norm`

Y0 and Y1 are Hermitian CVXPY variables of size $4\times4$ for a qubit.
block is the $8\times8$ matrix $\begin{pmatrix}Y_0&-J\-J^\dagger&Y_1\end{pmatrix}$, and block >> 0 is the positive-semidefinite constraint.
The partial trace over the output is written without any special atom. For each output basis vector $|k\rangle$ we form the constant matrix $A_k=\mathbb 1_{\text{in}}\otimes\langle k|$, and then
$$
\mathrm{Tr}_{\text{out}},Y=\sum_k A_k,Y,A_k^{\mathsf T}.
$$
This keeps the model purely affine in the variables and works for complex Hermitian matrices.
ptrace_out(Y0) << t0*I encodes $\mathrm{Tr}_{\text{out}}Y_0\preceq t_0\mathbb 1$, the standard way to express a spectral-norm bound inside an SDP.
The objective is $\tfrac12(t_0+t_1)$. The solver is Clarabel, an interior-point method that reaches high accuracy on small SDPs, with SCS as a fallback.

6.3 Speed

Each SDP involves only $4\times4$ and $8\times8$ matrices. The optimal values are obtained to roughly $10^{-8}$ accuracy by an interior-point solver in a few milliseconds each, so the whole study of about 280 SDPs finishes in seconds. The Bloch-sphere scan uses only $2\times2$ eigenvalue computations. The 3D map therefore needs no approximation scheme, only a modest $13\times13$ grid that is enough for a smooth surface.

6.4 Comparison tools

trace_norm() computes $|M|_1$ from the eigenvalues of the Hermitian matrix $M$. best_single_input() scans pure input states $|\psi(\theta,\phi)\rangle$ over the Bloch sphere and returns $\max_\psi|\Phi_0(\psi)-\Phi_1(\psi)|_1$, the best any ancilla-free strategy can do.

7. Results

Execution result image

Console output

Solver: CLARABEL
================================================================
(A) identity vs depolarizing(p)   [analytic: 3p/2]
     p        SDP   analytic    abs.err
  0.00   0.000000   0.000000   3.77e-11
  0.10   0.150000   0.150000   2.78e-10
  0.20   0.300000   0.300000   1.75e-10
  0.30   0.450000   0.450000   1.79e-10
  0.40   0.600000   0.600000   1.14e-10
  0.50   0.750000   0.750000   2.19e-09
  0.60   0.900000   0.900000   5.60e-09
  0.70   1.050000   1.050000   1.17e-08
  0.80   1.200000   1.200000   8.65e-09
  0.90   1.350000   1.350000   3.62e-09
  1.00   1.500000   1.500000   2.08e-09
max error = 1.17e-08
================================================================
(B) identity vs Z-rotation(theta)   [analytic: 2|sin(theta/2)|]
 theta/pi        SDP   analytic    abs.err
    0.000   0.000000   0.000000   3.77e-11
    0.333   1.000000   1.000000   6.22e-10
    0.667   1.732051   1.732051   5.52e-10
    1.000   2.000000   2.000000   2.20e-09
    1.333   1.732051   1.732051   5.52e-10
    1.667   1.000000   1.000000   6.22e-10
    2.000   0.000000   0.000000   6.39e-12
max error = 2.20e-09
================================================================
(C) depolarizing vs amplitude damping on a 13x13 grid
max diamond distance = 2.000000 at p = 0.000, gamma = 1.000
Helstrom success probability range: 0.5000 .. 1.0000
================================================================
(D) fixed pair p = 0.6, gamma = 0.3
diamond norm (SDP, entangled input)  = 0.665150
best single-input trace norm (sphere) = 0.600936
================================================================
(E) single input vs ancilla-assisted (diamond norm)
id vs full depolarizing            single = 1.000000   diamond = 1.500000
id vs Z-dephasing                  single = 2.000000   diamond = 2.000000
depol(0.6) vs ampl.damp(0.3)       single = 0.600941   diamond = 0.665150
================================================================
Total elapsed time: 20.2 s

8. Reading the figure

Panel (A): identity versus depolarizing. The SDP values lie on the line $3p/2$. At $p=1$ the value $3/2$ is the trace distance between the maximally entangled state and the maximally mixed state on two qubits, $|,|\Phi^+\rangle\langle\Phi^+|-\mathbb 1/4|_1=\tfrac34+3\cdot\tfrac14$. The match with the analytic line shows that the SDP, the Choi convention and the partial trace are all correct.

Panel (B): identity versus $Z$-rotation. The SDP reproduces $2|\sin(\theta/2)|$. The curve peaks at $\theta=\pi$, where $U_\theta^\dagger$ has eigenvalues $\pm i$, so the convex hull of the eigenvalues contains the origin. There the two unitaries are perfectly distinguishable and the norm reaches its maximum value $2$. At $\theta=2\pi$ the unitary equals $-\mathbb 1$, which is the identity channel up to a global phase, so the norm returns to $0$.

Panel (C): the 3D landscape. For depolarizing versus amplitude damping, the surface is zero only where the two channels coincide, namely the corner $p=\gamma=0$ where both are the identity. It rises toward the corner $p=0,\ \gamma=1$, where the norm reaches its maximum of $2$. There the identity channel is compared with the channel that resets every state to $|0\rangle$, and these two outputs can be made orthogonal. The surface is a continuous, convex-looking sheet, as expected, since the diamond norm is a norm and both channels depend affinely on $p$ and $\gamma$ (the Kraus-form amplitude damping channel is affine in $\gamma$ at the level of its action on density matrices).

Panel (D): single-input distinguishability on the Bloch sphere. Each point of the sphere is a pure input state, coloured by $|\Phi_0(\psi)-\Phi_1(\psi)|_1$ for $p=0.6,\ \gamma=0.3$. The maximum over the sphere is smaller than the diamond norm printed in the title. The remaining difference can be obtained only by sending half of an entangled state and measuring the pair jointly.

Panel (E): optimal success probability. The contour map shows $\tfrac12+\tfrac14|\Phi_0-\Phi_1|_\diamond$ over the $(p,\gamma)$ plane. Success ranges from $1/2$, which is a coin flip, at the coincident corner up to $1$ at the corner where the channels are perfectly distinguishable. The star marks the pair studied in panels (D) and (F).

Panel (F): the entanglement advantage. Three cases appear side by side.

Identity versus full depolarizing: the best ancilla-free value is $1$, but the diamond norm is $3/2$. This is the clean example of a gap created by entanglement.
Identity versus $Z$-dephasing with $p=1$: both values equal $2$. A single input such as $|+\rangle$ is already perfectly sufficient, since $|+\rangle$ and $|-\rangle$ are orthogonal, so entanglement brings no improvement.
Depolarizing versus amplitude damping at $(0.6,0.3)$: the diamond norm exceeds the best single-input value, a smaller but genuine advantage.

9. Takeaways

The diamond norm is the correct figure of merit for channel discrimination, because it accounts for entangled inputs and adaptive measurements.
A semidefinite program with only two small positive-semidefinite blocks computes it exactly for qubit channels.
Closed-form cases (depolarizing noise and unitary rotations) agree with the SDP to numerical precision, which makes them good unit tests for any implementation.
Optimal discrimination can require an ancilla. The gap between the best single-input value and the diamond norm is the quantitative benefit of entanglement.

The same code extends to larger systems by changing the dimension arguments d_in and d_out and supplying the Choi matrix of the difference of two channels.