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.

Computing Quantum Channel Capacities in Python

October 6, 2026

Maximizing the Holevo Information of the Amplitude Damping Channel

How many classical bits, or how many qubits, can one use of a noisy quantum channel carry? In this post we take a concrete channel, the qubit amplitude damping channel, and compute its classical capacity by maximizing the Holevo information. We also compute its quantum capacity from the coherent information. The result is a table, six graphs including three 3D plots, and the full code.

1. The Example Problem

The amplitude damping channel models spontaneous emission: an excited qubit $|1\rangle$ decays to $|0\rangle$ with probability $\gamma$. It has two 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},\qquad
\mathcal N_\gamma(\rho)=K_0\rho K_0^\dagger+K_1\rho K_1^\dagger .
$$

Writing a qubit state as $\rho=\tfrac12(I+xX+yY+zZ)$, the channel acts on the Bloch vector as

$$
(x,y,z);\longmapsto;\bigl(\sqrt{1-\gamma},x,\ \sqrt{1-\gamma},y,\ (1-\gamma)z+\gamma\bigr).
$$

The Bloch sphere is squeezed into an ellipsoid and shifted towards the north pole $|0\rangle$.

Task. For $\gamma\in[0,1]$, compute

the classical capacity $C(\mathcal N_\gamma)$, by maximizing the Holevo information,
the quantum capacity $Q(\mathcal N_\gamma)$, by maximizing the coherent information.

2. Theory in Brief

Classical capacity

For an ensemble ${p_k,\rho_k}$ the Holevo information is

$$
\chi\bigl({p_k,\rho_k}\bigr)=S\Bigl(\sum_k p_k,\mathcal N(\rho_k)\Bigr)-\sum_k p_k,S\bigl(\mathcal N(\rho_k)\bigr),
$$

where $S(\rho)=-\mathrm{Tr},\rho\log_2\rho$. The Holevo capacity is its maximum over all ensembles:

$$
C(\mathcal N)=\max_{(p_k,\rho_k)}\chi\bigl((p_k,\rho_k)\bigr).
$$

For the amplitude damping channel this quantity is additive, so the single-letter maximization gives the true classical capacity. For a qubit, four pure states suffice in the ensemble. A qubit with Bloch vector of length $r$ has entropy

$$
S=h_2!\left(\frac{1+r}{2}\right),\qquad h_2(x)=-x\log_2x-(1-x)\log_2(1-x).
$$

A symmetric ensemble

The channel commutes with rotations about the $z$-axis. Averaging an ensemble over such rotations keeps every output entropy unchanged and cannot decrease the entropy of the average output. So it suffices to look at a north pole with weight $1-q$ plus a uniform ring at polar angle $\theta$ with weight $q$. Define

$$
z_\theta=(1-\gamma)\cos\theta+\gamma,\qquad
r_\theta=\sqrt{(1-\gamma)\sin^2\theta+z_\theta^{,2}},\qquad
\bar z=(1-q)+q,z_\theta .
$$

Then

$$
\chi(\theta,q)=h_2!\left(\frac{1+\bar z}{2}\right)-q,h_2!\left(\frac{1+r_\theta}{2}\right),
\qquad
C=\max_{\theta,q}\chi(\theta,q).
$$

A natural but suboptimal choice

Restricting to the orthogonal inputs ${|0\rangle,|1\rangle}$, with $|1\rangle$ sent with probability $p$, gives a lower bound:

$$
\chi_{\rm diag}(p,\gamma)=h_2\bigl((1-\gamma)p\bigr)-p,h_2(\gamma).
$$

Non-orthogonal inputs do better, and we will see by how much.

Quantum capacity

The amplitude damping channel is degradable for $\gamma\le\tfrac12$ and antidegradable for $\gamma\ge\tfrac12$. Therefore

$$
Q(\mathcal N_\gamma)=
\begin{cases}
\displaystyle\max_{p\in[0,1]}\Bigl[h_2\bigl((1-\gamma)p\bigr)-h_2(\gamma p)\Bigr], & \gamma\le \tfrac12,\[2mm]
0, & \gamma\ge\tfrac12 .
\end{cases}
$$

3. The Code

Everything is in a single script. Run it as one cell.

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

np.set_printoptions(precision=6, suppress=True)
EPS = 1e-15

# ------------------------------------------------------------
# 1. Basic tools: entropy, Kraus operators, Bloch representation
# ------------------------------------------------------------
def h2(x):
    x = np.clip(np.asarray(x, dtype=float), EPS, 1.0 - EPS)
    return -x * np.log2(x) - (1.0 - x) * np.log2(1.0 - x)

def vn_entropy(rho):
    ev = np.clip(np.linalg.eigvalsh(rho), EPS, None)
    return float(-np.sum(ev * np.log2(ev)))

def kraus_ad(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 [K0, K1]

def channel_rho(rho, g):
    return sum(K @ rho @ K.conj().T for K in kraus_ad(g))

def bloch_to_rho(v):
    x, y, z = v
    return 0.5 * np.array([[1 + z, x - 1j * y], [x + 1j * y, 1 - z]], dtype=complex)

def ad_bloch(v, g):
    v = np.atleast_2d(v)
    s = np.sqrt(1.0 - g)
    return np.column_stack([s * v[:, 0], s * v[:, 1], (1.0 - g) * v[:, 2] + g])

# ------------------------------------------------------------
# 2. Holevo quantity: slow (density matrices) vs fast (Bloch vectors)
# ------------------------------------------------------------
def holevo_naive(vs, p, g):
    outs = [channel_rho(bloch_to_rho(v), g) for v in vs]
    avg = sum(pk * o for pk, o in zip(p, outs))
    return vn_entropy(avg) - sum(pk * vn_entropy(o) for pk, o in zip(p, outs))

def holevo_bloch(vs, p, g):
    out = ad_bloch(vs, g)
    avg = (p[:, None] * out).sum(axis=0)
    s_avg = h2((1.0 + np.linalg.norm(avg)) / 2.0)
    s_each = h2((1.0 + np.linalg.norm(out, axis=1)) / 2.0)
    return float(s_avg - np.dot(p, s_each))

rng = np.random.default_rng(0)
vs_test = rng.normal(size=(4, 3))
vs_test = vs_test / np.linalg.norm(vs_test, axis=1, keepdims=True)
p_test = rng.dirichlet(np.ones(4))
g_test = 0.3
n_rep = 3000

t0 = time.perf_counter()
for _ in range(n_rep):
    a = holevo_naive(vs_test, p_test, g_test)
t_naive = time.perf_counter() - t0

t0 = time.perf_counter()
for _ in range(n_rep):
    b = holevo_bloch(vs_test, p_test, g_test)
t_fast = time.perf_counter() - t0

print("=== Consistency and speed check ===")
print(f"Holevo (density matrix) : {a:.12f}")
print(f"Holevo (Bloch vector)   : {b:.12f}")
print(f"Absolute difference     : {abs(a - b):.2e}")
print(f"Time naive / fast       : {t_naive:.3f} s / {t_fast:.3f} s  (speed-up x{t_naive / t_fast:.1f})")

# ------------------------------------------------------------
# 3. General optimization over 4-state ensembles (Davies bound for qubits)
# ------------------------------------------------------------
def unpack(x, K):
    th, ph, w = x[:K], x[K:2 * K], x[2 * K:3 * K]
    vs = np.column_stack([np.sin(th) * np.cos(ph), np.sin(th) * np.sin(ph), np.cos(th)])
    e = np.exp(w - w.max())
    return vs, e / e.sum()

def max_holevo_general(g, K=4, n_restarts=8, seed=1):
    rg = np.random.default_rng(seed)
    best_val, best_x = -1.0, None
    for _ in range(n_restarts):
        x0 = np.concatenate([rg.uniform(0, np.pi, K), rg.uniform(0, 2 * np.pi, K), rg.normal(size=K)])
        obj = lambda x: -holevo_bloch(*unpack(x, K), g)
        res = minimize(obj, x0, method="Nelder-Mead",
                       options={"maxiter": 4000, "xatol": 1e-10, "fatol": 1e-13})
        if -res.fun > best_val:
            best_val, best_x = -res.fun, res.x
    return best_val, best_x

# ------------------------------------------------------------
# 4. Symmetric ensemble: north pole + uniform ring at polar angle theta
# ------------------------------------------------------------
def chi_ring(theta, q, g):
    zr = (1.0 - g) * np.cos(theta) + g
    rr = np.sqrt((1.0 - g) * np.sin(theta) ** 2 + zr ** 2)
    zavg = (1.0 - q) + q * zr
    return h2((1.0 + np.abs(zavg)) / 2.0) - q * h2((1.0 + rr) / 2.0)

TH = np.linspace(0.0, np.pi, 181)
QQ = np.linspace(0.0, 1.0, 101)
TH_G, QQ_G = np.meshgrid(TH, QQ)

def max_holevo_ring(g):
    grid = chi_ring(TH_G, QQ_G, g)
    i, j = np.unravel_index(np.argmax(grid), grid.shape)
    x0 = np.array([TH_G[i, j], QQ_G[i, j]])
    res = minimize(lambda x: -chi_ring(x[0], x[1], g), x0, method="L-BFGS-B",
                   bounds=[(0.0, np.pi), (0.0, 1.0)])
    if -res.fun >= grid[i, j]:
        return -res.fun, res.x[0], res.x[1]
    return grid[i, j], x0[0], x0[1]

# ------------------------------------------------------------
# 5. Restricted lower bound {|0>,|1>} and quantum capacity
# ------------------------------------------------------------
def chi_diag(p, g):
    return h2((1.0 - g) * p) - p * h2(g)

def coh_info(p, g):
    return h2((1.0 - g) * p) - h2(g * p)

def C_diag(g):
    r = minimize_scalar(lambda p: -chi_diag(p, g), bounds=(0.0, 1.0),
                        method="bounded", options={"xatol": 1e-12})
    return -r.fun

def Q_cap(g):
    if g >= 0.5:
        return 0.0, 0.0
    r = minimize_scalar(lambda p: -coh_info(p, g), bounds=(0.0, 1.0),
                        method="bounded", options={"xatol": 1e-12})
    return max(-r.fun, 0.0), r.x

# ------------------------------------------------------------
# 6. Table of results
# ------------------------------------------------------------
gammas_tab = np.array([0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8])
rows = []
t0 = time.perf_counter()
for g in gammas_tab:
    c_gen, _ = max_holevo_general(g)
    c_ring, th_s, q_s = max_holevo_ring(g)
    rows.append((g, c_gen, c_ring, C_diag(g), Q_cap(g)[0], th_s, q_s))
t_tab = time.perf_counter() - t0

print("\n=== Capacities of the amplitude damping channel (in bits per channel use) ===")
print(f"{'gamma':>6} | {'C general':>10} | {'C ring':>10} | {'C {|0>,|1>}':>11} | {'Q':>8} | {'theta*/pi':>9} | {'q*':>6}")
print("-" * 78)
for g, cg, cr, cd, qq, th_s, q_s in rows:
    print(f"{g:6.2f} | {cg:10.6f} | {cr:10.6f} | {cd:11.6f} | {qq:8.6f} | {th_s / np.pi:9.4f} | {q_s:6.4f}")
print(f"\nTotal optimization time for the table: {t_tab:.2f} s")

# ------------------------------------------------------------
# 7. Curves over the whole range of gamma
# ------------------------------------------------------------
gam = np.linspace(0.0, 1.0, 101)
C_ring_curve = np.zeros_like(gam)
th_curve = np.zeros_like(gam)
q_curve = np.zeros_like(gam)
for k, g in enumerate(gam):
    C_ring_curve[k], th_curve[k], q_curve[k] = max_holevo_ring(g)
C_diag_curve = np.array([C_diag(g) for g in gam])
Q_curve = np.array([Q_cap(g)[0] for g in gam])
p_Q_curve = np.array([Q_cap(g)[1] for g in gam])

# ------------------------------------------------------------
# 8. One figure with six panels
# ------------------------------------------------------------
g_demo = 0.3
c_demo, th_demo, q_demo = max_holevo_ring(g_demo)

fig = plt.figure(figsize=(20, 12))
gs = gridspec.GridSpec(2, 3, figure=fig, wspace=0.28, hspace=0.30)

# (a) capacities versus gamma
ax = fig.add_subplot(gs[0, 0])
ax.plot(gam, C_ring_curve, lw=2.5, color="tab:blue", label=r"Classical capacity $C$ (optimized)")
ax.plot(gam, C_diag_curve, lw=2, ls="--", color="tab:orange", label=r"Orthogonal inputs $\{|0\rangle,|1\rangle\}$ only")
ax.plot(gam, Q_curve, lw=2.5, color="tab:green", label=r"Quantum capacity $Q$")
ax.axvline(0.5, color="gray", ls=":", lw=1)
ax.set_xlabel(r"Damping parameter $\gamma$")
ax.set_ylabel("Capacity [bits / use]")
ax.set_title("(a) Capacities of the amplitude damping channel")
ax.legend(loc="upper right", fontsize=9)
ax.grid(alpha=0.3)

# (b) 3D surface of the Holevo quantity
ax = fig.add_subplot(gs[0, 1], projection="3d")
Z = chi_ring(TH_G, QQ_G, g_demo)
ax.plot_surface(TH_G / np.pi, QQ_G, Z, cmap="viridis", alpha=0.9, linewidth=0, antialiased=True)
ax.scatter([th_demo / np.pi], [q_demo], [c_demo + 0.02], color="red", s=80, depthshade=False)
ax.set_xlabel(r"Ring angle $\theta/\pi$")
ax.set_ylabel(r"Ring weight $q$")
ax.set_zlabel(r"$\chi$ [bits]")
ax.set_title(rf"(b) Holevo quantity, $\gamma={g_demo}$ (red = maximum)")
ax.view_init(elev=28, azim=-125)

# (c) 3D surface of the coherent information
ax = fig.add_subplot(gs[0, 2], projection="3d")
P_g, G_g = np.meshgrid(np.linspace(0, 1, 80), np.linspace(0, 1, 80))
Ic = coh_info(P_g, G_g)
ax.plot_surface(P_g, G_g, Ic, cmap="plasma", alpha=0.9, linewidth=0, antialiased=True)
mask = gam < 0.5
ax.plot(p_Q_curve[mask], gam[mask], Q_curve[mask] + 0.01, color="black", lw=3)
ax.set_xlabel(r"Input weight $p$")
ax.set_ylabel(r"Damping $\gamma$")
ax.set_zlabel(r"$I_c$ [bits]")
ax.set_title(r"(c) Coherent information (black = ridge $=Q$)")
ax.view_init(elev=28, azim=-60)

# (d) optimal ensemble parameters
ax = fig.add_subplot(gs[1, 0])
ax.plot(gam, th_curve / np.pi, lw=2.5, color="tab:purple", label=r"Optimal ring angle $\theta^*/\pi$")
ax.plot(gam, q_curve, lw=2.5, color="tab:red", label=r"Optimal ring weight $q^*$")
ax.set_xlabel(r"Damping parameter $\gamma$")
ax.set_ylabel("Parameter value")
ax.set_title("(d) Optimal input ensemble")
ax.legend(fontsize=9)
ax.grid(alpha=0.3)

# (e) Bloch sphere picture
ax = fig.add_subplot(gs[1, 1], projection="3d")
u, v = np.mgrid[0:2 * np.pi:60j, 0:np.pi:30j]
ax.plot_wireframe(np.cos(u) * np.sin(v), np.sin(u) * np.sin(v), np.cos(v), color="lightgray", alpha=0.35, lw=0.5)
ex = np.sqrt(1 - g_demo) * np.cos(u) * np.sin(v)
ey = np.sqrt(1 - g_demo) * np.sin(u) * np.sin(v)
ez = (1 - g_demo) * np.cos(v) + g_demo
ax.plot_surface(ex, ey, ez, color="tab:cyan", alpha=0.35, linewidth=0)
phi = np.linspace(0, 2 * np.pi, 200)
ring_in = np.column_stack([np.sin(th_demo) * np.cos(phi), np.sin(th_demo) * np.sin(phi), np.cos(th_demo) * np.ones_like(phi)])
ring_out = ad_bloch(ring_in, g_demo)
ax.plot(ring_in[:, 0], ring_in[:, 1], ring_in[:, 2], color="tab:red", lw=3, label="Optimal input ring")
ax.plot(ring_out[:, 0], ring_out[:, 1], ring_out[:, 2], color="tab:blue", lw=3, label="Output ring")
ax.scatter([0], [0], [1], color="tab:red", s=70, depthshade=False)
ax.scatter([0], [0], [1], color="tab:blue", s=25, depthshade=False)
ax.set_xlabel("x"); ax.set_ylabel("y"); ax.set_zlabel("z")
ax.set_title(rf"(e) Bloch sphere and channel image, $\gamma={g_demo}$")
ax.legend(loc="upper left", fontsize=8)
ax.view_init(elev=20, azim=35)

# (f) agreement between the two optimizers
ax = fig.add_subplot(gs[1, 2])
diff = np.array([abs(r[1] - r[2]) for r in rows]) + 1e-12
ax.semilogy(gammas_tab, diff, "o-", lw=2, color="tab:brown")
ax.set_xlabel(r"Damping parameter $\gamma$")
ax.set_ylabel(r"$|C_{\mathrm{general}}-C_{\mathrm{ring}}|$")
ax.set_title("(f) General 4-state search vs. symmetric ring ensemble")
ax.grid(alpha=0.3, which="both")

plt.show()

4. Execution Results

Console output

=== Consistency and speed check ===
Holevo (density matrix) : 0.521417030487
Holevo (Bloch vector)   : 0.521417030487
Absolute difference     : 4.44e-16
Time naive / fast       : 0.794 s / 0.149 s  (speed-up x5.3)

=== Capacities of the amplitude damping channel (in bits per channel use) ===
 gamma |  C general |     C ring | C {|0>,|1>} |        Q | theta*/pi |     q*
------------------------------------------------------------------------------
  0.00 |   1.000000 |   1.000000 |    1.000000 | 1.000000 |    1.0000 | 0.5000
  0.10 |   0.840497 |   0.840497 |    0.762848 | 0.709418 |    0.4702 | 1.0000
  0.20 |   0.731645 |   0.731645 |    0.618231 | 0.506215 |    0.4571 | 1.0000
  0.30 |   0.638329 |   0.638329 |    0.503692 | 0.327955 |    0.4485 | 1.0000
  0.40 |   0.552957 |   0.552957 |    0.406787 | 0.161480 |    0.4425 | 1.0000
  0.50 |   0.471729 |   0.471729 |    0.321928 | 0.000000 |    0.4384 | 1.0000
  0.60 |   0.392017 |   0.392017 |    0.245986 | 0.000000 |    0.4359 | 1.0000
  0.80 |   0.226742 |   0.226742 |    0.113594 | 0.000000 |    0.4362 | 1.0000

Total optimization time for the table: 18.29 s

Graph output

5. Code Walkthrough

Section 1: Basic tools

h2 is the binary entropy $h_2(x)$. The argument is clipped to $[10^{-15},1-10^{-15}]$ so that $0\log 0$ never produces nan. vn_entropy computes $S(\rho)$ from the eigenvalues of a density matrix. kraus_ad and channel_rho implement $\mathcal N_\gamma(\rho)=\sum_i K_i\rho K_i^\dagger$ literally. bloch_to_rho converts a Bloch vector to a $2\times2$ density matrix. ad_bloch applies the channel directly to Bloch vectors, using the affine map

$$
(x,y,z)\mapsto(\sqrt{1-\gamma},x,\ \sqrt{1-\gamma},y,\ (1-\gamma)z+\gamma).
$$

Section 2: Slow and fast Holevo information

holevo_naive follows the definition literally: it builds density matrices, applies the Kraus operators and diagonalizes. holevo_bloch does the same on Bloch vectors. The output entropy of each state is just $h_2((1+r)/2)$ with $r$ the vector length, so no matrices and no eigenvalue decompositions are needed. Because the optimizer calls this function tens of thousands of times, the fast version matters. The script runs both on a random four-state ensemble, prints that the values agree to machine precision, and prints the measured speed-up.

Section 3: General numerical maximization

An ensemble of $K=4$ pure states is parametrized by $3K$ real numbers: polar angles, azimuthal angles, and unnormalized weights. unpack turns them into unit vectors and a probability vector via a softmax, so every point of the search space is a valid ensemble and no constraints are needed. max_holevo_general runs Nelder–Mead from eight random starting points and keeps the best value, which protects against local maxima. This is a brute-force check that does not assume anything about the optimal structure.

Section 4: The symmetric ring ensemble

chi_ring is the closed-form $\chi(\theta,q)$ from Section 2. It is fully vectorized, so it evaluates on a whole $181\times101$ grid in one call. max_holevo_ring finds the best grid point and then refines it with the bounded optimizer L-BFGS-B. The grid prevents the local optimizer from starting in the wrong basin. This two-parameter search is far cheaper than the twelve-parameter one and is used for the full curves.

Section 5: Lower bound and quantum capacity

chi_diag is the orthogonal-input formula $h_2((1-\gamma)p)-p,h_2(\gamma)$, maximized over $p$ with a bounded scalar optimizer. coh_info is $I_c(p,\gamma)=h_2((1-\gamma)p)-h_2(\gamma p)$. Q_cap returns the maximum and the maximizer for $\gamma<\tfrac12$ and $0$ otherwise, in line with the degradable/antidegradable split.

Section 6: The table

For eight values of $\gamma$ the script prints the general 4-state result, the ring result, the orthogonal-input lower bound, the quantum capacity, and the optimal ring parameters $(\theta^*/\pi,q^*)$. The first two columns should agree to the printed precision.

Section 7: Full curves

The ring optimizer is run on 101 values of $\gamma$. This is cheap precisely because of the vectorized grid. The same loop collects the optimal $(\theta^*,q^*)$ and the optimal input weight for $Q$.

Section 8: The figure

All six panels sit in one figure, drawn with a GridSpec. Panels (b), (c) and (e) use projection="3d". Panel (e) draws the unit sphere as a light wireframe and the channel’s image as a translucent ellipsoid.

6. Reading the Results

Panel (a): the three capacities

The blue curve is the classical capacity $C$, the dashed orange curve is the orthogonal-input lower bound, and the green curve is the quantum capacity $Q$. At $\gamma=0$ the channel is the identity and all three equal 1 bit. As $\gamma$ grows, $Q$ drops fastest and hits exactly zero at $\gamma=\tfrac12$ (dotted line), where the channel becomes antidegradable. The classical capacity remains positive past that point. Even at $\gamma=0.8$ some classical information gets through, while no quantum information can.

The blue curve lies strictly above the orange one for every $\gamma>0$. Using $|0\rangle$ and $|1\rangle$ as the signal states is a good scheme but not the best one.

Panel (b): the Holevo landscape (3D)

For $\gamma=0.3$, the surface shows $\chi$ over the ring angle $\theta/\pi$ and ring weight $q$, with the maximum marked in red. The landscape is smooth with a single clear peak. The peak sits at the edge $q=1$ and at an angle $\theta^*$ below $\pi/2$. The optimal signal is therefore two pure states symmetric about the $z$-axis, tilted toward $|0\rangle$, with overlap $\cos\theta^*>0$. The tilt towards $|0\rangle$ is sensible: $|0\rangle$ is the fixed point of the decay and is the least damaged state.

Panel (c): the coherent information (3D)

The surface $I_c(p,\gamma)$ forms a ridge over the $(p,\gamma)$ plane. The black curve follows the ridge top, which is exactly $Q(\gamma)$ with its maximizer $p^*$. The ridge shrinks as $\gamma\to\tfrac12$ and the surface drops to zero there. Beyond $\gamma=\tfrac12$ the maximum over $p$ would be non-positive, which is why the quantum capacity is set to zero.

Panel (d): the optimal ensemble

For $\gamma>0$ the optimal ring weight is $q^*=1$: the best ensemble contains no pole state at all, only the ring. The optimal angle $\theta^*/\pi$ is around $0.47$ for small $\gamma$ and settles near $0.44$ for larger $\gamma$. At $\gamma=0$ the optimum is the degenerate orthogonal pair (north and south poles, i.e. $\theta=\pi$ with $q=\tfrac12$), which explains the jump at the left edge of the plot.

Panel (e): the geometry (3D)

The gray wireframe is the Bloch sphere. The cyan ellipsoid is its image under the channel at $\gamma=0.3$, shrunk horizontally by $\sqrt{1-\gamma}$, vertically by $1-\gamma$, and lifted towards $|0\rangle$. The red ring is the optimal input ensemble and the blue ring is its output. The output ring is smaller and sits higher, but its points stay far enough from the center of the sphere to remain distinguishable, which is what makes the Holevo quantity large.

Panel (f): sanity check

The unrestricted 4-state search and the symmetric ring ensemble agree to within the optimizer’s tolerance for every tested $\gamma$ (plotted on a log scale, with a $10^{-12}$ floor so that exact agreement is still visible). This supports the symmetry argument: nothing is lost by restricting to the ring family.

7. Conclusion

Maximizing the Holevo information numerically gives the classical capacity of the amplitude damping channel. The optimal signal states are non-orthogonal and tilted towards the fixed point $|0\rangle$, and this beats the naive ${|0\rangle,|1\rangle}$ encoding for every $\gamma>0$. The quantum capacity, from the coherent information, vanishes at $\gamma=\tfrac12$, while the classical capacity stays positive all the way up to $\gamma=1$. Working with Bloch vectors instead of density matrices makes the optimization much faster, and the symmetric ring parametrization reduces a twelve-parameter problem to two parameters without losing anything.

Optimal Quantum Cloning

October 5, 2026

Finding the Channel That Maximizes Fidelity with Python

The no-cloning theorem says that no physical process can turn an unknown quantum state $|\psi\rangle$ into $|\psi\rangle\otimes|\psi\rangle$ exactly. It does not forbid approximate copies. If we accept imperfect copies, the natural question is how good they can be. This article answers that for a qubit by treating the cloner as a convex optimization problem over quantum channels, and then solving a concrete instance numerically.

1. The Problem

A $1\to 2$ cloning machine is a completely positive trace-preserving (CPTP) map

$$
\Lambda:\ \mathcal{L}(\mathbb{C}^2)\ \to\ \mathcal{L}(\mathbb{C}^2\otimes\mathbb{C}^2),
$$

taking one input qubit to two output qubits, $A$ and $B$. For a pure input $|\psi\rangle$, the single-copy fidelities are

$$
F_A(\psi)=\langle\psi|,\mathrm{Tr}_B,\Lambda(|\psi\rangle\langle\psi|),|\psi\rangle,\qquad
F_B(\psi)=\langle\psi|,\mathrm{Tr}_A,\Lambda(|\psi\rangle\langle\psi|),|\psi\rangle .
$$

We want a universal cloner that works equally well for every input state, so we maximize the average over the Haar measure on pure states:

$$
\bar F_\lambda(\Lambda)=\int d\psi,\Big[\lambda,F_A(\psi)+(1-\lambda),F_B(\psi)\Big],\qquad \lambda\in[0,1].
$$

The symmetric case $\lambda=\tfrac12$ is the classic Bužek–Hillery problem. Other values of $\lambda$ give asymmetric cloners, which trade quality of copy $A$ against copy $B$.

Known analytic answer (our benchmark)

For $\lambda=\tfrac12$ the optimum is

$$
F_{\mathrm{opt}}=\frac56\approx 0.8333,
$$

reached by the transformation

$$
|0\rangle\ \mapsto\ \sqrt{\tfrac23},|00\rangle|\uparrow\rangle+\sqrt{\tfrac13},|\Psi^+\rangle|\downarrow\rangle,\qquad
|\Psi^+\rangle=\tfrac{1}{\sqrt2}(|01\rangle+|10\rangle),
$$

together with its mirror image for $|1\rangle$. Each output qubit’s Bloch vector is shrunk by the factor $\eta=2F-1=\tfrac23$.

We will not hard-code this. We solve the optimization numerically and check that it rediscovers $5/6$, then map the whole asymmetric trade-off curve.

2. Turning the Problem into Linear Algebra

Stinespring form

Every channel $\Lambda$ can be written as an isometry $V:\mathbb{C}^2\to\mathbb{C}^2_A\otimes\mathbb{C}^2_B\otimes\mathbb{C}^{r}_E$ followed by a partial trace over the environment $E$:

$$
\Lambda(\rho)=\mathrm{Tr}_E!\left[V\rho V^\dagger\right],\qquad V^\dagger V=\mathbb{1}_2 .
$$

The Choi rank of a map from $\mathbb{C}^2$ to $\mathbb{C}^4$ is at most $8$, so $r=8$ covers every possible cloner.

Replacing the Haar integral by six states

With $P_\psi=|\psi\rangle\langle\psi|$, the fidelity is

$$
F_A(\psi)=\mathrm{Tr}!\left[\big(P_\psi\otimes\mathbb{1}_B\otimes\mathbb{1}_E\big),V P_\psi V^\dagger\right],
$$

which is a degree-2 polynomial in $P_\psi$. The six states at the vertices of the Bloch-sphere octahedron ($\pm x,\pm y,\pm z$) form a projective 2-design. The Haar average of any such polynomial therefore equals the plain average over these six states, exactly:

$$
\bar F_\lambda(V)=\frac16\sum_{s=1}^{6}\mathrm{Tr}!\left[V^\dagger Q_s V R_s\right],
$$

A fast fixed-point algorithm

Because $Q_s\succeq0$ and $R_s\succeq0$, the objective is a convex quadratic form in $V$. Maximizing it over the set of isometries ${V^\dagger V=\mathbb{1}}$ can be done by a power-iteration-like scheme. Define the gradient matrix

$$
G(V)=\sum_{s=1}^{6}Q_s,V,R_s,
$$

and update with the polar factor of $G$. If $G=U\Sigma W^\dagger$ is a thin singular value decomposition, then

$$
V_{n+1}=U,W^\dagger .
$$

This choice maximizes $\mathrm{Re},\mathrm{Tr}(G^\dagger V)$ over all isometries. Convexity guarantees

$$
\bar F_\lambda(V_{n+1})\ \ge\ \bar F_\lambda(V_n),
$$

so the objective rises monotonically. Each iteration costs a few microseconds. No semidefinite-programming solver is needed, which is also the speed-up for this problem: a generic SDP or gradient solver would take far longer for the 21-point trade-off sweep.

3. Reference Strategies

We compare the optimized cloner against two simple strategies.

CNOT cloner. It maps $\alpha|0\rangle+\beta|1\rangle\mapsto\alpha|00\rangle+\beta|11\rangle$ and is not universal. Its fidelity depends on the input state:

$$
F(\theta)=1-\tfrac12\sin^2\theta,\qquad \bar F=\tfrac23 .
$$

Measure-and-prepare. Optimal state estimation from one copy gives

$$
\bar F=\tfrac23 .
$$

4. Full Source Code

import numpy as np
import matplotlib.pyplot as plt

plt.style.use("dark_background")
rng = np.random.default_rng(7)

# ------------------------------------------------------------
# 1. Dimensions
# ------------------------------------------------------------
R_ENV = 8                      # environment dimension (covers every 1 -> 2 qubit channel)
DIM_OUT = 4                    # two output qubits A and B
ROWS = DIM_OUT * R_ENV         # row index of V is (a, b, e) -> (2a + b) * R_ENV + e
F_OPT = 5.0 / 6.0              # analytic optimum for the symmetric universal cloner

# ------------------------------------------------------------
# 2. Six octahedron states (an exact projective 2-design)
# ------------------------------------------------------------
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)
I2 = np.eye(2, dtype=complex)
IR = np.eye(R_ENV, dtype=complex)

bloch_axes = np.array([[1, 0, 0], [-1, 0, 0],
                       [0, 1, 0], [0, -1, 0],
                       [0, 0, 1], [0, 0, -1]], dtype=float)

s2 = 1.0 / np.sqrt(2.0)
octa_states = np.array([[s2, s2], [s2, -s2],
                        [s2, 1j * s2], [s2, -1j * s2],
                        [1.0, 0.0], [0.0, 1.0]], dtype=complex)


def projector_from_bloch(n):
    return 0.5 * (I2 + n[0] * sx + n[1] * sy + n[2] * sz)


R_list = np.array([projector_from_bloch(n) for n in bloch_axes])             # (6, 2, 2)
QA = np.array([np.kron(np.kron(P, I2), IR) for P in R_list])                  # P (x) I_B (x) I_E
QB = np.array([np.kron(np.kron(I2, P), IR) for P in R_list])                  # I_A (x) P (x) I_E
N_STATES = R_list.shape[0]


# ------------------------------------------------------------
# 3. Core routines
# ------------------------------------------------------------
def polar_isometry(G):
    U, _, Wh = np.linalg.svd(G, full_matrices=False)
    return U @ Wh


def random_isometry():
    G = rng.normal(size=(ROWS, 2)) + 1j * rng.normal(size=(ROWS, 2))
    return polar_isometry(G)


def gradient_and_objective(V, Q):
    G = np.einsum('sij,jk,skl->il', Q, V, R_list)
    f = float(np.real(np.trace(V.conj().T @ G))) / N_STATES
    return G, f


def optimize_cloner(lam, n_restarts=6, n_iter=3000, tol=1e-13):
    Q = lam * QA + (1.0 - lam) * QB
    best_V, best_f = None, -1.0
    all_hist = []
    for rs in range(n_restarts):
        V = random_isometry()
        hist = []
        prev = -1.0
        for it in range(n_iter):
            G, f = gradient_and_objective(V, Q)
            hist.append(f)
            if f - prev < tol:
                break
            prev = f
            V = polar_isometry(G)
        _, f_final = gradient_and_objective(V, Q)
        all_hist.append(np.array(hist))
        if f_final > best_f:
            best_f, best_V = f_final, V.copy()
    return best_V, best_f, all_hist


def single_copy_fidelities(V, psi):
    # psi: (N, 2) array of pure states; returns F_A and F_B for each state
    psiT = psi.T
    Psi = (V @ psiT).reshape(2, 2, R_ENV, psi.shape[0])
    phiA = np.einsum('an,aben->ben', psiT.conj(), Psi)
    phiB = np.einsum('bn,aben->aen', psiT.conj(), Psi)
    FA = np.sum(np.abs(phiA) ** 2, axis=(0, 1))
    FB = np.sum(np.abs(phiB) ** 2, axis=(0, 1))
    return FA, FB


def haar_states(n):
    z = rng.normal(size=(n, 2)) + 1j * rng.normal(size=(n, 2))
    return z / np.linalg.norm(z, axis=1, keepdims=True)


# ------------------------------------------------------------
# 4. Reference cloner: CNOT copy (|0> -> |00>, |1> -> |11>)
# ------------------------------------------------------------
V_cnot = np.zeros((ROWS, 2), dtype=complex)
V_cnot[0 * R_ENV + 0, 0] = 1.0
V_cnot[3 * R_ENV + 0, 1] = 1.0

# ------------------------------------------------------------
# 5. Symmetric optimal cloner (lambda = 0.5)
# ------------------------------------------------------------
V_sym, f_sym, hist_sym = optimize_cloner(0.5, n_restarts=6, n_iter=3000)

FA_oct, FB_oct = single_copy_fidelities(V_sym, octa_states)

haar_batch = haar_states(20000)
FA_h, FB_h = single_copy_fidelities(V_sym, haar_batch)
FA_cn_h, FB_cn_h = single_copy_fidelities(V_cnot, haar_batch)

Psi0 = (V_sym @ np.array([1.0, 0.0], dtype=complex)).reshape(2, 2, R_ENV)
rhoA0 = np.einsum('abe,cbe->ac', Psi0, Psi0.conj())
eta = float(np.real(rhoA0[0, 0] - rhoA0[1, 1]))

M = V_sym.reshape(DIM_OUT, R_ENV, 2).transpose(1, 0, 2).reshape(R_ENV, DIM_OUT * 2)
sv = np.linalg.svd(M, compute_uv=False)
kraus_rank = int(np.sum(sv > 1e-8))

print("=== Symmetric optimal cloner (lambda = 0.5) ===")
print(f"Optimized average fidelity      : {f_sym:.12f}")
print(f"Analytic optimum 5/6            : {F_OPT:.12f}")
print(f"Absolute error                  : {abs(f_sym - F_OPT):.3e}")
print(f"Isometry error ||V^dag V - I||  : {np.linalg.norm(V_sym.conj().T @ V_sym - np.eye(2)):.3e}")
print(f"F_A (octahedron average)        : {np.mean(FA_oct):.12f}")
print(f"F_B (octahedron average)        : {np.mean(FB_oct):.12f}")
print(f"Bloch shrinking factor eta      : {eta:.12f}   (analytic 2/3 = {2.0 / 3.0:.12f})")
print(f"Numerical Kraus rank            : {kraus_rank}")
print()
print("=== Haar-random check with 20000 states ===")
print(f"Optimal cloner : mean F_A = {np.mean(FA_h):.9f}, mean F_B = {np.mean(FB_h):.9f}, spread of F_A = {np.ptp(FA_h):.3e}")
print(f"CNOT cloner    : mean F_A = {np.mean(FA_cn_h):.9f}, min = {np.min(FA_cn_h):.4f}, max = {np.max(FA_cn_h):.4f}")

# ------------------------------------------------------------
# 6. Asymmetric trade-off curve
# ------------------------------------------------------------
lam_grid = np.linspace(0.0, 1.0, 21)
FA_curve = np.zeros_like(lam_grid)
FB_curve = np.zeros_like(lam_grid)
for i, lam in enumerate(lam_grid):
    V_l, _, _ = optimize_cloner(lam, n_restarts=4, n_iter=4000)
    fa, fb = single_copy_fidelities(V_l, octa_states)
    FA_curve[i] = np.mean(fa)
    FB_curve[i] = np.mean(fb)

print()
print("=== Asymmetric trade-off (selected weights) ===")
print(f"{'lambda':>8} {'F_A':>12} {'F_B':>12}")
for idx in [0, 5, 10, 15, 20]:
    print(f"{lam_grid[idx]:8.2f} {FA_curve[idx]:12.8f} {FB_curve[idx]:12.8f}")

# ------------------------------------------------------------
# 7. Visualization (one figure)
# ------------------------------------------------------------
theta = np.linspace(0.0, np.pi, 50)
phi = np.linspace(0.0, 2.0 * np.pi, 100)
TH, PH = np.meshgrid(theta, phi, indexing='ij')
psi_grid = np.stack([np.cos(TH / 2.0).ravel(),
                     (np.exp(1j * PH) * np.sin(TH / 2.0)).ravel()], axis=1)
F_opt_grid = single_copy_fidelities(V_sym, psi_grid)[0].reshape(TH.shape)
F_cnot_grid = single_copy_fidelities(V_cnot, psi_grid)[0].reshape(TH.shape)

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

# (1) convergence of the fixed-point iteration
ax1 = fig.add_subplot(2, 3, 1)
for k, h in enumerate(hist_sym):
    ax1.semilogy(np.maximum(F_OPT - h, 1e-16), lw=1.6, label=f"restart {k + 1}")
ax1.set_xlabel("iteration")
ax1.set_ylabel(r"$5/6-\bar F$")
ax1.set_title("Convergence of the polar-decomposition iteration")
ax1.legend(fontsize=8)
ax1.grid(alpha=0.3)

# (2) 3D: fidelity over the Bloch sphere (theta, phi)
ax2 = fig.add_subplot(2, 3, 2, projection='3d')
ax2.plot_surface(PH, TH, F_opt_grid, color='deepskyblue', alpha=0.6, linewidth=0)
ax2.plot_surface(PH, TH, F_cnot_grid, cmap='autumn', alpha=0.9, linewidth=0)
ax2.set_xlabel(r"$\varphi$")
ax2.set_ylabel(r"$\theta$")
ax2.set_zlabel(r"$F_A$")
ax2.set_zlim(0.4, 1.02)
ax2.set_title("Single-copy fidelity: optimal (flat, blue) vs CNOT (dome, warm)")
ax2.view_init(elev=24, azim=-58)

# (3) 3D: trade-off curve (F_A, F_B, lambda)
ax3 = fig.add_subplot(2, 3, 3, projection='3d')
ax3.plot(FA_curve, FB_curve, lam_grid, color='white', lw=1.2)
sc = ax3.scatter(FA_curve, FB_curve, lam_grid, c=lam_grid, cmap='plasma', s=45)
ax3.scatter([np.mean(FA_oct)], [np.mean(FB_oct)], [0.5], color='cyan', s=180, marker='*')
ax3.set_xlabel(r"$F_A$")
ax3.set_ylabel(r"$F_B$")
ax3.set_zlabel(r"weight $\lambda$")
ax3.set_title("Asymmetric cloning trade-off (star = symmetric optimum)")
ax3.view_init(elev=22, azim=-62)
fig.colorbar(sc, ax=ax3, shrink=0.6, pad=0.1, label=r"$\lambda$")

# (4) histogram over Haar-random inputs
ax4 = fig.add_subplot(2, 3, 4)
ax4.hist(FA_h, bins=80, range=(0.4, 1.0), color='deepskyblue', alpha=0.85, label="Optimal cloner")
ax4.hist(FA_cn_h, bins=80, range=(0.4, 1.0), color='orange', alpha=0.7, label="CNOT cloner")
ax4.set_yscale('log')
ax4.set_xlabel(r"single-copy fidelity $F_A$")
ax4.set_ylabel("count (20000 Haar-random inputs)")
ax4.set_title("Universality: fidelity distribution")
ax4.legend()
ax4.grid(alpha=0.3)

# (5) bar comparison of strategies
ax5 = fig.add_subplot(2, 3, 5)
labels = ["Measure &\nprepare\n(analytic)", "CNOT\ncloner\n(numerical)", "Optimized\ncloner\n(numerical)", "Bound 5/6\n(analytic)"]
vals = [2.0 / 3.0, float(np.mean(FA_cn_h)), float(np.mean(FA_h)), F_OPT]
bars = ax5.bar(labels, vals, color=['gray', 'orange', 'deepskyblue', 'limegreen'])
for b, v in zip(bars, vals):
    ax5.text(b.get_x() + b.get_width() / 2.0, v + 0.01, f"{v:.4f}", ha='center')
ax5.set_ylim(0.0, 1.0)
ax5.set_ylabel("average single-copy fidelity")
ax5.set_title("Strategy comparison")
ax5.grid(alpha=0.3, axis='y')

# (6) fidelities versus weight lambda
ax6 = fig.add_subplot(2, 3, 6)
ax6.plot(lam_grid, FA_curve, 'o-', color='deepskyblue', label=r"$F_A$")
ax6.plot(lam_grid, FB_curve, 's-', color='orange', label=r"$F_B$")
ax6.axhline(F_OPT, color='limegreen', ls='--', label="5/6")
ax6.axhline(0.5, color='gray', ls=':')
ax6.set_xlabel(r"weight $\lambda$")
ax6.set_ylabel("fidelity")
ax6.set_title("Fidelities of the optimized asymmetric cloners")
ax6.legend()
ax6.grid(alpha=0.3)

fig.suptitle("Optimal quantum cloning of a qubit: fidelity-maximizing channel", fontsize=16)
plt.tight_layout(rect=[0, 0, 1, 0.96])
plt.show()

5. Code Walkthrough

Section 1 – Dimensions

The output space is two qubits, so DIM_OUT = 4. The environment dimension is R_ENV = 8, which is the largest Choi rank a map from $\mathbb{C}^2$ to $\mathbb{C}^4$ can have. The isometry $V$ is therefore a $32\times 2$ complex matrix. Its row index encodes $(a,b,e)$ as $(2a+b)\cdot 8+e$, which matches the Kronecker-product ordering used later.

Section 2 – The octahedron design

octa_states lists the six states $|\pm x\rangle,|\pm y\rangle,|\pm z\rangle$. R_list stores their projectors $R_s$. QA and QB hold $P_s\otimes\mathbb{1}\otimes\mathbb{1}$ and $\mathbb{1}\otimes P_s\otimes\mathbb{1}$, so that $Q_s=\lambda Q^A_s+(1-\lambda)Q^B_s$ can be formed for any weight with a single line of array arithmetic.

Section 3 – Core routines

polar_isometry computes $UW^\dagger$ from a thin SVD. It is the projection of an arbitrary matrix onto the nearest isometry.
gradient_and_objective evaluates $G=\sum_s Q_sVR_s$ with one einsum call and returns $\bar F=\frac16\mathrm{Re},\mathrm{Tr}(V^\dagger G)$.
optimize_cloner runs the monotone iteration $V\leftarrow\mathrm{polar}(G(V))$ from several random isometries and keeps the best one. Random restarts guard against the possibility of non-global stationary points. The loop stops when the increase falls below $10^{-13}$.
single_copy_fidelities is a fully vectorized fidelity evaluator. It uses the identity
$$
F_A(\psi)=\sum_{b,e}\Big|\sum_a \langle\psi|a\rangle,\Psi_{abe}\Big|^2,\qquad \Psi=V|\psi\rangle,
$$
which avoids building density matrices and handles thousands of states in one call. This function runs the Haar-random verification and the Bloch-sphere grid with no Python loops.

Section 4 – The CNOT cloner

The CNOT cloner is encoded as an isometry with a trivial environment: $|0\rangle\to|00\rangle$ sits at row $0$ and $|1\rangle\to|11\rangle$ at row $3\cdot R_{\mathrm{ENV}}$. It serves as a non-universal baseline evaluated with exactly the same machinery.

Section 5 – The symmetric optimum

With $\lambda=\tfrac12$ the optimizer is run from six random starts. The script then reports the following:

the optimized average fidelity compared with $5/6$;
the isometry defect $|V^\dagger V-\mathbb{1}|$, which verifies that the channel is trace preserving;
the Bloch shrinking factor $\eta$ for the input $|0\rangle$, via $\rho_A=\mathrm{Tr}_{B,E}|\Psi\rangle\langle\Psi|$;
the numerical Kraus rank, taken from the singular values of the matrix whose rows are the vectorized Kraus operators;
a check against $20000$ Haar-random pure states, which confirms that the six-state average really equals the Haar average.

The theory predicts a Kraus rank of $2$, because the optimal cloner can be written as
$$
\Lambda(\rho)=\tfrac23,\Pi_{\mathrm{sym}}(\rho\otimes\mathbb{1})\Pi_{\mathrm{sym}},
$$
with $\Pi_{\mathrm{sym}}$ the projector onto the symmetric subspace of two qubits.

Section 6 – Asymmetric trade-off

The same optimizer is repeated for 21 weights $\lambda\in[0,1]$. At $\lambda=1$ the best strategy is to hand the input to $A$ untouched ($F_A=1$) and give $B$ only a maximally mixed qubit ($F_B=\tfrac12$). The intermediate weights trace out the Pareto frontier between the two copies.

Section 7 – Visualization

All six panels live in a single figure (see the next section for how to read it).

6. Execution Results

=== Symmetric optimal cloner (lambda = 0.5) ===
Optimized average fidelity      : 0.833333333333
Analytic optimum 5/6            : 0.833333333333
Absolute error                  : 6.052e-13
Isometry error ||V^dag V - I||  : 9.956e-16
F_A (octahedron average)        : 0.833333333333
F_B (octahedron average)        : 0.833333333333
Bloch shrinking factor eta      : 0.666664184649   (analytic 2/3 = 0.666666666667)
Numerical Kraus rank            : 2

=== Haar-random check with 20000 states ===
Optimal cloner : mean F_A = 0.833333322, mean F_B = 0.833333322, spread of F_A = 3.110e-06
CNOT cloner    : mean F_A = 0.667233172, min = 0.5000, max = 1.0000

=== Asymmetric trade-off (selected weights) ===
  lambda          F_A          F_B
    0.00   0.50000000   1.00000000
    0.25   0.60367259   0.98163706
    0.50   0.83333333   0.83333333
    0.75   0.98163706   0.60367259
    1.00   1.00000000   0.50000000

7. Reading the Results

Convergence (top left). The vertical axis is the gap $5/6-\bar F$ on a logarithmic scale. Because each update is a polar projection of the gradient of a convex quadratic form, every curve decreases monotonically. All restarts should approach the same floor near machine precision. That agreement is numerical evidence that the symmetric optimum is a global one and that the algorithm does not stall in spurious local maxima.

Fidelity over the Bloch sphere (top middle, 3D). The horizontal plane is the parameter space $(\varphi,\theta)$ of input states, and the height is $F_A$. The optimized cloner appears as a perfectly flat sheet at $5/6$: it copies every state equally well, which is exactly what universal means. The CNOT cloner forms a dome, $1-\tfrac12\sin^2\theta$. It is perfect at the poles ($|0\rangle$ and $|1\rangle$, the basis states it was designed for) and drops to $\tfrac12$ on the equator. Its average over the sphere is only $\tfrac23$. The flatness of the blue sheet is a numerical rediscovery of the covariance of the optimal channel, which the optimizer was never told to impose.

Trade-off curve (top right, 3D). Each point is one optimized cloner with $x=F_A$, $y=F_B$ and height equal to the weight $\lambda$. The curve runs from $(F_A,F_B)=(\tfrac12,1)$ at $\lambda=0$ to $(1,\tfrac12)$ at $\lambda=1$, and passes through the symmetric optimum $(\tfrac56,\tfrac56)$ marked by the star. Its concavity matters physically: improving one copy beyond $5/6$ always costs more than a one-for-one loss in the other, which is a quantitative form of the no-cloning theorem.

Fidelity histogram (bottom left). For 20000 Haar-random inputs, the optimal cloner puts essentially all of its weight in a single bin at $5/6$ (note the logarithmic vertical axis). The CNOT cloner is spread over $[\tfrac12,1]$. The spread, not just the mean, separates a universal machine from a basis-dependent one.

Strategy comparison (bottom middle). Measure-and-prepare and the CNOT cloner both give $\tfrac23$. The optimized quantum cloner reaches $\tfrac56$, which matches the analytic bound. Quantum-coherent copying beats the best classical strategy of measuring once and re-preparing.

Fidelities versus the weight (bottom right). $F_A$ rises and $F_B$ falls as $\lambda$ increases, and the curves cross at $\lambda=\tfrac12$ at the height $5/6$. They saturate at the extremes: $F_A=1$ means the second copy carries no information beyond a random guess of fidelity $\tfrac12$.

8. Summary

Optimal approximate cloning is a convex quadratic maximization over isometries, which is a semidefinite program in disguise.
Replacing the Haar integral with an exact 2-design (six octahedron states) reduces the problem to a few small matrix products.
The fixed-point iteration $V\leftarrow\mathrm{polar}\big(\sum_s Q_sVR_s\big)$ increases the objective monotonically and recovers the Bužek–Hillery value $\bar F=5/6$, the shrinking factor $\eta=2/3$ and a Kraus rank of $2$, without any hand-built ansatz.
The same code, with only the weight $\lambda$ changed, maps the whole asymmetric trade-off curve between the two copies.

Designing the Optimal Quantum POVM

October 4, 2026

Maximizing Mutual Information and Accessible Information with Python

Suppose Alice sends Bob one of several non-orthogonal quantum states. Bob wants to learn as much as possible about which one she sent. Which measurement should he perform?

This is the problem of accessible information. In this article we solve it numerically for a classic example, the trine ensemble, and check the result against the known analytic answer. Along the way we look at the geometry of the optimal measurement on the Bloch sphere.

1. Problem Setting

An ensemble is a set of states $\rho_i$ with prior probabilities $p_i$:

$$
\mathcal{E}={p_i,\rho_i}_{i=1}^{N}.
$$

Bob measures with a POVM ${E_k}$, where

$$
E_k\ge 0,\qquad \sum_k E_k=\mathbb{1}.
$$

The conditional probability of outcome $k$ given state $i$ is

$$
p(k|i)=\mathrm{Tr}(\rho_i E_k),
$$

and the outcome distribution is

$$
q_k=\sum_i p_i,p(k|i).
$$

The mutual information between the sent index $i$ and the outcome $k$ is

$$
I(X!:!Y)=\sum_{i,k}p_i,p(k|i)\log_2\frac{p(k|i)}{q_k}.
$$

The accessible information is the maximum over all POVMs:

$$
I_{\mathrm{acc}}=\max_I(X!:!Y).
$$

It is bounded above by the Holevo quantity:

$$
I_{\mathrm{acc}}\le\chi=S!\left(\sum_i p_i\rho_i\right)-\sum_i p_i S(\rho_i),
$$

where $S(\rho)=-\mathrm{Tr},\rho\log_2\rho$ is the von Neumann entropy.

Two facts shape the numerical strategy:

By Davies’ theorem, an optimal POVM can always be chosen with rank-one elements, and at most $d^2$ of them. For a qubit ($d=2$), four outcomes are enough.
The objective is non-convex in the POVM, so a single local optimization can get stuck. We use many random restarts.

2. The Example: the Trine Ensemble

The trine ensemble consists of three pure qubit states whose Bloch vectors lie in the $x$–$z$ plane, $120^\circ$ apart, each with probability $1/3$:

$$
\vec r_i=\bigl(\sin\phi_i,;0,;\cos\phi_i\bigr),\qquad \phi_i=\frac{2\pi i}{3},\quad i=0,1,2,
$$

$$
\rho_i=\frac12\bigl(\mathbb{1}+\vec r_i\cdot\vec\sigma\bigr).
$$

Because $\sum_i\vec r_i=0$, the average state is maximally mixed, $\bar\rho=\mathbb{1}/2$. Since the states are pure,

$$
\chi=S(\bar\rho)-0=1\ \text{bit}.
$$

Three non-orthogonal states cannot be perfectly distinguished, so the accessible information must be strictly smaller than $1$ bit.

A closed-form candidate

Consider the anti-trine POVM

$$
E_k=\frac23\cdot\frac12\bigl(\mathbb{1}-\vec r_k\cdot\vec\sigma\bigr),\qquad k=0,1,2.
$$

Each element projects onto the state orthogonal to $\rho_k$, and $\sum_k E_k=\mathbb{1}$ because $\sum_k\vec r_k=0$. The conditional probabilities are

$$
p(k|i)=\frac13\bigl(1-\vec r_i\cdot\vec r_k\bigr)=
\begin{cases}
0 & k=i,\[2pt]
\tfrac12 & k\neq i,
\end{cases}
$$

so $q_k=1/3$ and

$$
I=H(Y)-H(Y|X)=\log_2 3-1=\log_2\frac32\approx 0.585\ \text{bits}.
$$

This is the known value of the accessible information of the trine. In the rest of the article we let the computer rediscover it, without telling the optimizer anything about the structure.

3. Strategy for the Numerical Optimization

Parametrizing a valid POVM

Optimizers work best without constraints, so we build a POVM that is valid by construction. Take $K$ arbitrary complex vectors $|v_k\rangle\in\mathbb{C}^2$, form

$$
S=\sum_k|v_k\rangle\langle v_k|,
$$

and define

$$
E_k=S^{-1/2},|v_k\rangle\langle v_k|,S^{-1/2}.
$$

Then $E_k\ge0$ and $\sum_kE_k=S^{-1/2}SS^{-1/2}=\mathbb{1}$ automatically. The $4K$ real numbers (real and imaginary parts of the vectors) can be optimized freely with BFGS.

Bloch representation of a POVM element

For a qubit, every rank-one element can be written

$$
E_k=a_k\bigl(\mathbb{1}+\vec m_k\cdot\vec\sigma\bigr),\qquad |\vec m_k|=1,
$$

with the completeness conditions

$$
\sum_k a_k=1,\qquad \sum_k a_k\vec m_k=\vec 0.
$$

The pair $(a_k,\vec m_k)$ is easy to visualize: $\vec m_k$ is a direction on the Bloch sphere and $a_k$ is its weight.

A two-parameter family for a 3D landscape

To see the optimum as a peak, we define a rotated and “softened” anti-trine family:

$$
E_k(\alpha,\eta)=\frac23\cdot\frac12\bigl(\mathbb{1}+\eta,R_y(\alpha)(-\vec r_k)\cdot\vec\sigma\bigr),\qquad 0\le\eta\le1.
$$

Here $\alpha$ rotates the measurement within the $x$–$z$ plane, and $\eta$ controls the sharpness ($\eta=0$ is a trivial measurement, $\eta=1$ is rank-one). Every member is a valid POVM, so $I(\alpha,\eta)$ is a legitimate surface to plot.

4. The Complete Source Code

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

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

# ------------------------------------------------------------
# 1. Pauli matrices and the trine ensemble
# ------------------------------------------------------------
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)
PAULI = np.array([SX, SY, SZ])

def bloch_to_rho(r):
    return 0.5 * (I2 + np.tensordot(r, PAULI, axes=1))

angles = 2 * np.pi * np.arange(3) / 3
bloch_states = np.stack([np.sin(angles), np.zeros(3), np.cos(angles)], axis=1)
priors = np.ones(3) / 3
rhos = np.array([bloch_to_rho(r) for r in bloch_states])

def entropy_bits(rho):
    w = np.linalg.eigvalsh(rho)
    w = w[w > 1e-12]
    return float(-np.sum(w * np.log2(w)))

rho_avg = np.tensordot(priors, rhos, axes=1)
holevo = entropy_bits(rho_avg) - sum(p * entropy_bits(r) for p, r in zip(priors, rhos))

# ------------------------------------------------------------
# 2. Mutual information of a POVM
# ------------------------------------------------------------
def mutual_information(povm, rhos, priors):
    povm = np.asarray(povm)
    cond = np.real(np.einsum("kab,iba->ik", povm, rhos))
    cond = np.clip(cond, 0.0, None)
    pk = priors @ cond
    ratio = np.where(cond > 1e-15, cond / np.maximum(pk[None, :], 1e-15), 1.0)
    return float(np.sum(priors[:, None] * cond * np.log2(ratio)))

# ------------------------------------------------------------
# 3. Unconstrained parametrization of a K-outcome POVM
# ------------------------------------------------------------
def params_to_povm(x, K):
    x = x.reshape(K, 4)
    v = x[:, 0:2] + 1j * x[:, 2:4]
    G = np.einsum("ka,kb->kab", v, v.conj())
    S = G.sum(axis=0)
    w, U = np.linalg.eigh(S)
    w = np.maximum(w, 1e-12)
    S_inv_sqrt = (U * (1.0 / np.sqrt(w))) @ U.conj().T
    return np.einsum("ab,kbc,cd->kad", S_inv_sqrt, G, S_inv_sqrt)

def optimize_povm(K, n_restarts=20):
    best, values = None, []
    for _ in range(n_restarts):
        x0 = rng.normal(size=4 * K)
        history = [mutual_information(params_to_povm(x0, K), rhos, priors)]
        def objective(x):
            return -mutual_information(params_to_povm(x, K), rhos, priors)
        def callback(xk):
            history.append(-objective(xk))
        res = minimize(objective, x0, method="BFGS", callback=callback,
                       options={"gtol": 1e-10, "maxiter": 500})
        values.append(-res.fun)
        if best is None or -res.fun > best["value"]:
            best = {"value": -res.fun, "x": res.x, "history": history}
    return best, np.array(values)

best3, vals3 = optimize_povm(3)
best4, vals4 = optimize_povm(4)
povm4 = params_to_povm(best4["x"], 4)

# ------------------------------------------------------------
# 4. Merge outcomes that share the same Bloch direction
# ------------------------------------------------------------
def povm_to_bloch(povm):
    weights, vecs = [], []
    for E in povm:
        a = np.real(np.trace(E)) / 2
        if a > 1e-9:
            m = np.real(np.array([np.trace(E @ s) for s in PAULI])) / (2 * a)
        else:
            m = np.zeros(3)
        weights.append(a)
        vecs.append(m)
    return np.array(weights), np.array(vecs)

def merge_povm(povm, tol=1e-3):
    w, m = povm_to_bloch(povm)
    groups = []
    for k in range(len(w)):
        if w[k] < 1e-6:
            continue
        for g in groups:
            if np.linalg.norm(g["m"] - m[k]) < tol * 10:
                tot = g["w"] + w[k]
                g["m"] = (g["m"] * g["w"] + m[k] * w[k]) / tot
                g["w"] = tot
                break
        else:
            groups.append({"w": w[k], "m": m[k].copy()})
    return np.array([g["w"] for g in groups]), np.array([g["m"] for g in groups])

w_merged, m_merged = merge_povm(povm4)
povm_merged = [w * bloch_to_rho(m) * 2 for w, m in zip(w_merged, m_merged)]
I_merged = mutual_information(povm_merged, rhos, priors)

# ------------------------------------------------------------
# 5. Reference measurements
# ------------------------------------------------------------
anti_dirs = -bloch_states
anti_trine = [(2.0 / 3.0) * bloch_to_rho(n) for n in anti_dirs]
trine_povm = [(2.0 / 3.0) * bloch_to_rho(n) for n in bloch_states]
I_anti = mutual_information(anti_trine, rhos, priors)
I_trine = mutual_information(trine_povm, rhos, priors)

theta = np.linspace(0, np.pi, 61)
phi = np.linspace(0, 2 * np.pi, 121)
TH, PH = np.meshgrid(theta, phi, indexing="ij")
NX, NY, NZ = np.sin(TH) * np.cos(PH), np.sin(TH) * np.sin(PH), np.cos(TH)
dots = np.stack([NX * r[0] + NY * r[1] + NZ * r[2] for r in bloch_states], axis=0)
p_plus = np.clip(0.5 * (1 + dots), 1e-15, 1)
p_minus = np.clip(0.5 * (1 - dots), 1e-15, 1)
q_plus = np.tensordot(priors, p_plus, axes=1)
q_minus = np.tensordot(priors, p_minus, axes=1)
I_map = np.zeros_like(TH)
for i in range(3):
    I_map += priors[i] * (p_plus[i] * np.log2(p_plus[i] / q_plus)
                          + p_minus[i] * np.log2(p_minus[i] / q_minus))
idx = np.unravel_index(np.argmax(I_map), I_map.shape)
I_proj_best = I_map[idx]
n_proj_best = np.array([NX[idx], NY[idx], NZ[idx]])

# Rotation angle alpha and shrink factor eta family
def rotated_povm(alpha, eta):
    c, s = np.cos(alpha), np.sin(alpha)
    out = []
    for n in anti_dirs:
        nr = np.array([c * n[0] + s * n[2], n[1], -s * n[0] + c * n[2]])
        out.append((2.0 / 3.0) * bloch_to_rho(eta * nr))
    return out

alphas = np.linspace(0, 2 * np.pi, 121)
etas = np.linspace(0, 1, 41)
I_surf = np.array([[mutual_information(rotated_povm(a, e), rhos, priors)
                    for a in alphas] for e in etas])
A, E = np.meshgrid(alphas, etas)

# Random POVMs as a baseline
n_random = 2000
I_random = np.array([mutual_information(params_to_povm(rng.normal(size=12), 3), rhos, priors)
                     for _ in range(n_random)])

# ------------------------------------------------------------
# 6. Console output
# ------------------------------------------------------------
print("=== Ensemble: trine states (equal priors) ===")
print(f"Holevo bound chi                    : {holevo:.6f} bits")
print(f"Theoretical optimum log2(3/2)       : {np.log2(1.5):.6f} bits")
print()
print("=== Numerical optimization ===")
print(f"3-outcome POVM (best of {len(vals3)})        : {best3['value']:.6f} bits")
print(f"4-outcome POVM (best of {len(vals4)})        : {best4['value']:.6f} bits")
print(f"4-outcome POVM after merging         : {I_merged:.6f} bits ({len(w_merged)} distinct outcomes)")
print(f"Spread over restarts (3 outcomes)    : min {vals3.min():.6f} / max {vals3.max():.6f}")
print()
print("=== Reference measurements ===")
print(f"Best projective measurement          : {I_proj_best:.6f} bits")
print(f"Trine-aligned POVM                   : {I_trine:.6f} bits")
print(f"Anti-trine POVM                      : {I_anti:.6f} bits")
print(f"Best of {n_random} random 3-outcome POVMs  : {I_random.max():.6f} bits")
print()
print("=== Optimal POVM (merged): weights and Bloch directions ===")
for k, (w, m) in enumerate(zip(w_merged, m_merged)):
    print(f"E_{k}: weight a = {w:.4f},  direction m = {np.round(m, 4)}")
print(f"Sum of weights = {w_merged.sum():.6f}")
print()
print("=== Completeness check ===")
print("Sum of E_k =")
print(np.round(sum(povm_merged), 6))
print()
print("=== Conditional probabilities p(k|i) of the optimal POVM ===")
cond = np.real(np.einsum("kab,iba->ik", np.array(povm_merged), rhos))
print(np.round(cond, 4))

# ------------------------------------------------------------
# 7. Figure (all panels in one output)
# ------------------------------------------------------------
fig = plt.figure(figsize=(20, 12))
fig.suptitle("Information-maximizing POVM for the trine ensemble", fontsize=18, y=0.99)

sphere_u = np.linspace(0, 2 * np.pi, 40)
sphere_v = np.linspace(0, np.pi, 20)
SPX = np.outer(np.cos(sphere_u), np.sin(sphere_v))
SPY = np.outer(np.sin(sphere_u), np.sin(sphere_v))
SPZ = np.outer(np.ones_like(sphere_u), np.cos(sphere_v))

# (a) mutual information of projective measurements over the Bloch sphere
ax1 = fig.add_subplot(2, 3, 1, projection="3d")
norm = colors.Normalize(vmin=I_map.min(), vmax=I_map.max())
ax1.plot_surface(NX, NY, NZ, facecolors=cm.viridis(norm(I_map)), rstride=1, cstride=1,
                 linewidth=0, antialiased=False, shade=False)
for r in bloch_states:
    ax1.quiver(0, 0, 0, *(1.25 * r), color="white", linewidth=3, arrow_length_ratio=0.08)
    ax1.quiver(0, 0, 0, *(1.25 * r), color="tab:blue", linewidth=1.5, arrow_length_ratio=0.08)
ax1.scatter(*n_proj_best, color="red", s=60)
ax1.set_box_aspect((1, 1, 1))
ax1.set_title("(a) I of projective measurement vs. axis n")
ax1.set_xlabel("x"); ax1.set_ylabel("y"); ax1.set_zlabel("z")
mappable = cm.ScalarMappable(norm=norm, cmap="viridis")
mappable.set_array([])
fig.colorbar(mappable, ax=ax1, shrink=0.55, pad=0.1, label="I (bits)")

# (b) signal states and optimal POVM on the Bloch sphere
ax2 = fig.add_subplot(2, 3, 2, projection="3d")
ax2.plot_wireframe(SPX, SPY, SPZ, color="lightgray", linewidth=0.4)
for k, r in enumerate(bloch_states):
    ax2.quiver(0, 0, 0, *r, color="tab:blue", linewidth=2.5, arrow_length_ratio=0.1,
               label="signal states" if k == 0 else None)
for k, (w, m) in enumerate(zip(w_merged, m_merged)):
    ax2.quiver(0, 0, 0, *m, color="tab:red", linewidth=2.5, arrow_length_ratio=0.1,
               label="optimal POVM directions" if k == 0 else None)
ax2.set_box_aspect((1, 1, 1))
ax2.set_title("(b) Signal states and optimal POVM")
ax2.set_xlabel("x"); ax2.set_ylabel("y"); ax2.set_zlabel("z")
ax2.legend(loc="upper left", fontsize=8)

# (c) 3D surface of I(alpha, eta)
ax3 = fig.add_subplot(2, 3, 3, projection="3d")
surf = ax3.plot_surface(A, E, I_surf, cmap="plasma", linewidth=0, antialiased=True)
ax3.set_xlabel(r"rotation $\alpha$ (rad)")
ax3.set_ylabel(r"sharpness $\eta$")
ax3.set_zlabel("I (bits)")
ax3.set_title(r"(c) $I(\alpha,\eta)$ for rotated, shrunk anti-trine POVMs")
ax3.view_init(elev=28, azim=-60)
fig.colorbar(surf, ax=ax3, shrink=0.55, pad=0.1, label="I (bits)")

# (d) bar chart comparison
ax4 = fig.add_subplot(2, 3, 4)
labels = ["Holevo\nbound", "Optimal\nPOVM", "Anti-trine", "Best\nprojective", "Trine-\naligned"]
values = [holevo, best3["value"], I_anti, I_proj_best, I_trine]
bar_colors = ["#999999", "#d62728", "#ff7f0e", "#1f77b4", "#2ca02c"]
bars = ax4.bar(labels, values, color=bar_colors)
for b, v in zip(bars, values):
    ax4.text(b.get_x() + b.get_width() / 2, v + 0.015, f"{v:.3f}", ha="center", fontsize=10)
ax4.set_ylabel("bits")
ax4.set_ylim(0, 1.15)
ax4.set_title("(d) Information gain by measurement")
ax4.grid(axis="y", alpha=0.3)

# (e) optimization trajectories
ax5 = fig.add_subplot(2, 3, 5)
ax5.plot(best3["history"], marker="o", label="3 outcomes (best run)")
ax5.plot(best4["history"], marker="s", label="4 outcomes (best run)")
ax5.axhline(np.log2(1.5), color="k", linestyle="--", label=r"$\log_2(3/2)$")
ax5.set_xlabel("BFGS iteration")
ax5.set_ylabel("I (bits)")
ax5.set_title("(e) Convergence of the optimization")
ax5.legend()
ax5.grid(alpha=0.3)

# (f) random POVM baseline
ax6 = fig.add_subplot(2, 3, 6)
ax6.hist(I_random, bins=40, color="#8c8cff", edgecolor="white")
ax6.axvline(np.log2(1.5), color="red", linewidth=2, label=r"optimum $\log_2(3/2)$")
ax6.axvline(I_random.mean(), color="k", linestyle=":", label="random mean")
ax6.set_xlabel("I (bits)")
ax6.set_ylabel("count")
ax6.set_title(f"(f) {n_random} random 3-outcome POVMs")
ax6.legend()
ax6.grid(alpha=0.3)

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

5. Code Walkthrough

Section 1: Ensemble construction

bloch_to_rho converts a Bloch vector $\vec r$ into $\rho=\tfrac12(\mathbb 1+\vec r\cdot\vec\sigma)$ with a single tensordot. The three trine vectors are generated from the angles $0,,2\pi/3,,4\pi/3$, and rhos is an array of shape (3, 2, 2). entropy_bits computes the von Neumann entropy from the eigenvalues, discarding values below $10^{-12}$ so that $0\log 0$ never produces nan. The Holevo quantity $\chi$ is computed directly from its definition and serves as the upper bound.

Section 2: Mutual information

mutual_information is the heart of the program. The line

1	cond = np.real(np.einsum("kab,iba->ik", povm, rhos))

computes all conditional probabilities $p(k|i)=\mathrm{Tr}(\rho_iE_k)$ in one shot. The index pattern kab,iba->ik contracts $E_k^{ab}\rho_i^{ba}$, which is exactly the trace of a matrix product. The remaining lines compute $q_k$ and the sum $\sum p_i,p(k|i)\log_2\frac{p(k|i)}{q_k}$. The np.where guard enforces the convention $0\log 0=0$, so zero-probability entries (which occur at the optimum!) do not generate warnings or nan.

Section 3: POVM parametrization and optimization

params_to_povm implements $E_k=S^{-1/2}|v_k\rangle\langle v_k|S^{-1/2}$. The matrix $S^{-1/2}$ is built from the eigendecomposition of the Hermitian matrix $S$, with eigenvalues floored at $10^{-12}$ for numerical safety. Because completeness holds by construction, the optimizer needs no constraints and BFGS can be used directly.

optimize_povm runs BFGS from n_restarts random starting points, records the objective after every iteration through the callback, and keeps the best run. Restarts matter here: the landscape has stationary points (for example one corresponding to projective measurements), and some runs may settle there instead of at the global optimum.

We run it for $K=3$ and $K=4$ outcomes. Since Davies’ theorem says four elements suffice for a qubit, the $K=4$ run is the general case, and $K=3$ confirms that the structure of the optimum is simpler than the upper limit suggests.

Section 4: Merging duplicate outcomes

When $K=4$, the optimizer may split one physical outcome into two elements pointing in the same direction. Merging them does not change the mutual information (coarse-graining outcomes with identical likelihood ratios loses nothing). povm_to_bloch extracts the weight $a_k=\mathrm{Tr}E_k/2$ and the direction $\vec m_k$ from each element, and merge_povm groups elements whose directions agree within a tolerance and sums their weights. The merged POVM is rebuilt and its information is recomputed as a consistency check.

Section 5: Reference measurements and the landscapes

Three baselines are computed:

Anti-trine POVM: the analytic candidate from Section 2 above.
Trine-aligned POVM: the same construction with $+\vec r_k$ instead of $-\vec r_k$, a natural but inferior guess.
Best projective measurement: for a two-outcome projective measurement along axis $\hat n$, the likelihoods are $p(\pm|i)=\tfrac12(1\pm\hat n\cdot\vec r_i)$. This is evaluated on a $61\times121$ grid over the whole sphere.

Speed-up note. The grid evaluation is fully vectorized with NumPy broadcasting; the only Python-level loop runs over the three states, not over the 7,381 grid points. The POVM objective uses einsum instead of nested loops, and the unconstrained parametrization lets us use a fast quasi-Newton method instead of a constrained solver. As a result, the whole program, including 40 optimization runs, a $41\times121$ surface, and 2,000 random POVMs, finishes within a few seconds to tens of seconds.

The family rotated_povm(alpha, eta) rotates the anti-trine directions by $\alpha$ in the $x$–$z$ plane and shrinks them by $\eta$, which gives the surface $I(\alpha,\eta)$. Finally, 2,000 random 3-outcome POVMs provide a “no optimization” baseline.

Section 6: Console output

All numerical results are printed: the Holevo bound, the theoretical value $\log_2(3/2)$, the optimizer results, the reference values, the merged POVM in Bloch form, a completeness check $\sum_kE_k=\mathbb 1$, and the matrix of conditional probabilities $p(k|i)$.

Section 7: The figure

All six panels are drawn into a single fig and displayed once with plt.show(). Panels (a), (b), and (c) are 3D.

6. Execution Results

Figure output

Console output

=== Ensemble: trine states (equal priors) ===
Holevo bound chi                    : 1.000000 bits
Theoretical optimum log2(3/2)       : 0.584963 bits

=== Numerical optimization ===
3-outcome POVM (best of 20)        : 0.584963 bits
4-outcome POVM (best of 20)        : 0.584963 bits
4-outcome POVM after merging         : 0.584963 bits (3 distinct outcomes)
Spread over restarts (3 outcomes)    : min 0.459148 / max 0.584963

=== Reference measurements ===
Best projective measurement          : 0.459148 bits
Trine-aligned POVM                   : 0.333333 bits
Anti-trine POVM                      : 0.584963 bits
Best of 2000 random 3-outcome POVMs  : 0.532563 bits

=== Optimal POVM (merged): weights and Bloch directions ===
E_0: weight a = 0.3333,  direction m = [-0. -0. -1.]
E_1: weight a = 0.3333,  direction m = [-0.866  0.     0.5  ]
E_2: weight a = 0.3333,  direction m = [0.866 0.    0.5  ]
Sum of weights = 1.000000

=== Completeness check ===
Sum of E_k =
[[1.+0.j 0.+0.j]
 [0.-0.j 1.+0.j]]

=== Conditional probabilities p(k|i) of the optimal POVM ===
[[0.  0.5 0.5]
 [0.5 0.  0.5]
 [0.5 0.5 0. ]]

7. Interpreting the Results

Console output

The optimizer should report

$$
I_{\mathrm{acc}}\approx 0.584963\ \text{bits}=\log_2\frac32,
$$

for both the 3-outcome and the 4-outcome searches. The numerical optimum therefore reproduces the analytic value, even though the optimizer started from random vectors and knew nothing about the anti-trine structure.

Several facts stand out:

The gap to the Holevo bound is large. The Holevo quantity is $\chi=1$ bit, but only about $0.585$ bits are accessible. The Holevo bound is not achievable for this ensemble, because the states are non-orthogonal and no measurement can fully resolve them.
The 4-outcome optimum collapses to 3 outcomes. After merging, only three distinct directions remain. Allowing the maximum number of outcomes allowed by Davies’ theorem brings no extra information here.
The optimal conditional probabilities are $0$ and $\tfrac12$. Each outcome rules out exactly one state: outcome $k$ never occurs when state $k$ was sent, and otherwise occurs with probability $1/2$. This matrix is $p(k|i)=\tfrac12(1-\delta_{ik})$ up to a permutation of outcomes.
The optimal weights are $a_k=1/3$ with directions equal to the anti-trine vectors $-\vec r_k$ (in some order), so $\sum_ka_k=1$ and $\sum_ka_k\vec m_k=\vec0$.
Restarts matter. The spread over restarts shows that some runs end below the optimum, typically at the plateau of the best projective measurement (about $0.459$ bits). Without restarts, a single unlucky run could have been mistaken for the answer.

Panel (a): projective measurements over the Bloch sphere

The color shows the mutual information obtained by a two-outcome projective measurement along each axis $\hat n$. The map is symmetric under $\hat n\to-\hat n$ and has a pattern with the three-fold symmetry of the trine. The bright regions lie in the $x$–$z$ plane near the signal directions, and the dark regions lie near the $\pm y$ poles: measuring along $y$ is orthogonal to every signal Bloch vector, so it returns pure noise and $I=0$. The maximum over all projective measurements, marked by the red dot, is only about $0.459$ bits, clearly below $\log_2(3/2)$. Projective measurements are not enough: a genuinely generalized measurement beats the best von Neumann measurement.

Panel (b): signal states and the optimal POVM

The blue arrows are the three trine states and the red arrows are the directions of the optimal POVM elements. All six vectors lie in the $x$–$z$ plane, and each red arrow points opposite to one blue arrow. Geometrically, the optimal measurement asks, for each state, the question “is the state not this one?”. The three red vectors form an equilateral triangle, so $\sum_ka_k\vec m_k=0$ is satisfied and the elements add up to the identity.

Panel (c): the 3D landscape $I(\alpha,\eta)$

This surface is the most informative picture of the optimization problem.

At $\eta=0$ the surface is flat at $I=0$: a POVM with $E_k=\mathbb 1/3$ ignores the state entirely.
Along $\eta$, information grows monotonically; sharper measurements always help in this family.
Along $\alpha$, the surface oscillates with period $2\pi/3$, reflecting the three-fold symmetry of the ensemble.
At $\eta=1$ the peaks (value $\log_2(3/2)\approx0.585$) sit at $\alpha=0,,2\pi/3,,4\pi/3$, where the measurement is exactly the anti-trine.
Midway between the peaks, at $\alpha=\pi/3$ and its translates, the rotated anti-trine coincides with the trine-aligned POVM. There $I=1/3$ bit, the bottom of the ridge. The same physical hardware, rotated by $60^\circ$ in the Bloch sphere, loses almost half of the information.

So a sharp, correctly oriented measurement is a global maximum within this family, and the landscape shows how quickly the information degrades when the measurement is misaligned.

Panel (d): bar chart

The five bars give the ranking: Holevo bound ($1$) $>$ optimal POVM $=$ anti-trine ($\approx0.585$) $>$ best projective ($\approx0.459$) $>$ trine-aligned ($1/3$). The optimal POVM and the anti-trine bars coincide, confirming the analytic solution. The distance between the gray Holevo bar and the red optimum bar is the part of the information that is encoded in the ensemble but cannot be extracted by any measurement.

Panel (e): convergence

Both curves start from a low value determined by the random initial vectors and climb to the dashed line at $\log_2(3/2)$ within a few dozen BFGS iterations. The 4-outcome run converges to the same level as the 3-outcome run, matching the merging result: the extra degree of freedom is not used.

Panel (f): random POVMs

The histogram of 2,000 random 3-outcome POVMs lies entirely to the left of the red line. Random POVMs typically capture a fraction of the available information, and even the best of 2,000 draws stays below the optimum. This is a sanity check that the optimizer is doing real work, and that the value $\log_2(3/2)$ is a ceiling and not just one more point in the distribution.

8. Conclusion

We computed the accessible information of the trine ensemble by optimizing over POVMs, using a parametrization that enforces completeness automatically. The optimum is

$$
I_{\mathrm{acc}}=\log_2\frac32\approx0.585\ \text{bits},
$$

attained by the three-outcome anti-trine POVM

$$
E_k=\frac23\cdot\frac12\bigl(\mathbb 1-\vec r_k\cdot\vec\sigma\bigr),
$$

while the best projective measurement gives only about $0.459$ bits and the Holevo bound of $1$ bit is out of reach. The same pipeline (rank-one POVM parametrization, multi-start BFGS, and Bloch-sphere visualization) works for any qubit ensemble: change bloch_states and priors and the program finds the information-maximizing measurement for you. For higher dimensions, the parametrization generalizes directly by using $d$-dimensional complex vectors and up to $d^2$ outcomes.

Unambiguous Quantum State Discrimination

October 3, 2026

Maximizing the Error-Free Success Probability

Introduction

Suppose someone hands you a qubit and tells you it is in one of two known pure states, $|\psi_1\rangle$ or $|\psi_2\rangle$. If the states are not orthogonal, no measurement can tell them apart perfectly. Two philosophies exist for dealing with this limitation:

Minimum-error discrimination always gives an answer, but accepts that the answer is sometimes wrong. The optimum is given by the Helstrom bound.
Unambiguous state discrimination (USD) never gives a wrong answer, but accepts that the measurement sometimes returns “I don’t know”. The goal is to maximize the probability of a conclusive and correct result.

USD is the right tool whenever a wrong answer is far more costly than no answer at all, as in quantum key distribution or quantum money verification. In this article we formulate USD as a small constrained optimization problem, derive the closed-form optimum, verify it by brute-force search, simulate the measurement with Monte Carlo sampling, and visualize everything, including 3D landscapes.

Problem Formulation

Let the two states be prepared with prior probabilities $\eta_1$ and $\eta_2 = 1-\eta_1$, and let their overlap be

$$
s = \langle \psi_1 | \psi_2 \rangle, \qquad 0 \le |s| < 1 .
$$

We look for a three-outcome POVM ${\Pi_1, \Pi_2, \Pi_?}$ with

$$
\Pi_1 + \Pi_2 + \Pi_? = I, \qquad \Pi_1,\Pi_2,\Pi_? \succeq 0 .
$$

Outcome “1” means “the state was $\psi_1$”, outcome “2” means “the state was $\psi_2$”, and outcome “?” means “inconclusive”. Error-free discrimination requires

$$
\langle \psi_2 | \Pi_1 | \psi_2 \rangle = 0, \qquad \langle \psi_1 | \Pi_2 | \psi_1 \rangle = 0 .
$$

The objective is the total probability of a conclusive result:

$$
P_{\mathrm{succ}} = \eta_1 \langle \psi_1 | \Pi_1 | \psi_1 \rangle + \eta_2 \langle \psi_2 | \Pi_2 | \psi_2 \rangle .
$$

The problem is therefore: maximize $P_{\mathrm{succ}}$ over POVMs subject to the zero-error constraints.

Reduction to a Two-Variable Optimization

The zero-error constraints force $\Pi_1$ to be proportional to the projector onto the vector orthogonal to $|\psi_2\rangle$, and $\Pi_2$ to be proportional to the projector onto the vector orthogonal to $|\psi_1\rangle$:

$$
\Pi_1 = p_1 |\psi_2^{\perp}\rangle\langle\psi_2^{\perp}|, \qquad
\Pi_2 = p_2 |\psi_1^{\perp}\rangle\langle\psi_1^{\perp}|, \qquad p_1,p_2 \in [0,1].
$$

Since $|\langle \psi_2^{\perp}|\psi_1\rangle|^2 = 1-|s|^2$, the objective becomes

$$
P_{\mathrm{succ}}(p_1,p_2) = (1-|s|^2),(\eta_1 p_1 + \eta_2 p_2).
$$

The remaining requirement is positivity of the inconclusive operator

$$
\Pi_? = I - p_1 |\psi_2^{\perp}\rangle\langle\psi_2^{\perp}| - p_2 |\psi_1^{\perp}\rangle\langle\psi_1^{\perp}| \succeq 0 .
$$

For a $2\times 2$ Hermitian matrix, positive semidefiniteness is equivalent to a non-negative trace and a non-negative determinant:

$$
\operatorname{tr}\Pi_? = 2 - p_1 - p_2 \ge 0, \qquad
\det \Pi_? = 1 - p_1 - p_2 + p_1 p_2 ,(1-|s|^2) \ge 0 .
$$

We have reduced the quantum problem to maximizing a linear function over a convex region in the $(p_1,p_2)$ plane.

Closed-Form Solution

Let $r = \sqrt{\eta_2/\eta_1}$. The optimum lies on the curve $\det\Pi_? = 0$, and solving the Lagrange conditions gives two regimes.

Regime I (moderate overlap), $|s| \le \min\left(r,, 1/r\right)$:

$$
p_1 = \frac{1 - r|s|}{1-|s|^2}, \qquad
p_2 = \frac{1 - |s|/r}{1-|s|^2}, \qquad
P_{\mathrm{succ}} = 1 - 2\sqrt{\eta_1\eta_2},|s| .
$$

Regime II (large overlap, strongly biased priors), $|s| > \min\left(r,,1/r\right)$. It is best to only identify the more probable state:

$$
P_{\mathrm{succ}} = \max(\eta_1,\eta_2),\bigl(1-|s|^2\bigr).
$$

For equal priors the answer simplifies to the famous Ivanovic–Dieks–Peres limit

$$
P_{\mathrm{succ}} = 1 - |\langle\psi_1|\psi_2\rangle| .
$$

For comparison, the minimum-error (Helstrom) success probability is

$$
P_{\mathrm{Hel}} = \frac{1}{2}\left(1 + \sqrt{1 - 4\eta_1\eta_2,|s|^2}\right),
$$

which is always larger than $P_{\mathrm{succ}}$. The difference is the price paid for certainty.

Concrete Example

We use two real qubit states in the plane,

$$
|\psi_1\rangle = \begin{pmatrix}1\0\end{pmatrix}, \qquad
|\psi_2\rangle = \begin{pmatrix}\cos\theta\ \sin\theta\end{pmatrix}, \qquad s = \cos\theta,
$$

and the orthogonal partners

$$
|\psi_2^{\perp}\rangle = \begin{pmatrix}-\sin\theta\ \cos\theta\end{pmatrix}, \qquad
|\psi_1^{\perp}\rangle = \begin{pmatrix}0\1\end{pmatrix}.
$$

Three scenarios are examined:

Case	Angle $\theta$	Prior $\eta_1$	Expected regime
A	$60^\circ$	$0.5$	Regime I, $P_{\mathrm{succ}} = 0.5$
B	$60^\circ$	$0.3$	Regime I, $P_{\mathrm{succ}} \approx 0.5417$
C	$30^\circ$	$0.2$	Regime II, $P_{\mathrm{succ}} = 0.2$

For each case we (1) evaluate the closed-form optimum, (2) verify it with a brute-force search over the feasible region, and (3) run a Monte Carlo simulation of the optimal POVM.

Python Source Code

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

plt.style.use("dark_background")


# ------------------------------------------------------------
# Theory: closed-form USD optimum and Helstrom bound
# ------------------------------------------------------------
def usd_success(eta1, x):
    """Optimal unambiguous-discrimination success probability (vectorized)."""
    eta1 = np.asarray(eta1, dtype=float)
    x = np.asarray(x, dtype=float)
    eta2 = 1.0 - eta1
    thr = np.sqrt(np.minimum(eta1, eta2) / np.maximum(eta1, eta2))
    low = 1.0 - 2.0 * np.sqrt(eta1 * eta2) * x
    high = np.maximum(eta1, eta2) * (1.0 - x ** 2)
    return np.where(x <= thr, low, high)


def helstrom_success(eta1, x):
    """Minimum-error (Helstrom) success probability (vectorized)."""
    eta1 = np.asarray(eta1, dtype=float)
    x = np.asarray(x, dtype=float)
    return 0.5 * (1.0 + np.sqrt(1.0 - 4.0 * eta1 * (1.0 - eta1) * x ** 2))


def usd_optimal_weights(eta1, s):
    """Optimal POVM weights (p1, p2) for the two zero-error elements."""
    eta1 = float(eta1)
    eta2 = 1.0 - eta1
    x = abs(float(s))
    r = np.sqrt(eta2 / eta1)
    if x <= min(r, 1.0 / r):
        p1 = (1.0 - r * x) / (1.0 - x ** 2)
        p2 = (1.0 - x / r) / (1.0 - x ** 2)
        return max(p1, 0.0), max(p2, 0.0)
    if eta1 >= eta2:
        return 1.0, 0.0
    return 0.0, 1.0


# ------------------------------------------------------------
# Brute-force verification over the (p1, p2) square
# ------------------------------------------------------------
def brute_force_success(eta1, s, n=1201):
    eta2 = 1.0 - eta1
    p = np.linspace(0.0, 1.0, n)
    P1, P2 = np.meshgrid(p, p, indexing="ij")
    tr = 2.0 - P1 - P2
    det = 1.0 - P1 - P2 + P1 * P2 * (1.0 - s ** 2)
    feas = (tr >= -1e-12) & (det >= -1e-12)
    obj = (1.0 - s ** 2) * (eta1 * P1 + eta2 * P2)
    obj_masked = np.where(feas, obj, -np.inf)
    idx = np.unravel_index(np.argmax(obj_masked), obj_masked.shape)
    return float(obj_masked[idx]), float(P1[idx]), float(P2[idx]), feas, obj, p


# ------------------------------------------------------------
# Explicit POVM and Monte Carlo simulation
# ------------------------------------------------------------
def build_povm(theta, p1, p2):
    psi1 = np.array([1.0, 0.0])
    psi2 = np.array([np.cos(theta), np.sin(theta)])
    a = np.array([-np.sin(theta), np.cos(theta)])  # orthogonal to psi2
    b = np.array([0.0, 1.0])                       # orthogonal to psi1
    E1 = p1 * np.outer(a, a)
    E2 = p2 * np.outer(b, b)
    E0 = np.eye(2) - E1 - E2
    return psi1, psi2, E1, E2, E0


def simulate(theta, eta1, n_shots, rng):
    s = np.cos(theta)
    p1, p2 = usd_optimal_weights(eta1, s)
    psi1, psi2, E1, E2, E0 = build_povm(theta, p1, p2)
    truth = np.where(rng.random(n_shots) < eta1, 1, 2)
    outcome = np.empty(n_shots, dtype=int)
    for k, psi in ((1, psi1), (2, psi2)):
        probs = np.array([psi @ E1 @ psi, psi @ E2 @ psi, psi @ E0 @ psi])
        probs = np.clip(probs, 0.0, None)
        probs = probs / probs.sum()
        mask = truth == k
        outcome[mask] = rng.choice(3, size=int(mask.sum()), p=probs)
    # outcome 0 -> declare "1", 1 -> declare "2", 2 -> inconclusive
    correct = ((outcome == 0) & (truth == 1)) | ((outcome == 1) & (truth == 2))
    error = ((outcome == 0) & (truth == 2)) | ((outcome == 1) & (truth == 1))
    inconclusive = outcome == 2
    min_eig = float(np.linalg.eigvalsh(E0).min())
    return correct.mean(), error.mean(), inconclusive.mean(), p1, p2, min_eig


# ------------------------------------------------------------
# Run the three cases
# ------------------------------------------------------------
rng = np.random.default_rng(2024)
N_SHOTS = 400_000
cases = [
    ("A", 60.0, 0.5),
    ("B", 60.0, 0.3),
    ("C", 30.0, 0.2),
]

results = []
print("=" * 118)
print("Unambiguous Quantum State Discrimination: theory vs brute force vs Monte Carlo")
print("=" * 118)
print(f"{'Case':<5}{'theta':>7}{'eta1':>7}{'|s|':>8}{'p1':>9}{'p2':>9}"
      f"{'Theory':>10}{'BruteF.':>10}{'MC succ':>10}{'MC error':>10}"
      f"{'MC inconc.':>12}{'Helstrom':>10}{'min eig':>11}")
print("-" * 118)

for name, th_deg, eta1 in cases:
    theta = np.deg2rad(th_deg)
    s = np.cos(theta)
    theory = float(usd_success(eta1, abs(s)))
    brute, bp1, bp2, _, _, _ = brute_force_success(eta1, s, n=1201)
    succ, err, inc, p1, p2, min_eig = simulate(theta, eta1, N_SHOTS, rng)
    hel = float(helstrom_success(eta1, abs(s)))
    results.append(dict(name=name, theory=theory, succ=succ, err=err, inc=inc, hel=hel))
    print(f"{name:<5}{th_deg:>7.1f}{eta1:>7.2f}{abs(s):>8.4f}{p1:>9.4f}{p2:>9.4f}"
          f"{theory:>10.4f}{brute:>10.4f}{succ:>10.4f}{err:>10.4f}"
          f"{inc:>12.4f}{hel:>10.4f}{min_eig:>11.2e}")

print("=" * 118)
print(f"Shots per case: {N_SHOTS:,}")

# ------------------------------------------------------------
# Visualization (single combined figure)
# ------------------------------------------------------------
fig = plt.figure(figsize=(22, 13))
gs = GridSpec(2, 3, figure=fig)
fig.suptitle("Unambiguous Quantum State Discrimination of Two Qubit States",
             fontsize=20, fontweight="bold")

# (0,0) 3D surface: optimal USD success probability
x_vals = np.linspace(0.0, 0.99, 70)
e_vals = np.linspace(0.02, 0.98, 70)
X, E = np.meshgrid(x_vals, e_vals)
P_usd = usd_success(E, X)
P_hel = helstrom_success(E, X)

ax1 = fig.add_subplot(gs[0, 0], projection="3d")
surf1 = ax1.plot_surface(X, E, P_usd, cmap="viridis", edgecolor="none", alpha=0.95)
ax1.set_xlabel(r"$|\langle\psi_1|\psi_2\rangle|$", labelpad=8)
ax1.set_ylabel(r"$\eta_1$", labelpad=8)
ax1.set_zlabel(r"$P_{\mathrm{succ}}$", labelpad=8)
ax1.set_title("Optimal USD success probability", fontsize=13)
ax1.view_init(elev=28, azim=-125)
fig.colorbar(surf1, ax=ax1, shrink=0.6, pad=0.1)

# (0,1) 3D surface: price of certainty (Helstrom - USD)
ax2 = fig.add_subplot(gs[0, 1], projection="3d")
surf2 = ax2.plot_surface(X, E, P_hel - P_usd, cmap="magma", edgecolor="none", alpha=0.95)
ax2.set_xlabel(r"$|\langle\psi_1|\psi_2\rangle|$", labelpad=8)
ax2.set_ylabel(r"$\eta_1$", labelpad=8)
ax2.set_zlabel(r"$P_{\mathrm{Hel}} - P_{\mathrm{succ}}$", labelpad=8)
ax2.set_title("Price of certainty (Helstrom minus USD)", fontsize=13)
ax2.view_init(elev=28, azim=-125)
fig.colorbar(surf2, ax=ax2, shrink=0.6, pad=0.1)

# (0,2) geometry of states and POVM directions (case A)
ax3 = fig.add_subplot(gs[0, 2])
th = np.deg2rad(60.0)
circle_t = np.linspace(0.0, 2.0 * np.pi, 400)
ax3.plot(np.cos(circle_t), np.sin(circle_t), color="gray", lw=1, ls=":")
vecs = [
    (r"$|\psi_1\rangle$", (1.0, 0.0), "#00e5ff"),
    (r"$|\psi_2\rangle$", (np.cos(th), np.sin(th)), "#ffb300"),
    (r"$|\psi_2^{\perp}\rangle$", (-np.sin(th), np.cos(th)), "#69f0ae"),
    (r"$|\psi_1^{\perp}\rangle$", (0.0, 1.0), "#ff4081"),
]
for label, (vx, vy), col in vecs:
    ax3.annotate("", xy=(vx, vy), xytext=(0.0, 0.0),
                 arrowprops=dict(arrowstyle="-|>", color=col, lw=2.5))
    ax3.text(1.17 * vx, 1.17 * vy, label, color=col, fontsize=15,
             ha="center", va="center")
ax3.set_xlim(-1.4, 1.4)
ax3.set_ylim(-0.4, 1.4)
ax3.set_aspect("equal")
ax3.grid(alpha=0.25)
ax3.set_title(r"State geometry (case A, $\theta=60^\circ$)", fontsize=13)

# (1,0) success probability curves vs overlap
ax4 = fig.add_subplot(gs[1, 0])
xs = np.linspace(0.0, 0.999, 500)
for eta, col in zip((0.5, 0.3, 0.1), ("#00e5ff", "#ffb300", "#ff4081")):
    ax4.plot(xs, usd_success(eta, xs), color=col, lw=2.5,
             label=rf"USD, $\eta_1={eta}$")
    ax4.plot(xs, helstrom_success(eta, xs), color=col, lw=1.8, ls="--",
             label=rf"Helstrom, $\eta_1={eta}$")
ax4.set_xlabel(r"$|\langle\psi_1|\psi_2\rangle|$")
ax4.set_ylabel("Success probability")
ax4.set_title("USD (solid) vs minimum-error (dashed)", fontsize=13)
ax4.grid(alpha=0.25)
ax4.legend(fontsize=9, loc="lower left")

# (1,1) feasible region in the (p1, p2) plane (case B)
s_b = np.cos(np.deg2rad(60.0))
eta_b = 0.3
_, _, _, feas, obj, pgrid = brute_force_success(eta_b, s_b, n=401)
Z = np.ma.masked_where(~feas, obj)
ax5 = fig.add_subplot(gs[1, 1])
mesh = ax5.pcolormesh(pgrid, pgrid, Z.T, shading="auto", cmap="plasma")
p1_line = np.linspace(0.0, 1.0, 400)
p2_line = (1.0 - p1_line) / (1.0 - (1.0 - s_b ** 2) * p1_line)
ax5.plot(p1_line, p2_line, color="white", ls="--", lw=1.8, label=r"$\det\Pi_?=0$")
bp1, bp2 = usd_optimal_weights(eta_b, s_b)
ax5.scatter([bp1], [bp2], marker="*", s=420, color="#00e5ff",
            edgecolor="white", zorder=5, label="Optimum")
ax5.set_xlabel(r"$p_1$")
ax5.set_ylabel(r"$p_2$")
ax5.set_title(r"Feasible region and objective (case B)", fontsize=13)
ax5.legend(loc="upper right", fontsize=10)
fig.colorbar(mesh, ax=ax5, shrink=0.85, label=r"$P_{\mathrm{succ}}$")

# (1,2) Monte Carlo outcome statistics
ax6 = fig.add_subplot(gs[1, 2])
idx = np.arange(len(results))
w = 0.2
series = [
    ("Simulated success", [r["succ"] for r in results], "#00e5ff"),
    ("Theoretical success", [r["theory"] for r in results], "#ffb300"),
    ("Simulated inconclusive", [r["inc"] for r in results], "#9575cd"),
    ("Simulated error", [r["err"] for r in results], "#ff4081"),
]
for k, (lab, vals, col) in enumerate(series):
    bars = ax6.bar(idx + (k - 1.5) * w, vals, width=w, color=col, label=lab)
    for bar, v in zip(bars, vals):
        ax6.text(bar.get_x() + bar.get_width() / 2.0, v + 0.01, f"{v:.3f}",
                 ha="center", va="bottom", fontsize=8, rotation=90)
ax6.set_xticks(idx)
ax6.set_xticklabels([f"Case {r['name']}" for r in results])
ax6.set_ylim(0.0, 1.2)
ax6.set_ylabel("Probability")
ax6.set_title(f"Monte Carlo outcomes ({N_SHOTS:,} shots per case)", fontsize=13)
ax6.legend(fontsize=9, loc="upper right")
ax6.grid(alpha=0.25, axis="y")

fig.subplots_adjust(left=0.04, right=0.97, top=0.92, bottom=0.07, wspace=0.28, hspace=0.30)
plt.show()

Code Walkthrough

Theory functions

usd_success implements the two-regime formula. It computes the threshold $\sqrt{\eta_{\min}/\eta_{\max}}$ with NumPy broadcasting and selects between

$$
1 - 2\sqrt{\eta_1\eta_2},|s| \qquad\text{and}\qquad \max(\eta_1,\eta_2)(1-|s|^2)
$$

through np.where. Because it accepts arrays, the same function produces both the 3D surface and the 2D curves without a single Python loop. helstrom_success evaluates the minimum-error bound for the comparison.

usd_optimal_weights returns the pair $(p_1,p_2)$. In Regime I it applies the closed-form expressions. In Regime II it sets the weight of the more probable state to $1$ and the other to $0$, which means the measurement only ever tries to identify the likelier state.

Brute-force verification

brute_force_success scans a $1201\times1201$ grid over $(p_1,p_2)$. Instead of calling an eigenvalue solver 1.4 million times, it uses the fact that a $2\times 2$ Hermitian matrix is positive semidefinite exactly when its trace and determinant are both non-negative. With $|\langle\psi_2^{\perp}|\psi_1^{\perp}\rangle| = |s|$, these are

$$
\operatorname{tr}\Pi_? = 2-p_1-p_2, \qquad \det\Pi_? = 1-p_1-p_2+p_1p_2(1-s^2).
$$

Both quantities are computed for the whole grid at once, so the search finishes almost instantly. Infeasible points are assigned $-\infty$, and np.argmax finds the best feasible point. The grid result agrees with the closed form up to the grid resolution.

POVM construction and Monte Carlo simulation

build_povm creates the explicit operators

$$
\Pi_1 = p_1|\psi_2^{\perp}\rangle\langle\psi_2^{\perp}|, \quad
\Pi_2 = p_2|\psi_1^{\perp}\rangle\langle\psi_1^{\perp}|, \quad
\Pi_? = I-\Pi_1-\Pi_2 .
$$

simulate draws the true state of every shot according to the priors. For each true state, the Born-rule outcome probabilities $\langle\psi|\Pi_j|\psi\rangle$ are computed once, and all shots of that state are sampled in a single vectorized rng.choice call. Each shot is then classified as correct, erroneous, or inconclusive with boolean masks. The smallest eigenvalue of $\Pi_?$ is reported as a sanity check that the POVM is valid. A value of order $10^{-16}$ or $0$ means the inconclusive operator sits exactly on the boundary of positivity, which is the signature of optimality.

Visualization

The single figure contains six panels:

A 3D surface of the optimal success probability over overlap and prior.
A 3D surface of the gap between the Helstrom bound and USD.
The geometry of the states and the POVM directions.
Success-probability curves for several priors.
The feasible region in the $(p_1,p_2)$ plane.
The Monte Carlo statistics for the three cases.

The vectorized design (grid evaluation with broadcasting, trace and determinant instead of eigendecomposition, and mask-based sampling) keeps the whole script fast enough to run in a few seconds.

Execution Results

======================================================================================================================
Unambiguous Quantum State Discrimination: theory vs brute force vs Monte Carlo
======================================================================================================================
Case   theta   eta1     |s|       p1       p2    Theory   BruteF.   MC succ  MC error  MC inconc.  Helstrom    min eig
----------------------------------------------------------------------------------------------------------------------
A       60.0   0.50  0.5000   0.6667   0.6667    0.5000    0.5000    0.4996    0.0000      0.5004    0.9330  -4.16e-17
B       60.0   0.30  0.5000   0.3150   0.8969    0.5417    0.5417    0.5413    0.0000      0.4587    0.9444   5.20e-17
C       30.0   0.20  0.8660   0.0000   1.0000    0.2000    0.2000    0.2006    0.0000      0.7994    0.8606   0.00e+00
======================================================================================================================
Shots per case: 400,000

Detailed Discussion of the Results

The success-probability surface

The first 3D panel shows $P_{\mathrm{succ}}$ as a function of the overlap $|s|$ and the prior $\eta_1$. Along the edge $|s|=0$ the surface sits at height $1$, because orthogonal states can be distinguished perfectly. As the overlap grows, the height falls toward $0$ along the edge $|s|\to 1$, where the two states become physically identical. The surface is highest along the central ridge $\eta_1 = 0.5$ at small overlap and bends differently near the edges $\eta_1\to 0$ and $\eta_1 \to 1$. There, the more probable state dominates and the Regime II formula $\max(\eta_1,\eta_2)(1-|s|^2)$ takes over. The kink between the two regimes is where the optimal strategy changes qualitatively.

The price of certainty

The second surface plots $P_{\mathrm{Hel}} - P_{\mathrm{succ}}$. It vanishes at $|s|=0$, since both strategies are perfect for orthogonal states. It also shrinks again near $|s|\to1$, where neither strategy can do much. The largest gap appears at intermediate overlaps and balanced priors. For case A, the theory gives $P_{\mathrm{Hel}}\approx 0.933$ against $P_{\mathrm{succ}}=0.5$. The Helstrom measurement always answers and is right about 93% of the time, whereas the unambiguous measurement answers only half the time but is never wrong.

The geometry panel

The geometry panel shows why the measurement works. The vector $|\psi_2^{\perp}\rangle$ is perpendicular to $|\psi_2\rangle$, so the projector onto it can never fire when the state is $\psi_2$. A click therefore proves the state was $\psi_1$. The same argument applies to $|\psi_1^{\perp}\rangle$ and $\psi_2$. The price is that these two directions are not orthogonal to each other, so they cannot be combined into a complete measurement. The leftover element $\Pi_?$ must absorb the remainder, and that is where the inconclusive outcomes come from.

The curves

The solid curves are the USD optimum and the dashed curves are the Helstrom bound for the same priors. Each dashed curve lies above its solid partner. For $\eta_1 = 0.5$ the solid curve is the straight line $1-|s|$. For biased priors the solid curves follow the Regime I line at first and then switch to the downward parabola $\max(\eta_1,\eta_2)(1-|s|^2)$ after the threshold overlap. The curve is continuous at the threshold, and its slope is continuous there as well, which makes the transition smooth.

The feasible region

The heatmap shows the objective over the feasible region for case B. The white dashed curve is the boundary $\det\Pi_?=0$, and the region below it is feasible. Because the objective is linear, its maximum must lie on this boundary, and the star marks the optimum, with $p_1\approx0.315$ and $p_2\approx 0.897$. The solution gives a larger weight to the less probable state’s identifier, which is a counterintuitive but correct feature of USD with unequal priors: a state that is rarely prepared needs a more aggressive unambiguous test to extract its share of the success probability.

The Monte Carlo statistics

The bar chart compares simulated and theoretical success probabilities for the three cases: $0.5$, about $0.5417$, and $0.2$. The simulated bars agree with the theoretical bars to within statistical fluctuations of order $1/\sqrt{N}$. The most important bar is the pink one. The simulated error rate is exactly zero in all three cases, because the forbidden outcomes have probability zero by construction. The rest of the probability is the inconclusive fraction, so the three bars for each case, success, inconclusive and error, add up to one.

Case C illustrates Regime II. With a small angle and a strongly biased prior, the optimal strategy abandons any attempt to identify the rare state. The measurement only tries to confirm the likely state, with $p_1=0$ and $p_2=1$ in the notation above, and it succeeds with probability $0.8\times(1-0.75)=0.2$. This is a good reminder that the “best” unambiguous measurement is not always symmetric.

Conclusion

Unambiguous state discrimination turns a fundamental quantum limitation into a clean optimization problem. By trading some certainty about whether we get an answer for absolute certainty about what the answer is, we obtain the closed-form optimum

$$
P_{\mathrm{succ}} = 1 - 2\sqrt{\eta_1\eta_2},|\langle\psi_1|\psi_2\rangle|
$$

in the balanced regime, and the simpler $\max(\eta_1,\eta_2)(1-|s|^2)$ in the strongly biased regime. Brute-force search over the feasible region and a Monte Carlo simulation of the explicit POVM both confirm the analysis, and the 3D landscapes show how the optimum depends on both the overlap and the priors. The same framework extends to more than two states, to mixed states, and to discrimination under noise, where the optimization is usually solved numerically with semidefinite programming.