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.

Minimum-Error Quantum State Discrimination as a Semidefinite Program

October 2, 2026

Trine States, Biased Priors, and a Helstrom Surface in Python

Suppose a quantum system is prepared in one of several known states $\rho_1,\dots,\rho_n$, where state $i$ occurs with prior probability $p_i$. Your job is to guess which one it was. Because non-orthogonal states cannot be distinguished perfectly, the best you can do is to minimize the probability of a wrong guess. This is minimum-error state discrimination. The optimal measurement is a POVM ${\Pi_i}$, and finding it is a textbook semidefinite program (SDP).

In this article we will:

Formulate the problem as an SDP and derive its dual.
Solve the famous trine ensemble and confirm the analytic optimum $2/3$.
Solve an asymmetric ensemble of four mixed qubit states with unequal priors, and compare the SDP with the Pretty Good Measurement (PGM).
Sweep the two-state problem over prior and overlap, and compare the SDP surface with the closed-form Helstrom bound.

1. Problem Formulation

A POVM is a set of positive semidefinite operators that sum to the identity. If we measure with ${\Pi_i}$ and answer “$i$” when outcome $i$ occurs, the success probability is

$$
P_{\mathrm{succ}} = \sum_{i=1}^{n} p_i ,\mathrm{tr}!\left(\rho_i \Pi_i\right).
$$

Maximizing it gives the SDP

$$
\begin{aligned}
\max_{\Pi_1,\dots,\Pi_n} \quad & \sum_{i=1}^{n} p_i,\mathrm{tr}(\rho_i \Pi_i) \
\text{s.t.} \quad & \Pi_i \succeq 0, \quad i=1,\dots,n, \
& \sum_{i=1}^{n} \Pi_i = I .
\end{aligned}
$$

The minimum error probability is $P_{\mathrm{err}} = 1 - P_{\mathrm{succ}}$.

Dual problem

Introduce a Hermitian multiplier $Y$ for the completeness constraint. The dual SDP is

$$
\min_{Y = Y^\dagger} \ \mathrm{tr},Y \quad \text{s.t.} \quad Y \succeq p_i \rho_i \ \ \text{for all } i .
$$

Strong duality holds, since the primal is strictly feasible with $\Pi_i = I/n$. So the optimal values coincide. This yields a certificate of optimality (the Holevo–Yuen–Kennedy–Lax conditions). A POVM is optimal if and only if

$$
Y = \sum_i p_i \rho_i \Pi_i \ \text{is Hermitian}, \qquad Y - p_i\rho_i \succeq 0 \ \ \forall i .
$$

In the code we check the second condition numerically.

Two-state special case: the Helstrom bound

For $n=2$ the SDP has a closed-form solution:

$$
P_{\mathrm{succ}} = \frac{1}{2}\left(1 + \left\lVert p_1\rho_1 - p_2\rho_2 \right\rVert_1\right).
$$

$$
P_{\mathrm{succ}} = \frac{1}{2}\left(1 + \sqrt{1 - 4p(1-p)c^2}\right).
$$

Pretty Good Measurement (baseline)

With $\rho = \sum_i p_i\rho_i$, the PGM is

$$
\Pi_i^{\mathrm{PGM}} = \rho^{-1/2}, p_i\rho_i, \rho^{-1/2}.
$$

It is easy to compute and often near-optimal, but in general it is not optimal.

2. The Three Examples

Example 1: Trine states. Three pure qubit states with equal priors $p_i = 1/3$. Their Bloch vectors lie in the $x$–$z$ plane, $120^\circ$ apart:

$$
\mathbf{r}_k = \left(\sin\theta_k,\ 0,\ \cos\theta_k\right), \qquad \theta_k = 0,\ \tfrac{2\pi}{3},\ \tfrac{4\pi}{3}.
$$

The known optimum is $\Pi_k = \tfrac{2}{3}|\psi_k\rangle\langle\psi_k|$ with $P_{\mathrm{succ}} = 2/3$.

Example 2: Asymmetric mixed ensemble. Four qubit states $\rho_k = \tfrac12(I + \mathbf{r}_k\cdot\boldsymbol{\sigma})$ with $|\mathbf{r}_k|<1$ (mixed), complex phases (nonzero $y$-components), and unequal priors $(0.40, 0.30, 0.20, 0.10)$.

Example 3: Helstrom surface. Two pure states $|\psi_0\rangle = (1,0)^T$ and $|\psi_1\rangle = (c, \sqrt{1-c^2})^T$. We solve the SDP on a grid of priors $p\in[0.02,0.98]$ and overlaps $c\in[0,1]$, then compare with the Helstrom formula.

3. Full Source Code

import time
import numpy as np
import cvxpy as cp
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm

plt.style.use("dark_background")

# ------------------------------------------------------------
# Basic quantum objects
# ------------------------------------------------------------
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)


def rho_from_bloch(r):
    r = np.asarray(r, dtype=float)
    return 0.5 * (I2 + r[0] * SX + r[1] * SY + r[2] * SZ)


def bloch_from_op(M):
    t = float(np.real(np.trace(M)))
    return np.array([
        np.real(np.trace(M @ SX)),
        np.real(np.trace(M @ SY)),
        np.real(np.trace(M @ SZ)),
    ]) / t


# ------------------------------------------------------------
# SDP for minimum-error discrimination
# ------------------------------------------------------------
SOLVER = "CLARABEL" if "CLARABEL" in cp.installed_solvers() else "SCS"


def solve_msd(priors, rhos):
    n = len(rhos)
    d = rhos[0].shape[0]
    Pi = [cp.Variable((d, d), hermitian=True) for _ in range(n)]
    constraints = [P >> 0 for P in Pi]
    constraints.append(sum(Pi) == np.eye(d))
    objective = cp.Maximize(
        cp.real(sum(float(p) * cp.trace(cp.matmul(r, P))
                    for p, r, P in zip(priors, rhos, Pi)))
    )
    problem = cp.Problem(objective, constraints)
    try:
        problem.solve(solver=SOLVER)
    except Exception:
        problem.solve(solver="SCS")

    out = []
    for P in Pi:
        M = np.array(P.value, dtype=complex)
        M = 0.5 * (M + M.conj().T)
        w, V = np.linalg.eigh(M)
        w = np.clip(w, 0.0, None)
        out.append((V * w) @ V.conj().T)
    return out


def success_prob(priors, rhos, Pis):
    return float(sum(p * np.real(np.trace(r @ P))
                     for p, r, P in zip(priors, rhos, Pis)))


def pgm(priors, rhos):
    rho_avg = sum(p * r for p, r in zip(priors, rhos))
    w, V = np.linalg.eigh(rho_avg)
    inv_sqrt = (V / np.sqrt(w)) @ V.conj().T
    return [inv_sqrt @ (p * r) @ inv_sqrt for p, r in zip(priors, rhos)]


def optimality_certificate(priors, rhos, Pis):
    Y = sum(p * r @ P for p, r, P in zip(priors, rhos, Pis))
    Y = 0.5 * (Y + Y.conj().T)
    min_eig = min(np.linalg.eigvalsh(Y - p * r).min()
                  for p, r in zip(priors, rhos))
    return float(min_eig), float(np.real(np.trace(Y)))


def report(name, priors, rhos, Pis):
    ps = success_prob(priors, rhos, Pis)
    resid = np.linalg.norm(sum(Pis) - np.eye(rhos[0].shape[0]))
    min_eig, dual_val = optimality_certificate(priors, rhos, Pis)
    print("=" * 66)
    print(name)
    print("=" * 66)
    print(f"  solver                          : {SOLVER}")
    print(f"  success probability (SDP)       : {ps:.8f}")
    print(f"  error probability               : {1.0 - ps:.8f}")
    print(f"  ||sum(Pi) - I||_F               : {resid:.2e}")
    print(f"  dual value tr(Y)                : {dual_val:.8f}")
    print(f"  min eig of (Y - p_i rho_i)      : {min_eig:.2e}  (>= 0 => optimal)")
    for k, P in enumerate(Pis):
        print(f"  tr(Pi_{k}) = {np.real(np.trace(P)):.6f}")
    return ps


# ------------------------------------------------------------
# Example 1: trine states
# ------------------------------------------------------------
t0 = time.time()
thetas = [0.0, 2.0 * np.pi / 3.0, 4.0 * np.pi / 3.0]
rhos1 = [rho_from_bloch([np.sin(t), 0.0, np.cos(t)]) for t in thetas]
priors1 = [1.0 / 3.0] * 3
Pis1 = solve_msd(priors1, rhos1)
ps1 = report("Example 1: Trine states (equal priors)", priors1, rhos1, Pis1)
print(f"  analytic optimum 2/3            : {2.0 / 3.0:.8f}")
print(f"  |SDP - 2/3|                     : {abs(ps1 - 2.0 / 3.0):.2e}")
pgm1 = success_prob(priors1, rhos1, pgm(priors1, rhos1))
guess1 = max(priors1)

# ------------------------------------------------------------
# Example 2: asymmetric mixed ensemble with unequal priors
# ------------------------------------------------------------
bloch2 = [
    [0.00, 0.00, 0.95],
    [0.80, 0.30, -0.25],
    [-0.70, 0.50, -0.20],
    [-0.10, -0.85, -0.30],
]
priors2 = [0.40, 0.30, 0.20, 0.10]
rhos2 = [rho_from_bloch(r) for r in bloch2]
Pis2 = solve_msd(priors2, rhos2)
ps2 = report("Example 2: Four mixed qubit states (unequal priors)", priors2, rhos2, Pis2)
pgm2 = success_prob(priors2, rhos2, pgm(priors2, rhos2))
guess2 = max(priors2)
print(f"  PGM success probability         : {pgm2:.8f}")
print(f"  SDP gain over PGM               : {ps2 - pgm2:.2e}")
print(f"  best blind guess (max prior)    : {guess2:.8f}")

# ------------------------------------------------------------
# Example 3: two pure states, sweep over prior p and overlap c
# ------------------------------------------------------------
def helstrom(p, c):
    return 0.5 * (1.0 + np.sqrt(np.maximum(1.0 - 4.0 * p * (1.0 - p) * c ** 2, 0.0)))


p_grid = np.linspace(0.02, 0.98, 13)
c_grid = np.linspace(0.0, 1.0, 13)
Z = np.zeros((len(p_grid), len(c_grid)))
psi0 = np.array([1.0, 0.0], dtype=complex)
for i, p in enumerate(p_grid):
    for j, c in enumerate(c_grid):
        psi1 = np.array([c, np.sqrt(max(1.0 - c ** 2, 0.0))], dtype=complex)
        rhos = [np.outer(psi0, psi0.conj()), np.outer(psi1, psi1.conj())]
        Pis = solve_msd([p, 1.0 - p], rhos)
        Z[i, j] = success_prob([p, 1.0 - p], rhos, Pis)

PP, CC = np.meshgrid(p_grid, c_grid, indexing="ij")
ZH = helstrom(PP, CC)
ERR = np.abs(Z - ZH)
print("=" * 66)
print("Example 3: Two pure states, SDP vs Helstrom bound")
print("=" * 66)
print(f"  grid size                       : {len(p_grid)} x {len(c_grid)} = {Z.size} SDPs")
print(f"  max |SDP - Helstrom|            : {ERR.max():.2e}")
print(f"  mean |SDP - Helstrom|           : {ERR.mean():.2e}")
print(f"  total elapsed time              : {time.time() - t0:.1f} s")

# ------------------------------------------------------------
# Visualization (single figure)
# ------------------------------------------------------------
def draw_bloch(ax, rhos, Pis, title):
    u = np.linspace(0, 2 * np.pi, 40)
    v = np.linspace(0, np.pi, 20)
    xs = np.outer(np.cos(u), np.sin(v))
    ys = np.outer(np.sin(u), np.sin(v))
    zs = np.outer(np.ones_like(u), np.cos(v))
    ax.plot_wireframe(xs, ys, zs, color="white", alpha=0.08, linewidth=0.5)
    colors = plt.cm.tab10(np.arange(len(rhos)))
    for k, (r, P) in enumerate(zip(rhos, Pis)):
        b = bloch_from_op(r)
        ax.scatter([b[0]], [b[1]], [b[2]], color=colors[k], s=90,
                   depthshade=False, label=f"state {k}")
        t = float(np.real(np.trace(P)))
        if t > 1e-4:
            n = bloch_from_op(P)
            nn = np.linalg.norm(n)
            if nn > 1e-9:
                d = n / nn
                ax.quiver(0, 0, 0, 1.2 * d[0], 1.2 * d[1], 1.2 * d[2],
                          color=colors[k], linewidth=2.5, arrow_length_ratio=0.12)
                ax.text(1.4 * d[0], 1.4 * d[1], 1.4 * d[2],
                        f"$\\Pi_{k}$ (tr={t:.2f})", color=colors[k], fontsize=9)
    ax.set_xlim(-1.3, 1.3)
    ax.set_ylim(-1.3, 1.3)
    ax.set_zlim(-1.3, 1.3)
    ax.set_box_aspect((1, 1, 1))
    ax.set_xlabel("x")
    ax.set_ylabel("y")
    ax.set_zlabel("z")
    ax.set_title(title, fontsize=12)
    ax.legend(loc="upper left", fontsize=8)
    ax.view_init(elev=20, azim=35)


fig = plt.figure(figsize=(21, 12))
gs = fig.add_gridspec(2, 3, wspace=0.28, hspace=0.30)

ax1 = fig.add_subplot(gs[0, 0], projection="3d")
draw_bloch(ax1, rhos1, Pis1, "(a) Trine states and optimal POVM directions")

ax2 = fig.add_subplot(gs[0, 1], projection="3d")
draw_bloch(ax2, rhos2, Pis2, "(b) Four mixed states and optimal POVM directions")

ax3 = fig.add_subplot(gs[0, 2], projection="3d")
surf = ax3.plot_surface(PP, CC, Z, cmap="viridis", edgecolor="none", alpha=0.95)
ax3.set_xlabel("prior p")
ax3.set_ylabel("overlap c")
ax3.set_zlabel("$P_{succ}$")
ax3.set_title("(c) SDP optimum $P_{succ}(p, c)$ for two pure states", fontsize=12)
ax3.view_init(elev=25, azim=-130)
fig.colorbar(surf, ax=ax3, shrink=0.55, pad=0.12)

ax4 = fig.add_subplot(gs[1, 0])
labels = ["Trine", "Mixed 4-state"]
vals_guess = [guess1, guess2]
vals_pgm = [pgm1, pgm2]
vals_sdp = [ps1, ps2]
x = np.arange(2)
w = 0.25
b1 = ax4.bar(x - w, vals_guess, w, label="Blind guess", color="#7f7f7f")
b2 = ax4.bar(x, vals_pgm, w, label="PGM", color="#ff9f43")
b3 = ax4.bar(x + w, vals_sdp, w, label="SDP (optimal)", color="#1dd1a1")
for bars in (b1, b2, b3):
    for bar in bars:
        ax4.text(bar.get_x() + bar.get_width() / 2, bar.get_height() + 0.01,
                 f"{bar.get_height():.4f}", ha="center", va="bottom", fontsize=8)
ax4.set_xticks(x)
ax4.set_xticklabels(labels)
ax4.set_ylim(0, 1.1)
ax4.set_ylabel("success probability")
ax4.set_title("(d) Strategy comparison", fontsize=12)
ax4.legend(loc="lower right", fontsize=9)
ax4.grid(alpha=0.2)

ax5 = fig.add_subplot(gs[1, 1])
vmax = max(float(ERR.max()), 1e-10)
mesh = ax5.pcolormesh(c_grid, p_grid, np.maximum(ERR, 1e-12),
                      norm=LogNorm(vmin=1e-12, vmax=vmax),
                      shading="auto", cmap="magma")
ax5.set_xlabel("overlap c")
ax5.set_ylabel("prior p")
ax5.set_title("(e) $|P^{SDP}_{succ} - P^{Helstrom}_{succ}|$ (log scale)", fontsize=12)
fig.colorbar(mesh, ax=ax5)

ax6 = fig.add_subplot(gs[1, 2])
c_fine = np.linspace(0.0, 1.0, 200)
for idx, col in zip([6, 3, 1], ["#54a0ff", "#feca57", "#ff6b6b"]):
    p_val = p_grid[idx]
    ax6.plot(c_fine, helstrom(p_val, c_fine), color=col, linewidth=2,
             label=f"Helstrom, p={p_val:.2f}")
    ax6.plot(c_grid, Z[idx, :], "o", color=col, markersize=6,
             markeredgecolor="white")
ax6.set_xlabel("overlap c")
ax6.set_ylabel("success probability")
ax6.set_title("(f) SDP points (markers) vs Helstrom curves", fontsize=12)
ax6.legend(fontsize=9)
ax6.grid(alpha=0.2)

plt.show()

4. Code Walkthrough

4.1 Building blocks

The Pauli matrices SX, SY, SZ define qubit states through the Bloch-vector representation

$$
\rho(\mathbf{r}) = \tfrac{1}{2}\left(I + \mathbf{r}\cdot\boldsymbol{\sigma}\right), \qquad |\mathbf{r}|\le 1 .
$$

Pure states have $|\mathbf{r}|=1$ and mixed states have $|\mathbf{r}|<1$. The helper rho_from_bloch builds $\rho$ from $\mathbf{r}$. The inverse helper bloch_from_op extracts $r_a = \mathrm{tr}(M\sigma_a)/\mathrm{tr}(M)$ from any positive operator. We use it to draw both states and POVM elements on the Bloch sphere.

4.2 The SDP solver `solve_msd`

This function is a direct transcription of the SDP from Section 1:

Each $\Pi_i$ is a complex Hermitian CVXPY variable, cp.Variable((d, d), hermitian=True).
P >> 0 imposes $\Pi_i\succeq 0$.
sum(Pi) == np.eye(d) imposes the completeness relation $\sum_i\Pi_i=I$.
The objective is $\mathrm{Re}\sum_i p_i,\mathrm{tr}(\rho_i\Pi_i)$. Taking the real part is harmless because the trace of a product of two Hermitian matrices is real. It just tells CVXPY the objective is real-valued.

Interior-point solvers return matrices that are Hermitian and PSD only up to numerical tolerance. So the function symmetrizes each solution and clips tiny negative eigenvalues to zero. The returned POVM elements are then clean, valid operators for all downstream computations.

4.3 Verification tools

success_prob evaluates $\sum_i p_i,\mathrm{tr}(\rho_i\Pi_i)$ directly from the returned matrices.
pgm builds the Pretty Good Measurement $\rho^{-1/2}p_i\rho_i\rho^{-1/2}$ through an eigendecomposition of $\rho=\sum_i p_i\rho_i$.
optimality_certificate constructs $Y=\sum_i p_i\rho_i\Pi_i$ and evaluates $\min_i\lambda_{\min}(Y-p_i\rho_i)$. A value that is nonnegative up to round-off certifies that the primal solution is optimal. The code also returns $\mathrm{tr},Y$, which must equal the primal value when the gap is closed.
report prints the success probability, the completeness residual $\lVert\sum_i\Pi_i-I\rVert_F$, the dual value, and the certificate.

4.4 Example 1: Trine

The three Bloch vectors $(\sin\theta_k, 0, \cos\theta_k)$ with $\theta_k\in{0, 2\pi/3, 4\pi/3}$ are pure states whose pairwise overlaps satisfy $|\langle\psi_j|\psi_k\rangle|^2 = 1/4$. By symmetry the optimal POVM is $\Pi_k=\tfrac23|\psi_k\rangle\langle\psi_k|$, which gives $P_{\mathrm{succ}}=2/3$. The code prints $|P^{\mathrm{SDP}}_{\mathrm{succ}} - 2/3|$, so you can see directly that the numerical solution matches the analytic result. A blind guess gives only $1/3$, so the measurement doubles the success rate.

4.5 Example 2: Mixed, asymmetric, complex

Here nothing is symmetric: the priors are $(0.40, 0.30, 0.20, 0.10)$, the states are mixed, and the nonzero $y$-components make the density matrices genuinely complex. No closed form is available, so the SDP is the practical tool. The comparison against the PGM shows the benefit of solving the optimization exactly. The PGM is always feasible but generally leaves some success probability on the table. A POVM element can also shrink toward zero when a state has a small prior and is hard to distinguish from its neighbours. Panel (b) draws an arrow only for elements with non-negligible trace.

4.6 Example 3: Helstrom surface

For two pure states the SDP reduces to the Helstrom problem, which gives us a ground truth for checking the solver. We solve the SDP on a $13\times13$ grid of $(p, c)$. That is 169 small $2\times2$ SDPs, each cheap enough that the whole sweep finishes quickly. Keeping the grid coarse and the matrices tiny is the main speed decision. The surface remains smooth and the agreement with theory remains clearly visible. The grid includes the extreme points $c=0$ (orthogonal states, $P_{\mathrm{succ}}=1$) and $c=1$ (identical states, $P_{\mathrm{succ}}=\max(p,1-p)$).

4.7 Visualization

All six panels are placed in a single figure:

(a), (b) Bloch-sphere plots in 3D. Dots on or inside the sphere are the input states. Arrows show the direction of each optimal POVM element, with its trace $\mathrm{tr},\Pi_k$ as a label. In the trine case, each arrow should point along its own state.
(c) The 3D surface $P_{\mathrm{succ}}(p,c)$ from the SDP.
(d) Grouped bars for blind guess, PGM, and SDP on both ensembles.
(e) A log-scale heat map of . The dark-to-bright scale makes solver-level errors visible.
(f) Cross-sections of the surface at three priors. Markers are SDP solutions and solid lines are the Helstrom formula.

5. Execution Results

==================================================================
Example 1: Trine states (equal priors)
==================================================================
  solver                          : CLARABEL
  success probability (SDP)       : 0.66666666
  error probability               : 0.33333334
  ||sum(Pi) - I||_F               : 1.24e-16
  dual value tr(Y)                : 0.66666666
  min eig of (Y - p_i rho_i)      : -7.38e-06  (>= 0 => optimal)
  tr(Pi_0) = 0.666785
  tr(Pi_1) = 0.666608
  tr(Pi_2) = 0.666608
  analytic optimum 2/3            : 0.66666667
  |SDP - 2/3|                     : 6.80e-09
==================================================================
Example 2: Four mixed qubit states (unequal priors)
==================================================================
  solver                          : CLARABEL
  success probability (SDP)       : 0.61111539
  error probability               : 0.38888461
  ||sum(Pi) - I||_F               : 1.57e-15
  dual value tr(Y)                : 0.61111539
  min eig of (Y - p_i rho_i)      : -9.99e-10  (>= 0 => optimal)
  tr(Pi_0) = 1.000000
  tr(Pi_1) = 1.000000
  tr(Pi_2) = 0.000000
  tr(Pi_3) = 0.000000
  PGM success probability         : 0.51675067
  SDP gain over PGM               : 9.44e-02
  best blind guess (max prior)    : 0.40000000
==================================================================
Example 3: Two pure states, SDP vs Helstrom bound
==================================================================
  grid size                       : 13 x 13 = 169 SDPs
  max |SDP - Helstrom|            : 2.98e-08
  mean |SDP - Helstrom|           : 6.30e-09
  total elapsed time              : 8.4 s

6. Reading the Results

Panel (a): geometry of the trine. The three POVM arrows point along the three state vectors, and each element has trace $2/3$, matching $\Pi_k = \tfrac23|\psi_k\rangle\langle\psi_k|$. The elements sum to the identity because the three weighted projectors cancel in the Pauli components and add up in the $I$ component. The measurement is a clean three-outcome “cut” of the qubit Hilbert space, and it cannot do better than $2/3$ for this ensemble.

Panel (b): asymmetric ensemble. The arrows are no longer evenly spaced. Heavily weighted, well-separated states (high prior, large Bloch-vector norm) attract POVM elements with large traces. Weak or crowded states receive small elements, and a state may be given up entirely. This is the characteristic behavior of minimum-error discrimination with unequal priors: the measurement concentrates on outcomes that pay off in expected success probability.

Panel (c): the success surface. Along the edge $c=0$ the surface sits at $1$, because orthogonal states can be identified perfectly. Along $c=1$ the states are identical and measurement gives no information, so the surface falls to $\max(p,1-p)$. That is the best you can do by guessing the more likely state. The surface is lowest near $p=1/2$ and $c=1$, where both the prior and the overlap make the two hypotheses hardest to separate.

Panel (d): why the SDP matters. The blind guess is the floor and the SDP optimum is the ceiling. The PGM sits between them and is exactly optimal for the symmetric trine. In the asymmetric ensemble, where the prior and geometry break the symmetry, the PGM can fall short of the SDP. The SDP bar is the true optimum.

Panel (e): numerical accuracy. The error heat map shows the gap between the SDP and the analytic Helstrom bound. Its magnitude reflects the convergence tolerance of the interior-point solver rather than any modeling error, which confirms that the SDP formulation reproduces the closed-form theory across the whole $(p,c)$ domain.

Panel (f): slices. The markers lie on the Helstrom curves for every prior. As $p$ moves away from $1/2$, the curves lift and flatten near $\max(p,1-p)$, showing how a strong prior makes the measurement progressively less useful.

7. Summary

Minimum-error state discrimination is one of the cleanest applications of semidefinite programming in quantum information. The primal SDP optimizes over POVMs, the dual gives a certificate through $Y\succeq p_i\rho_i$, and strong duality ties them together. With only a few lines of CVXPY we reproduced the analytic trine optimum $2/3$, solved an asymmetric mixed-state problem that has no closed form, and validated the solver against the Helstrom bound over an entire parameter surface. The same code extends to higher dimensions, more states, and general mixed ensembles by changing only the list of density matrices and priors.

Reconstructing a Two-Qubit Density Matrix from Measurement Data

October 1, 2026

Quantum State Tomography in Python

A quantum state cannot be read off a single measurement. Each measurement returns only a random outcome, and the state is destroyed in the process. What we can do is prepare the same state many times, measure it in several different bases, and reconstruct the density matrix statistically from the outcome frequencies. This procedure is called quantum state tomography.

In this article we work through a concrete example: a noisy two-qubit entangled state measured in all nine Pauli-basis combinations. We reconstruct it with two methods, linear inversion and maximum likelihood estimation (MLE), and compare them with a 3D bar chart and a shot-number scaling study.

1. The Problem Setup

A density matrix $\rho$ of an $n$-qubit system is a $2^n \times 2^n$ matrix satisfying

$$
\rho = \rho^\dagger, \qquad \rho \succeq 0, \qquad \mathrm{Tr},\rho = 1 .
$$

A two-qubit state has $4^2 - 1 = 15$ free real parameters. Any such $\rho$ can be expanded in Pauli strings:

$$
\rho = \frac{1}{4} \sum_{i,j \in {I,X,Y,Z}} \langle \sigma_i \otimes \sigma_j \rangle , \sigma_i \otimes \sigma_j .
$$

Measuring in a basis with projectors $E_k$ yields outcome $k$ with probability given by the Born rule:

$$
p_k = \mathrm{Tr}\left( E_k \rho \right).
$$

The example state. We use a Bell-type state with a relative phase, mixed with white noise:

$$
\rho_{\mathrm{true}} = p , |\psi\rangle\langle\psi| + (1-p),\frac{I}{4}, \qquad |\psi\rangle = \frac{|00\rangle + e^{i\pi/3}|11\rangle}{\sqrt{2}}, \qquad p = 0.85 .
$$

The phase $e^{i\pi/3}$ makes the off-diagonal elements genuinely complex, so both the real and imaginary parts of $\rho$ matter.

The measurement. Each qubit is measured in the $X$, $Y$, or $Z$ basis, giving $3 \times 3 = 9$ settings. Each setting has 4 outcomes, so there are 36 projectors $E_k$ in total. Every setting is repeated $N$ times, and we record the relative frequencies $f_k$.

2. Two Reconstruction Methods

Linear inversion

Since $p_k = \mathrm{Tr}(E_k \rho)$ is linear in $\rho$, replacing $p_k$ with the observed frequency $f_k$ gives a linear least-squares problem:

It is fast and simple, but statistical noise often pushes the result outside the set of valid states. The reconstructed matrix can have negative eigenvalues, which is unphysical.

Maximum likelihood estimation

MLE maximizes the log-likelihood over the set of valid density matrices:

$$
\mathcal{L}(\rho) = \sum_k f_k \log \mathrm{Tr}(E_k \rho) .
$$

A classic fixed-point algorithm, the $R\rho R$ iteration, does this without any constraint handling. Define

$$
R(\rho) = \sum_k \frac{f_k}{\mathrm{Tr}(E_k \rho)} , E_k ,
$$

and iterate

$$
\rho \leftarrow \frac{R(\rho),\rho,R(\rho)}{\mathrm{Tr}\left[R(\rho),\rho,R(\rho)\right]} , \qquad \rho_0 = \frac{I}{4} .
$$

Because $R$ is Hermitian and positive semidefinite, every iterate is automatically a valid density matrix.

Evaluation metrics

We measure the quality of a reconstruction $\sigma$ against the truth $\rho$ with the Uhlmann fidelity and the trace distance:

$$
F(\rho,\sigma) = \left( \mathrm{Tr}\sqrt{\sqrt{\rho},\sigma,\sqrt{\rho}} \right)^2, \qquad T(\rho,\sigma) = \frac{1}{2}\left| \rho - \sigma \right|_1 .
$$

3. The Complete Code

The code below is a single script that runs everything: data simulation, both reconstructions, the scaling study, and one combined figure.

The MLE loop is the only expensive part. To speed it up, all 36 projectors are flattened into a single $36 \times 16$ matrix. Each iteration then needs only two matrix-vector products instead of Python loops over projectors.

import numpy as np
import matplotlib.pyplot as plt

np.set_printoptions(precision=3, suppress=True, linewidth=120)
rng = np.random.default_rng(2024)

# ------------------------------------------------------------
# Pauli matrices and measurement projectors
# ------------------------------------------------------------
PX = np.array([[0, 1], [1, 0]], dtype=complex)
PY = np.array([[0, -1j], [1j, 0]], dtype=complex)
PZ = np.array([[1, 0], [0, -1]], dtype=complex)
PAULI = {"X": PX, "Y": PY, "Z": PZ}


def eig_projectors(P):
    """Projectors onto the +1 and -1 eigenspaces (in this order)."""
    w, v = np.linalg.eigh(P)
    order = np.argsort(-w)
    return [np.outer(v[:, k], v[:, k].conj()) for k in order]


settings = [(a, b) for a in "XYZ" for b in "XYZ"]
n_set = len(settings)

E_list = []
for a, b in settings:
    Pa = eig_projectors(PAULI[a])
    Pb = eig_projectors(PAULI[b])
    for i in range(2):
        for j in range(2):
            E_list.append(np.kron(Pa[i], Pb[j]))

E = np.array(E_list)                       # shape (36, 4, 4)
Ef = E.reshape(len(E), 16)                 # rows: vec(E_k)
A = E.transpose(0, 2, 1).reshape(len(E), 16)  # A @ vec(rho) = Tr(E_k rho)


def born_probs(rho):
    return (A @ rho.reshape(-1)).real


# ------------------------------------------------------------
# Data simulation
# ------------------------------------------------------------
def simulate(rho, shots):
    """Return relative frequencies f_k (36,) with `shots` shots per setting."""
    p = np.clip(born_probs(rho), 0.0, None).reshape(n_set, 4)
    p = p / p.sum(axis=1, keepdims=True)
    counts = np.array([rng.multinomial(shots, row) for row in p])
    return (counts / shots).reshape(-1)


# ------------------------------------------------------------
# Reconstruction 1: linear inversion (least squares)
# ------------------------------------------------------------
def linear_inversion(freq):
    vec, *_ = np.linalg.lstsq(A, freq.astype(complex), rcond=None)
    rho = vec.reshape(4, 4)
    rho = (rho + rho.conj().T) / 2
    return rho / np.trace(rho).real


# ------------------------------------------------------------
# Reconstruction 2: maximum likelihood (R rho R iteration)
# ------------------------------------------------------------
def mle(freq, max_iter=1000, tol=1e-9):
    rho = np.eye(4, dtype=complex) / 4
    n_iter = 0
    for n_iter in range(1, max_iter + 1):
        p = np.clip(born_probs(rho), 1e-12, None)
        R = ((freq / p) @ Ef).reshape(4, 4)
        new = R @ rho @ R
        new = (new + new.conj().T) / 2
        new = new / np.trace(new).real
        converged = np.linalg.norm(new - rho) < tol
        rho = new
        if converged:
            break
    return rho, n_iter


# ------------------------------------------------------------
# Metrics
# ------------------------------------------------------------
def psd_sqrt(M):
    w, v = np.linalg.eigh((M + M.conj().T) / 2)
    w = np.clip(w, 0.0, None)
    return (v * np.sqrt(w)) @ v.conj().T


def fidelity(rho, sigma):
    s = psd_sqrt(rho)
    w = np.linalg.eigvalsh(s @ sigma @ s)
    return float(np.sum(np.sqrt(np.clip(w, 0.0, None))) ** 2)


def trace_distance(rho, sigma):
    return 0.5 * float(np.sum(np.abs(np.linalg.eigvalsh(rho - sigma))))


# ------------------------------------------------------------
# True state
# ------------------------------------------------------------
p_noise = 0.85
phase = np.pi / 3
psi = np.array([1, 0, 0, np.exp(1j * phase)], dtype=complex) / np.sqrt(2)
rho_true = p_noise * np.outer(psi, psi.conj()) + (1 - p_noise) * np.eye(4) / 4

# ------------------------------------------------------------
# Single experiment with 1000 shots per setting
# ------------------------------------------------------------
shots = 1000
freq = simulate(rho_true, shots)
rho_li = linear_inversion(freq)
rho_ml, n_iter = mle(freq)

eig_true = np.linalg.eigvalsh(rho_true)[::-1]
eig_li = np.linalg.eigvalsh(rho_li)[::-1]
eig_ml = np.linalg.eigvalsh(rho_ml)[::-1]

print("=== Experiment with", shots, "shots per setting ===")
print("Purity of true state        :", round(float(np.trace(rho_true @ rho_true).real), 4))
print("Eigenvalues (true)          :", eig_true)
print("Eigenvalues (linear inv.)   :", eig_li)
print("Eigenvalues (MLE)           :", eig_ml)
print("MLE iterations              :", n_iter)
print()
print("Fidelity        LI  : %.4f   MLE : %.4f" % (fidelity(rho_true, rho_li), fidelity(rho_true, rho_ml)))
print("Trace distance  LI  : %.4f   MLE : %.4f" % (trace_distance(rho_true, rho_li), trace_distance(rho_true, rho_ml)))
print()
print("True rho (real part):")
print(rho_true.real)
print("True rho (imaginary part):")
print(rho_true.imag)
print("MLE rho (real part):")
print(rho_ml.real)
print("MLE rho (imaginary part):")
print(rho_ml.imag)

# ------------------------------------------------------------
# Scaling study: accuracy versus number of shots
# ------------------------------------------------------------
shot_list = [10, 30, 100, 300, 1000, 3000, 10000, 30000]
n_trials = 15

fid_li = np.zeros((len(shot_list), n_trials))
fid_ml = np.zeros((len(shot_list), n_trials))
td_li = np.zeros((len(shot_list), n_trials))
td_ml = np.zeros((len(shot_list), n_trials))

for i, N in enumerate(shot_list):
    for t in range(n_trials):
        f = simulate(rho_true, N)
        r_li = linear_inversion(f)
        r_ml, _ = mle(f)
        fid_li[i, t] = fidelity(rho_true, r_li)
        fid_ml[i, t] = fidelity(rho_true, r_ml)
        td_li[i, t] = trace_distance(rho_true, r_li)
        td_ml[i, t] = trace_distance(rho_true, r_ml)

print()
print("=== Scaling study (mean over %d trials) ===" % n_trials)
print("%8s | %10s %10s | %10s %10s" % ("shots", "F (LI)", "F (MLE)", "T (LI)", "T (MLE)"))
for i, N in enumerate(shot_list):
    print("%8d | %10.4f %10.4f | %10.4f %10.4f" % (
        N, fid_li[i].mean(), fid_ml[i].mean(), td_li[i].mean(), td_ml[i].mean()))

# ------------------------------------------------------------
# Visualization (single combined figure)
# ------------------------------------------------------------
labels = ["00", "01", "10", "11"]
cmap = plt.cm.viridis
zmax = 1.1 * max(np.abs(rho_true).max(), np.abs(rho_li).max(), np.abs(rho_ml).max())


def bar3d(ax, M, title):
    xs, ys = np.meshgrid(np.arange(4), np.arange(4))
    xs = xs.ravel()
    ys = ys.ravel()
    dz = np.abs(M).ravel()
    colors = cmap(dz / zmax)
    ax.bar3d(xs - 0.4, ys - 0.4, np.zeros(16), 0.8, 0.8, dz, color=colors, shade=True)
    ax.set_xticks(range(4))
    ax.set_yticks(range(4))
    ax.set_xticklabels(labels, fontsize=8)
    ax.set_yticklabels(labels, fontsize=8)
    ax.set_xlabel("column", fontsize=9)
    ax.set_ylabel("row", fontsize=9)
    ax.set_zlim(0, zmax)
    ax.set_title(title, fontsize=12)
    ax.view_init(elev=28, azim=-55)


fig = plt.figure(figsize=(18, 11))
gs = fig.add_gridspec(2, 3)

ax1 = fig.add_subplot(gs[0, 0], projection="3d")
bar3d(ax1, rho_true, r"True state $|\rho_{ij}|$")

ax2 = fig.add_subplot(gs[0, 1], projection="3d")
bar3d(ax2, rho_li, r"Linear inversion $|\rho_{ij}|$ (N=%d)" % shots)

ax3 = fig.add_subplot(gs[0, 2], projection="3d")
bar3d(ax3, rho_ml, r"MLE $|\rho_{ij}|$ (N=%d)" % shots)

ax4 = fig.add_subplot(gs[1, 0])
idx = np.arange(4)
w = 0.27
ax4.bar(idx - w, eig_true, w, label="True", color="tab:gray")
ax4.bar(idx, eig_li, w, label="Linear inversion", color="tab:red")
ax4.bar(idx + w, eig_ml, w, label="MLE", color="tab:blue")
ax4.axhline(0.0, color="black", linewidth=0.8)
ax4.set_xticks(idx)
ax4.set_xticklabels([r"$\lambda_%d$" % (k + 1) for k in range(4)])
ax4.set_ylabel("eigenvalue")
ax4.set_title("Eigenvalue spectrum (N=%d)" % shots, fontsize=12)
ax4.legend()
ax4.grid(alpha=0.3)

ax5 = fig.add_subplot(gs[1, 1])
for data, name, col in [(fid_li, "Linear inversion", "tab:red"), (fid_ml, "MLE", "tab:blue")]:
    m = data.mean(axis=1)
    s = data.std(axis=1)
    ax5.semilogx(shot_list, m, "o-", color=col, label=name)
    ax5.fill_between(shot_list, m - s, np.minimum(m + s, 1.0), color=col, alpha=0.2)
ax5.set_xlabel("shots per setting N")
ax5.set_ylabel("fidelity $F$")
ax5.set_title("Fidelity versus number of shots", fontsize=12)
ax5.legend()
ax5.grid(alpha=0.3, which="both")

ax6 = fig.add_subplot(gs[1, 2])
for data, name, col in [(td_li, "Linear inversion", "tab:red"), (td_ml, "MLE", "tab:blue")]:
    m = data.mean(axis=1)
    s = data.std(axis=1)
    ax6.loglog(shot_list, m, "o-", color=col, label=name)
    ax6.fill_between(shot_list, np.maximum(m - s, 1e-4), m + s, color=col, alpha=0.2)
ref = td_ml[0].mean() * np.sqrt(shot_list[0] / np.array(shot_list, dtype=float))
ax6.loglog(shot_list, ref, "k--", label=r"$\propto N^{-1/2}$")
ax6.set_xlabel("shots per setting N")
ax6.set_ylabel("trace distance $T$")
ax6.set_title("Trace distance versus number of shots", fontsize=12)
ax6.legend()
ax6.grid(alpha=0.3, which="both")

plt.tight_layout()
plt.show()

4. Code Walkthrough

Measurement operators. The function eig_projectors diagonalizes a Pauli matrix and returns the projectors onto its $+1$ and $-1$ eigenspaces. For each of the nine settings $(a, b)$ we build four two-qubit projectors with a Kronecker product, $E = P_a^{(\pm)} \otimes P_b^{(\pm)}$. Stacking them yields the array E of shape $(36, 4, 4)$.

The vectorization trick. Because $\mathrm{Tr}(E\rho) = \sum_{mn} E_{nm}\rho_{mn}$, the Born probabilities for all 36 outcomes come from a single matrix-vector product, A @ rho.reshape(-1). Here A holds the transposed, flattened projectors. This is the basis of the fast implementation: no Python loop over projectors appears anywhere in the hot path. The same matrix also serves as the design matrix of the linear inversion.

Simulating data. simulate computes the Born probabilities for each setting, clips tiny negative values caused by floating-point error, renormalizes, and draws one multinomial sample of size shots per setting. Dividing by shots gives the relative frequencies $f_k$.

Linear inversion. linear_inversion solves the complex least-squares problem with np.linalg.lstsq. The 36 equations for 16 unknowns are overdetermined but full rank, because the Pauli-basis projectors span the whole operator space. The result is symmetrized to enforce Hermiticity and normalized to unit trace. Nothing forces positivity, which is exactly the weakness of this method.

MLE. In mle, each iteration computes the probabilities p, then builds $R = \sum_k (f_k/p_k) E_k$ as the single product (freq / p) @ Ef. The update $R\rho R$ followed by trace normalization keeps $\rho$ positive semidefinite at every step. A floor of $10^{-12}$ on p prevents division by zero, and the loop stops once successive iterates differ by less than the tolerance.

Metrics. psd_sqrt computes a matrix square root through an eigendecomposition with negative eigenvalues clipped to zero. This is more robust than a general-purpose matrix square root routine, and it lets the fidelity be evaluated even for the non-physical linear-inversion estimate. trace_distance sums the absolute eigenvalues of $\rho - \sigma$.

Scaling study. For shot counts from 10 to 30000 we repeat the whole experiment 15 times, reconstruct with both methods, and record mean and standard deviation of fidelity and trace distance.

5. Console Output

=== Experiment with 1000 shots per setting ===
Purity of true state        : 0.7919
Eigenvalues (true)          : [0.887 0.038 0.038 0.037]
Eigenvalues (linear inv.)   : [0.888 0.06  0.043 0.009]
Eigenvalues (MLE)           : [0.888 0.059 0.04  0.013]
MLE iterations              : 601

Fidelity        LI  : 0.9868   MLE : 0.9910
Trace distance  LI  : 0.0339   MLE : 0.0304

True rho (real part):
[[0.462 0.    0.    0.212]
 [0.    0.038 0.    0.   ]
 [0.    0.    0.038 0.   ]
 [0.212 0.    0.    0.462]]
True rho (imaginary part):
[[ 0.     0.     0.    -0.368]
 [ 0.     0.     0.     0.   ]
 [ 0.     0.     0.     0.   ]
 [ 0.368  0.     0.    -0.   ]]
MLE rho (real part):
[[ 0.46  -0.     0.008  0.211]
 [-0.     0.036  0.017  0.   ]
 [ 0.008  0.017  0.036  0.01 ]
 [ 0.211  0.     0.01   0.468]]
MLE rho (imaginary part):
[[ 0.     0.008  0.003 -0.367]
 [-0.008  0.     0.012 -0.008]
 [-0.003 -0.012  0.    -0.003]
 [ 0.367  0.008  0.003  0.   ]]

=== Scaling study (mean over 15 trials) ===
   shots |     F (LI)    F (MLE) |     T (LI)    T (MLE)
      10 |     1.0250     0.8615 |     0.3810     0.2334
      30 |     0.9862     0.9132 |     0.2496     0.1832
     100 |     0.9664     0.9446 |     0.1272     0.1025
     300 |     0.9688     0.9607 |     0.0782     0.0691
    1000 |     0.9902     0.9902 |     0.0415     0.0398
    3000 |     0.9979     0.9980 |     0.0210     0.0201
   10000 |     0.9994     0.9995 |     0.0112     0.0109
   30000 |     0.9997     0.9997 |     0.0076     0.0075

6. Result Image

7. Reading the Results

The figure has six panels. The top row shows the magnitudes $|\rho_{ij}|$ of the density matrix as 3D bar charts, and the bottom row summarizes spectra and accuracy.

Top row: the density matrices. The true state has a distinctive pattern. The four diagonal bars sit at the populations, with $|00\rangle$ and $|11\rangle$ tall and $|01\rangle$, $|10\rangle$ low. The two corner bars at positions $(00, 11)$ and $(11, 00)$ are the coherences, which encode the entanglement. Both reconstructions reproduce this pattern. Linear inversion shows small spurious bars in elements that should be zero, which is statistical noise from finite sampling. MLE suppresses much of this, because the positivity constraint links elements together.

Bottom left: the eigenvalue spectrum. The true state has one large eigenvalue and three equal small ones, reflecting a pure-state component plus a uniform noise floor. The spectrum of the linear-inversion estimate is spread apart, with the largest eigenvalue pushed up and the smallest pushed down, often below zero. A negative eigenvalue means the estimate is not a legitimate quantum state. The MLE spectrum stays non-negative by construction, and the smallest eigenvalues are pulled toward the boundary.

Bottom middle: fidelity versus shots. Fidelity rises toward 1 as the number of shots grows. At small $N$ the two methods differ visibly, and the MLE curve is generally higher and has a narrower band, indicating a more stable estimator. At large $N$ the statistical noise becomes small enough that the two curves merge.

Bottom right: trace distance versus shots. On log-log axes the error falls roughly along the dashed $N^{-1/2}$ reference line. This is the standard statistical limit: the noise in each estimated frequency shrinks as $1/\sqrt{N}$, and the reconstruction error follows it. Linear inversion and MLE approach the same asymptotic behavior, but MLE has a smaller prefactor in the low-shot regime, where positivity carries real information.

Takeaways. Linear inversion is a one-line solve and is perfectly adequate when data are plentiful. Its estimate can be unphysical when data are scarce. MLE costs an iterative loop, but the vectorized implementation above keeps the cost small, and it always returns a valid density matrix. For real experiments with limited shot budgets, MLE or a related constrained estimator is the standard choice.

8. Extending the Experiment

The same code structure scales to more qubits. Only the list of measurement settings and the projector construction change, while the linear-inversion and MLE routines remain untouched. Since the number of settings grows as $3^n$ and the matrix dimension as $2^n$, the vectorized formulation becomes even more valuable. Natural next steps include replacing MLE with a Bayesian estimator to obtain error bars on every matrix element, or using compressed-sensing methods that exploit the low rank of nearly pure states to reduce the number of required settings.

Optimal Allocation of Space Weather Observation Resources

September 30, 2026

A Concrete Example Solved with Python

Introduction

A solar flare can disturb HF radio within minutes. A coronal mass ejection (CME) can trigger a geomagnetic storm a day or two later, degrading GNSS positioning, stressing power grids, and endangering satellites. Forecasters cannot watch everything all the time. Coronagraph time, magnetograph scans, ground magnetometer sampling, ionosonde sounding, GNSS-TEC processing, and L1 solar wind monitoring all have finite budgets.

This article poses a concrete allocation problem, derives its exact solution with Lagrange multipliers, and solves it in Python. Along the way we compare a generic solver with a fast semi-analytic algorithm and visualize everything in a single multi-panel figure that includes 3D graphs.

Problem Setting

We plan the next $T = 24$ hours in one-hour slots and have $I = 6$ observation resources:

Index	Resource	Role
1	Coronagraph	CME detection and speed estimation
2	Magnetograph	Active-region magnetic-field evolution
3	Ground magnetometers	Geomagnetic disturbance monitoring
4	Ionosondes	Ionospheric state
5	GNSS-TEC receivers	Total electron content mapping
6	L1 solar wind monitor	In-situ solar wind measurement

Let $x_{t,i} \ge 0$ be the observation effort given to resource $i$ in slot $t$. Each resource has:

a value weight $w_i$,
a saturation rate $a_i$ (how quickly extra effort stops adding information),
a unit cost $c_i$.

Each slot has a risk level $r_t$, meaning how much an accurate observation is worth at that hour. We model one flare-related risk pulse around hour 6 and one CME-arrival pulse around hour 15:

$$
r_t = r_0 + s\left[A_1 \exp!\left(-\frac{(t-t_1)^2}{2\sigma_1^2}\right) + A_2 \exp!\left(-\frac{(t-t_2)^2}{2\sigma_2^2}\right)\right]
$$

with $r_0 = 0.15$, $s = 1$, $(A_1, t_1, \sigma_1) = (1.0, 6, 1.5)$, and $(A_2, t_2, \sigma_2) = (1.6, 15, 2.0)$.

Mathematical Formulation

The information gained from resource $i$ in slot $t$ saturates exponentially with effort. Maximizing the total risk-weighted information gain gives

$$
\max_{x}; U(x) = \sum_{t=1}^{T}\sum_{i=1}^{I} r_t, w_i \left(1 - e^{-a_i x_{t,i}}\right)
$$

subject to a total cost budget and per-slot effort caps:

$$
\sum_{t=1}^{T}\sum_{i=1}^{I} c_i, x_{t,i} \le B, \qquad 0 \le x_{t,i} \le x_{\max}.
$$

Each term of $U$ is concave in $x_{t,i}$ and all constraints are linear, so this is a convex optimization problem. The KKT conditions are therefore both necessary and sufficient.

Closed-form structure via the Lagrangian

Introduce a multiplier $\lambda \ge 0$ for the budget constraint:

$$
\mathcal{L}(x,\lambda) = \sum_{t,i} r_t w_i\left(1-e^{-a_i x_{t,i}}\right) - \lambda\left(\sum_{t,i} c_i x_{t,i} - B\right).
$$

Setting the derivative with respect to an interior $x_{t,i}$ to zero gives

$$
r_t, w_i, a_i, e^{-a_i x_{t,i}} = \lambda, c_i .
$$

In words, at the optimum the marginal information gain per unit cost is equalized across all active (slot, resource) pairs, and that common value is $\lambda$. Solving for $x$ and applying the bounds gives

A pair $(t,i)$ receives effort only if $r_t w_i a_i / c_i > \lambda$. The total cost $\sum c_i x^{*}_{t,i}(\lambda)$ is monotonically decreasing in $\lambda$, so the $\lambda$ that exactly exhausts the budget can be found by bisection. This is a water-filling type solution.

Source Code

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

plt.style.use("dark_background")

# ------------------------------------------------------------------
# 1. Problem data
# ------------------------------------------------------------------
names = ["Coronagraph", "Magnetograph", "Magnetometers",
         "Ionosondes", "GNSS-TEC", "L1 Monitor"]
w = np.array([1.00, 0.80, 0.90, 0.60, 0.70, 1.20])   # value weights
a = np.array([0.90, 0.70, 1.10, 0.80, 0.90, 1.30])   # saturation rates
c = np.array([3.00, 2.50, 1.50, 1.20, 1.00, 4.00])   # unit costs

I = len(names)
T = 24
hours = np.arange(T)
BUDGET = 60.0
XMAX = 4.0
R0 = 0.15


def risk_profile(scale=1.0):
    flare = 1.0 * np.exp(-0.5 * ((hours - 6.0) / 1.5) ** 2)
    cme = 1.6 * np.exp(-0.5 * ((hours - 15.0) / 2.0) ** 2)
    return R0 + scale * (flare + cme)


# ------------------------------------------------------------------
# 2. Objective and gradient
# ------------------------------------------------------------------
def utility(X, r):
    return float((r[:, None] * w[None, :] * (1.0 - np.exp(-a[None, :] * X))).sum())


def gradient(X, r):
    return r[:, None] * (w * a)[None, :] * np.exp(-a[None, :] * X)


# ------------------------------------------------------------------
# 3. Generic solver (SLSQP) - reference implementation
# ------------------------------------------------------------------
def solve_slsqp(r, budget):
    cvec = np.tile(c, T)
    x0 = np.full(T * I, budget / (T * c.sum()))

    def fun(z):
        return -utility(z.reshape(T, I), r)

    def jac(z):
        return -gradient(z.reshape(T, I), r).ravel()

    cons = [{"type": "ineq",
             "fun": lambda z: budget - cvec @ z,
             "jac": lambda z: -cvec}]
    res = minimize(fun, x0, jac=jac, bounds=[(0.0, XMAX)] * (T * I),
                   constraints=cons, method="SLSQP",
                   options={"maxiter": 500, "ftol": 1e-12})
    X = np.clip(res.x.reshape(T, I), 0.0, XMAX)
    cost = float((X * c[None, :]).sum())
    if cost > budget:
        X = X * (budget / cost)
    return X


# ------------------------------------------------------------------
# 4. Fast solver: Lagrange multiplier + bisection (vectorized)
# ------------------------------------------------------------------
def allocation_for_lambda(r, lam):
    arg = np.outer(r, w * a / c) / lam
    X = np.log(np.maximum(arg, 1e-300)) / a[None, :]
    return np.clip(X, 0.0, XMAX)


def solve_fast(r, budget, iters=100):
    lo, hi = 1e-9, 1e3
    X_lo = allocation_for_lambda(r, lo)
    if float((X_lo * c[None, :]).sum()) <= budget:
        return X_lo, lo
    for _ in range(iters):
        mid = np.sqrt(lo * hi)
        X_mid = allocation_for_lambda(r, mid)
        if float((X_mid * c[None, :]).sum()) > budget:
            lo = mid
        else:
            hi = mid
    return allocation_for_lambda(r, hi), hi


# ------------------------------------------------------------------
# 5. Solve the base case and compare strategies
# ------------------------------------------------------------------
r = risk_profile(1.0)

t0 = time.perf_counter()
X_slsqp = solve_slsqp(r, BUDGET)
t_slsqp = time.perf_counter() - t0

n_rep = 50
t0 = time.perf_counter()
for _ in range(n_rep):
    X_opt, lam_opt = solve_fast(r, BUDGET)
t_fast = (time.perf_counter() - t0) / n_rep

X_uni = np.full((T, I), BUDGET / (T * c.sum()))
X_rp = np.outer(r / r.sum(), np.ones(I)) * BUDGET / c.sum()

U_uni = utility(X_uni, r)
U_rp = utility(X_rp, r)
U_slsqp = utility(X_slsqp, r)
U_opt = utility(X_opt, r)

cost_opt = float((X_opt * c[None, :]).sum())
mg = r[:, None] * (w * a / c)[None, :] * np.exp(-a[None, :] * X_opt)
interior = (X_opt > 1e-9) & (X_opt < XMAX - 1e-9)
kkt_dev = float(np.max(np.abs(mg[interior] / lam_opt - 1.0))) if interior.any() else 0.0

print("=" * 66)
print(" Space Weather Observation Resource Allocation - Results")
print("=" * 66)
print(f" Budget B                          : {BUDGET:.2f}")
print(f" Cost used (Lagrangian solution)   : {cost_opt:.4f}")
print(f" Optimal multiplier lambda*        : {lam_opt:.6f}")
print(f" Max KKT deviation (interior)      : {kkt_dev:.3e}")
print("-" * 66)
print(f" Utility  Uniform                  : {U_uni:.4f}")
print(f" Utility  Risk-proportional        : {U_rp:.4f}")
print(f" Utility  SLSQP (generic)          : {U_slsqp:.4f}")
print(f" Utility  Lagrangian bisection     : {U_opt:.4f}")
print(f" Gain over uniform                 : {100.0 * (U_opt / U_uni - 1.0):.2f} %")
print(f" Gain over risk-proportional       : {100.0 * (U_opt / U_rp - 1.0):.2f} %")
print("-" * 66)
print(f" Time  SLSQP                       : {t_slsqp:.4f} s")
print(f" Time  Lagrangian bisection        : {t_fast:.6f} s")
print(f" Speed-up                          : {t_slsqp / max(t_fast, 1e-12):.1f} x")
print("-" * 66)
print(f" {'Resource':<15}{'Total effort':>14}{'Cost share [%]':>18}")
for i in range(I):
    eff = float(X_opt[:, i].sum())
    share = 100.0 * c[i] * eff / cost_opt
    print(f" {names[i]:<15}{eff:>14.3f}{share:>18.2f}")
print("=" * 66)

# ------------------------------------------------------------------
# 6. Sensitivity surface: utility vs. budget and storm intensity
# ------------------------------------------------------------------
budgets = np.linspace(10.0, 120.0, 16)
scales = np.linspace(0.5, 3.0, 16)
U_surf = np.zeros((len(scales), len(budgets)))
for j, s in enumerate(scales):
    r_s = risk_profile(s)
    for k, b in enumerate(budgets):
        X_s, _ = solve_fast(r_s, b)
        U_surf[j, k] = utility(X_s, r_s)
Bg, Sg = np.meshgrid(budgets, scales)

# ------------------------------------------------------------------
# 7. Single combined figure
# ------------------------------------------------------------------
short = ["Corona", "Magneto", "Ground", "Iono", "GNSS", "L1"]
Ig, Hg = np.meshgrid(np.arange(I), hours)

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

ax1 = fig.add_subplot(2, 3, 1, projection="3d")
surf1 = ax1.plot_surface(Ig, Hg, X_opt, cmap="plasma", edgecolor="none", alpha=0.95)
ax1.set_xticks(np.arange(I))
ax1.set_xticklabels(short, fontsize=7)
ax1.set_ylabel("Hour")
ax1.set_zlabel("Effort $x_{t,i}$")
ax1.set_title("3D optimal allocation surface")
ax1.view_init(elev=28, azim=-125)
fig.colorbar(surf1, ax=ax1, shrink=0.55, pad=0.1)

ax2 = fig.add_subplot(2, 3, 2)
im = ax2.imshow(X_opt.T, aspect="auto", origin="lower", cmap="magma",
                extent=[-0.5, T - 0.5, -0.5, I - 0.5])
ax2.set_yticks(np.arange(I))
ax2.set_yticklabels(names, fontsize=8)
ax2.set_xlabel("Hour")
ax2.set_title("Allocation heatmap (resource x hour)")
fig.colorbar(im, ax=ax2, shrink=0.8)

ax3 = fig.add_subplot(2, 3, 3)
ax3.stackplot(hours, (X_opt * c[None, :]).T, labels=names, alpha=0.9)
ax3.set_xlabel("Hour")
ax3.set_ylabel("Cost spent per slot")
ax3.set_title("Cost composition over time")
ax3.legend(fontsize=7, loc="upper left", ncol=2)

ax4 = fig.add_subplot(2, 3, 4)
labels = ["Uniform", "Risk-\nproportional", "SLSQP", "Lagrangian\nbisection"]
vals = [U_uni, U_rp, U_slsqp, U_opt]
bars = ax4.bar(labels, vals, color=["#7f8c8d", "#3498db", "#e67e22", "#2ecc71"])
for bar, v in zip(bars, vals):
    ax4.text(bar.get_x() + bar.get_width() / 2, v, f"{v:.2f}",
             ha="center", va="bottom", fontsize=9)
ax4.set_ylim(0, max(vals) * 1.15)
ax4.set_ylabel("Total utility $U$")
ax4.set_title("Strategy comparison")

ax5 = fig.add_subplot(2, 3, 5, projection="3d")
ax5.plot_surface(Bg, Sg, U_surf, cmap="viridis", edgecolor="none", alpha=0.95)
ax5.scatter([BUDGET], [1.0], [U_opt], color="red", s=60)
ax5.set_xlabel("Budget $B$")
ax5.set_ylabel("Storm scale $s$")
ax5.set_zlabel("Optimal utility")
ax5.set_title("3D sensitivity: budget x storm intensity")
ax5.view_init(elev=25, azim=-130)

ax6 = fig.add_subplot(2, 3, 6)
ax6.plot(hours, r, color="#f1c40f", lw=2.5, marker="o", label="Risk $r_t$")
ax6.set_xlabel("Hour")
ax6.set_ylabel("Risk $r_t$", color="#f1c40f")
ax6b = ax6.twinx()
ax6b.bar(hours, (X_opt * c[None, :]).sum(axis=1), alpha=0.45, color="#1abc9c")
ax6b.set_ylabel("Total cost per slot", color="#1abc9c")
ax6.set_title("Risk profile vs. resource spending")

fig.suptitle("Optimal Allocation of Space Weather Observation Resources", fontsize=17)
fig.tight_layout(rect=[0, 0, 1, 0.96])
plt.show()

Detailed Code Walkthrough

1. Problem data. The arrays w, a, and c hold $w_i$, $a_i$, and $c_i$. The L1 monitor has the highest weight and saturation rate, since in-situ measurements are the most informative, but also the highest cost. The GNSS-TEC receivers are cheap but lower in value. risk_profile builds $r_t$ as a constant floor plus two Gaussian pulses, and the scale argument plays the role of $s$ so that storm intensity can be varied later.

2. Objective and gradient. utility evaluates $U(x)$ with fully vectorized broadcasting over the $T \times I$ matrix. gradient returns

$$
\frac{\partial U}{\partial x_{t,i}} = r_t, w_i, a_i, e^{-a_i x_{t,i}},
$$

which the generic solver needs for fast convergence.

3. Generic solver. solve_slsqp flattens the $24 \times 6$ matrix into 144 variables. np.tile(c, T) produces a cost vector that matches row-major flattening, so index $t \cdot I + i$ maps to $c_i$. The budget is passed as a linear inequality constraint with its exact Jacobian, and the bounds enforce $0 \le x \le x_{\max}$. After convergence, the solution is clipped to the bounds and rescaled if it exceeds the budget by rounding error, which guarantees feasibility.

4. Fast solver. allocation_for_lambda implements

in one vectorized expression. np.outer(r, w*a/c) builds the whole matrix of $r_t w_i a_i / c_i$, and np.maximum(arg, 1e-300) protects the logarithm. solve_fast then bisects on $\log \lambda$. The interval $[10^{-9}, 10^{3}]$ spans 12 orders of magnitude, and 100 halvings shrink it far below floating-point resolution. Because total cost is decreasing in $\lambda$, the invariant “lower end infeasible, upper end feasible” is maintained. If even $\lambda = 10^{-9}$ fits within the budget, every effort is already saturated and that allocation is returned directly.

5. Comparison and diagnostics. Two intuitive baselines are computed alongside the optimum: uniform allocation and allocation proportional to risk. Both spend exactly the budget $B$, so the comparison is fair. The KKT check evaluates the marginal gain per unit cost

$$
\frac{r_t w_i a_i e^{-a_i x_{t,i}}}{c_i}
$$

on all interior pairs and confirms that it equals $\lambda^{*}$ up to floating-point error.

6. Sensitivity surface. The fast solver is called for a $16 \times 16$ grid of budgets and storm scales, 256 full optimizations in total. This is only practical because each solve takes a fraction of a millisecond.

7. Figure. All six panels are drawn into one figure and shown with a single plt.show().

Making It Fast

The generic SLSQP solver treats the problem as a black-box nonlinear program with 144 variables. Each iteration solves a dense quadratic subproblem whose cost grows roughly cubically with the number of variables. It also needs many iterations, and its runtime grows quickly as $T$ or $I$ increases.

The Lagrangian approach exploits the structure of the problem. The coupling between all $T \times I$ variables happens only through one scalar constraint, so the problem collapses to a one-dimensional root-finding problem in $\lambda$. Each bisection step is a single vectorized array operation, and the complexity is $O(T \cdot I \cdot \text{iters})$, which is linear in problem size. This is why the 256-point sensitivity surface costs almost nothing. The same idea scales to thousands of slots and hundreds of instruments.

Execution Results

==================================================================
 Space Weather Observation Resource Allocation - Results
==================================================================
 Budget B                          : 60.00
 Cost used (Lagrangian solution)   : 60.0000
 Optimal multiplier lambda*        : 0.266670
 Max KKT deviation (interior)      : 2.220e-16
------------------------------------------------------------------
 Utility  Uniform                  : 13.5547
 Utility  Risk-proportional        : 19.9097
 Utility  SLSQP (generic)          : 24.9295
 Utility  Lagrangian bisection     : 24.9295
 Gain over uniform                 : 83.92 %
 Gain over risk-proportional       : 25.21 %
------------------------------------------------------------------
 Time  SLSQP                       : 0.4205 s
 Time  Lagrangian bisection        : 0.003741 s
 Speed-up                          : 112.4 x
------------------------------------------------------------------
 Resource         Total effort    Cost share [%]
 Coronagraph             2.956             14.78
 Magnetograph            1.326              5.53
 Magnetometers           9.674             24.19
 Ionosondes              6.213             12.43
 GNSS-TEC               11.204             18.67
 L1 Monitor              3.661             24.41
==================================================================

Interpreting the Results

Console output. The cost used should match the budget $B$ almost exactly, which confirms that the budget constraint is active. The KKT deviation should be on the order of machine precision, showing that all active pairs share the same marginal value per unit cost $\lambda^{*}$. The utility list ranks the four strategies. The Lagrangian solution is the exact optimum of the convex problem, so SLSQP can match it but never beat it beyond numerical tolerance, and both must exceed the two heuristics. The timing lines quantify the benefit of exploiting problem structure. The per-resource table shows where the budget actually goes.

Panel 1 (3D allocation surface). Ridges appear at the two risk pulses, around hours 6 and 15. Along the resource axis the surface is highest for resources with a large ratio $w_i a_i / c_i$, because those cross the activation threshold first and are then pushed further by the logarithmic growth of $x^{*}$. Expensive resources appear mainly at the ridge peaks, and cheap resources also cover the quiet hours.

Panel 2 (heatmap). This is the same information from above. Bright vertical bands mark the flare and CME windows. The dark cells show pairs where the threshold is not met, meaning it is optimal to spend nothing there. This sparsity pattern is exactly the water-filling structure of the KKT solution.

Panel 3 (cost composition). The stacked areas show how the spending per hour is split among resources. The total height follows the risk profile, and the colored layers show that additional resources are “switched on” progressively as risk rises, rather than all scaling together.

Panel 4 (strategy comparison). Uniform allocation ignores when observations matter. Risk-proportional allocation ignores which resources are cost-effective and how quickly each saturates. The optimum accounts for both, which is why it clearly outperforms both baselines.

Panel 5 (3D sensitivity surface). Utility increases along both axes, but with different curvature. Along the budget axis the surface flattens, showing diminishing returns: the exponential saturation makes each extra unit of budget less valuable. Along the storm-scale axis, stronger events raise the achievable utility because more risk-weighted information is available to capture. The red marker shows the base case. The surface answers practical questions such as how much additional budget is needed to keep the same utility during a more active period.

Panel 6 (risk versus spending). The yellow curve is the risk profile and the teal bars are the total cost spent per hour. The bars follow the curve but are not proportional to it. Because of the concave utility, the optimal spending rises more gently than the risk itself, which avoids over-investing in a single hour.

Conclusion

The allocation of space weather observation resources can be written as a convex program with exponential-saturation utilities. Its KKT conditions yield a closed-form water-filling solution parametrized by a single Lagrange multiplier, which can be found by bisection in a fraction of a millisecond. Compared with uniform and risk-proportional allocation, the optimal plan concentrates effort where risk is high and where information is cheap, and the sensitivity surface shows how the optimal value responds to budget and storm intensity. The same framework extends naturally to additional instruments, finer time resolution, or richer risk forecasts.

Optimizing the Timing of Space Weather Warnings

September 29, 2026

A Worked Example in Python

Space weather forecasters face a very concrete trade-off every time a fast solar wind stream or a CME-driven shock approaches Earth: warn too early on noisy data and you get false alarms that erode public trust; wait too long to be sure and the warning becomes useless because the storm has already arrived. This article works through a simplified but physically grounded version of that problem, builds a cost-function optimization for it in Python, and visualizes the result — including a 3D cost surface.

1. The physical setup

Real-time solar wind monitors (like the ones sitting at the L1 Lagrange point, about 1.5 million km upstream of Earth) give operators a finite warning window before the same plasma reaches the magnetosphere. The nominal propagation time is

$$\tau = \frac{D_{L1}}{v_{sw}}$$

where $D_{L1} \approx 1.5\times10^{6},\text{km}$ and $v_{sw}$ is the solar wind speed. Faster wind means less warning time — a 400 km/s wind gives roughly an hour of lead time, while an 800 km/s shock front gives barely half that.

The classic trigger condition for a geomagnetic storm warning is a sustained southward turning of the interplanetary magnetic field ($B_z \ll 0$), because southward $B_z$ reconnects efficiently with Earth’s northward-pointing field and drives geomagnetic activity. In practice, raw $B_z$ measurements are noisy, so operators smooth the signal before comparing it against a threshold — but smoothing itself eats into the available lead time.

This gives us two knobs to tune:

Smoothing window length $w$ (minutes) — larger $w$ suppresses noise but adds detection delay.
Detection threshold $\theta$ (nT) — a less negative threshold triggers sooner but generates more false alarms.

2. Formulating the optimization problem

For a given $(w, \theta)$ pair, three quantities matter operationally:

$P_{FA}(w,\theta)$ — probability that a quiet (non-storm) period crosses the threshold and triggers a false alarm.
$P_{MD}(w,\theta)$ — probability that a real storm is missed (never crosses the threshold before impact).
$\overline{L}(w,\theta)$ — the mean effective lead time, in minutes, among correctly detected storms.

We combine these into a single operational cost function to minimize:

$$C(w,\theta) = \alpha, P_{FA}(w,\theta) ;+; \beta, P_{MD}(w,\theta) ;-; \gamma, \frac{\overline{L}(w,\theta)}{L_{max}}$$

Here $\alpha, \beta, \gamma$ are relative weights (missed storms are usually considered worse than false alarms, so $\beta > \alpha$), and $L_{max}$ is the full observation window used for normalization. The optimal warning strategy is then

To generate realistic southward-turning events for the simulation, each synthetic storm’s $B_z$ ramp follows a logistic (sigmoid) onset:

$$B_z(t) = B_{z,\min}\cdot \sigma!\left(\frac{t - t_0}{\tau_{ramp}/4}\right), \qquad \sigma(x) = \frac{1}{1+e^{-x}}$$

on top of an AR(1) noise process representing normal solar wind fluctuations.

3. Python implementation

The code below simulates 3,000 synthetic solar wind passages (1,500 quiet, 1,500 storm), performs a full grid search over smoothing windows and thresholds, computes the cost surface, and produces a combined diagnostic figure. Everything is vectorized with NumPy — the only Python-level loops are over the (small) parameter grids themselves, not over individual events, so the whole search runs in a few seconds even in the free Colab tier.

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D  # noqa: F401 (enables 3D projection)
from matplotlib import cm

# ---------------------------------------------------------------
# 1. Reproducibility
# ---------------------------------------------------------------
rng = np.random.default_rng(42)

# ---------------------------------------------------------------
# 2. Simulation parameters
# ---------------------------------------------------------------
N_QUIET = 1500
N_STORM = 1500
N_TOTAL = N_QUIET + N_STORM
T_OBS   = 180          # observation window before Earth impact [minutes]

D_L1 = 1.5e6           # Earth-L1 distance [km]

# ---------------------------------------------------------------
# 3. Vectorized AR(1) noise generator (one Python loop over time,
#    fully vectorized across all events simultaneously)
# ---------------------------------------------------------------
def generate_ar1_noise(n_events, n_steps, phi, sigma, rng):
    eps = rng.normal(0.0, sigma, size=(n_events, n_steps))
    noise = np.empty_like(eps)
    noise[:, 0] = eps[:, 0]
    for t in range(1, n_steps):
        noise[:, t] = phi * noise[:, t - 1] + eps[:, t]
    return noise

# ---------------------------------------------------------------
# 4. Synthetic quiet-time Bz and solar wind speed
# ---------------------------------------------------------------
Bz_quiet = generate_ar1_noise(N_QUIET, T_OBS, phi=0.9, sigma=0.9, rng=rng)
v_quiet  = np.clip(rng.normal(400, 40, size=(N_QUIET, T_OBS)), 300, 550)

# ---------------------------------------------------------------
# 5. Synthetic storm events: AR(1) noise + logistic southward ramp
# ---------------------------------------------------------------
Bz_storm_noise = generate_ar1_noise(N_STORM, T_OBS, phi=0.9, sigma=0.9, rng=rng)

t_axis     = np.arange(T_OBS)
onset_time = rng.uniform(20, 150, size=(N_STORM, 1))   # ramp start [min before impact]
ramp_dur   = rng.uniform(15, 50, size=(N_STORM, 1))    # ramp duration [min]
Bz_min     = rng.uniform(-22, -8, size=(N_STORM, 1))   # most negative Bz reached [nT]

sigmoid_arg = (t_axis[None, :] - onset_time) / (ramp_dur / 4.0)
ramp = Bz_min / (1.0 + np.exp(-sigmoid_arg))
Bz_storm = Bz_storm_noise + ramp

v_storm_base = rng.normal(420, 40, size=(N_STORM, T_OBS))
speed_gain   = rng.uniform(150, 350, size=(N_STORM, 1))
v_storm = np.clip(v_storm_base + speed_gain / (1.0 + np.exp(-sigmoid_arg)), 300, 900)

Bz_all = np.vstack([Bz_quiet, Bz_storm])
v_all  = np.vstack([v_quiet, v_storm])

# ---------------------------------------------------------------
# 6. Reference propagation times (background physics)
# ---------------------------------------------------------------
print("Reference L1 -> Earth propagation time (tau = D_L1 / v_sw):")
for v_ref in (400, 600, 800):
    tau_min = (D_L1 / v_ref) / 60.0
    print(f"  v_sw = {v_ref} km/s  ->  tau approx {tau_min:.1f} min")

# ---------------------------------------------------------------
# 7. Causal moving average (uses only past/current samples,
#    matching a real-time operational filter)
# ---------------------------------------------------------------
def causal_moving_average(x, w):
    n = x.shape[1]
    csum = np.cumsum(x, axis=1)
    csum = np.concatenate([np.zeros((x.shape[0], 1)), csum], axis=1)
    out = np.empty_like(x)
    for t in range(n):
        lo = max(0, t - w + 1)
        out[:, t] = (csum[:, t + 1] - csum[:, lo]) / (t - lo + 1)
    return out

# ---------------------------------------------------------------
# 8. Grid search over (window w, threshold theta)
# ---------------------------------------------------------------
W_GRID     = np.arange(1, 25, 2)               # smoothing window [min]
THETA_GRID = np.linspace(-3.0, -13.0, 21)       # Bz threshold [nT]
n_w, n_th  = len(W_GRID), len(THETA_GRID)

ALPHA, BETA, GAMMA = 1.0, 2.0, 1.5   # cost weights: false-alarm, missed-storm, lead-time reward

P_FA = np.zeros((n_w, n_th))
P_MD = np.zeros((n_w, n_th))
LEAD = np.zeros((n_w, n_th))

for i, w in enumerate(W_GRID):
    w_int = int(w)
    smoothed = causal_moving_average(Bz_all, w_int)
    delay = (w_int - 1) / 2.0   # average delay introduced by the causal filter

    crossed   = smoothed[:, :, None] <= THETA_GRID[None, None, :]   # (N, T, M)
    any_cross = crossed.any(axis=1)                                  # (N, M)
    first_idx = crossed.argmax(axis=1).astype(float)                 # (N, M)

    for j in range(n_th):
        trig_quiet = any_cross[:N_QUIET, j]
        trig_storm = any_cross[N_QUIET:, j]

        P_FA[i, j] = trig_quiet.mean()
        P_MD[i, j] = 1.0 - trig_storm.mean()

        if trig_storm.any():
            idx_detected = first_idx[N_QUIET:, j][trig_storm]
            lead_time = np.clip((T_OBS - idx_detected) - delay, 0, None)
            LEAD[i, j] = lead_time.mean()

COST = ALPHA * P_FA + BETA * P_MD - GAMMA * (LEAD / T_OBS)

i_opt, j_opt = np.unravel_index(np.argmin(COST), COST.shape)
w_opt, theta_opt = W_GRID[i_opt], THETA_GRID[j_opt]

print("\n=== Optimal warning-timing parameters ===")
print(f"optimal smoothing window   w*      = {w_opt} min")
print(f"optimal Bz threshold       theta*  = {theta_opt:.2f} nT")
print(f"minimum cost                C*      = {COST[i_opt, j_opt]:.4f}")
print(f"false alarm rate at optimum         = {P_FA[i_opt, j_opt]*100:.2f} %")
print(f"missed detection rate at optimum    = {P_MD[i_opt, j_opt]*100:.2f} %")
print(f"mean effective lead time at optimum = {LEAD[i_opt, j_opt]:.1f} min")

# ---------------------------------------------------------------
# 9. Visualization (single combined figure)
# ---------------------------------------------------------------
fig = plt.figure(figsize=(18, 11))
gs = fig.add_gridspec(2, 3, hspace=0.4, wspace=0.35)

# -- 3D cost surface --
ax1 = fig.add_subplot(gs[0, 0], projection='3d')
W_mesh, TH_mesh = np.meshgrid(W_GRID, THETA_GRID, indexing='ij')
surf = ax1.plot_surface(W_mesh, TH_mesh, COST, cmap=cm.viridis, edgecolor='none', alpha=0.9)
ax1.scatter([w_opt], [theta_opt], [COST[i_opt, j_opt]], color='red', s=60)
ax1.set_xlabel('window w [min]')
ax1.set_ylabel('threshold θ [nT]')
ax1.set_zlabel('cost C(w,θ)')
ax1.set_title('Cost surface C(w, θ)')
fig.colorbar(surf, ax=ax1, shrink=0.6, pad=0.1)

# -- top-view contour --
ax2 = fig.add_subplot(gs[0, 1])
im = ax2.contourf(THETA_GRID, W_GRID, COST, levels=20, cmap='viridis')
ax2.plot(theta_opt, w_opt, 'r*', markersize=16)
ax2.set_xlabel('threshold θ [nT]')
ax2.set_ylabel('window w [min]')
ax2.set_title('Cost contour (top view)')
fig.colorbar(im, ax=ax2)

# -- ROC-like curve at w_opt --
ax3 = fig.add_subplot(gs[0, 2])
det_rate = 1 - P_MD[i_opt, :]
ax3.plot(P_FA[i_opt, :], det_rate, 'o-', color='tab:blue')
ax3.plot([0, 1], [0, 1], 'k--', alpha=0.4)
ax3.set_xlabel('false alarm rate')
ax3.set_ylabel('detection rate')
ax3.set_title(f'ROC-like curve (w = {w_opt} min)')

# -- lead-time distribution at optimum --
w_int = int(w_opt)
smoothed_opt = causal_moving_average(Bz_all, w_int)
delay_opt = (w_int - 1) / 2.0
crossed_opt   = smoothed_opt <= theta_opt
any_cross_opt = crossed_opt.any(axis=1)
first_idx_opt = crossed_opt.argmax(axis=1).astype(float)
storm_trig = any_cross_opt[N_QUIET:]
lead_all = np.clip((T_OBS - first_idx_opt[N_QUIET:][storm_trig]) - delay_opt, 0, None)

ax4 = fig.add_subplot(gs[1, 0])
ax4.hist(lead_all, bins=25, color='tab:orange', edgecolor='k', alpha=0.8)
ax4.axvline(lead_all.mean(), color='red', linestyle='--', label=f'mean = {lead_all.mean():.1f} min')
ax4.set_xlabel('effective lead time [min]')
ax4.set_ylabel('count')
ax4.set_title('Lead-time distribution at optimum')
ax4.legend()

# -- example raw vs. smoothed Bz with threshold --
ax5 = fig.add_subplot(gs[1, 1:])
q_idx, s_idx = 0, N_QUIET + 3
ax5.plot(t_axis, Bz_all[q_idx], color='gray', alpha=0.6, label='quiet Bz (raw)')
ax5.plot(t_axis, Bz_all[s_idx], color='tab:blue', alpha=0.4, label='storm Bz (raw)')
ax5.plot(t_axis, smoothed_opt[s_idx], color='tab:red', linewidth=2,
         label=f'storm Bz (smoothed, w={w_int})')
ax5.axhline(theta_opt, color='k', linestyle='--', label=f'threshold θ={theta_opt:.1f} nT')
ax5.set_xlabel('time [min from window start]')
ax5.set_ylabel('Bz [nT]')
ax5.set_title('Example raw / smoothed Bz and detection threshold')
ax5.legend(loc='lower left', fontsize=8)

plt.suptitle('Space-Weather Warning Timing Optimization', fontsize=16, y=1.02)
plt.show()

4. Walking through the code

Data simulation (sections 4–5). Rather than pulling live satellite data (which would make the article non-reproducible), the script generates 1,500 “quiet” and 1,500 “storm” synthetic passages. Quiet-time $B_z$ is a pure AR(1) noise process centered on zero, mimicking the natural jitter of the interplanetary magnetic field. Storm events add a logistic southward ramp on top of the same noise process, with a randomized onset time, ramp duration, and minimum $B_z$ — so every simulated storm looks slightly different, just like real events. Solar wind speed is generated the same way, rising together with the $B_z$ ramp to represent the arrival of a fast stream or shock.

Causal smoothing (section 7). A real-time forecasting system cannot use a centered moving average, because that would require future data it doesn’t have yet. causal_moving_average uses a cumulative-sum trick to compute, at each time step, the average of only the current and preceding samples — this is what an operational monitor actually sees. The associated delay is $\approx (w-1)/2$ samples, which is subtracted from the raw lead time to get the effective lead time.

Grid search (section 8). For each smoothing window $w$ in W_GRID, the whole ensemble of 3,000 time series is smoothed once, and then compared against every threshold in THETA_GRID simultaneously using NumPy broadcasting (smoothed[:, :, None] <= THETA_GRID[None, None, :]). This produces a 3D boolean array (events × time × thresholds) in one shot, letting argmax find each event’s first threshold-crossing instant for every threshold value without any Python-level loop over individual events. Only the outer loop over $w$ (12 values) and an inner loop over $\theta$ (21 values) run in Python — 252 iterations total, each doing large vectorized array operations, which is why the whole grid search finishes in a few seconds rather than minutes. This vectorized structure is effectively the “fast” version of what would otherwise be a triple nested loop over events, time, and thresholds.

Cost evaluation. For every $(w,\theta)$ cell, the script computes the false-alarm rate among quiet events, the missed-detection rate among storm events, and the mean effective lead time among storms that were correctly detected, then combines them into $C(w,\theta)$ using the weights $\alpha=1.0$, $\beta=2.0$, $\gamma=1.5$ — reflecting the operational reality that missing a real storm is treated as roughly twice as costly as a false alarm, while longer lead time is explicitly rewarded.

Finding the optimum. np.argmin(COST) combined with np.unravel_index locates the $(w^{}, \theta^{})$ pair that minimizes the full cost surface — this is the numerical solution to the optimization problem stated in Section 2.

5. Reading the plots

Top-left (3D surface): the full cost landscape $C(w,\theta)$. You should see a bowl-shaped region — too small a window (noisy, high false-alarm cost) and too conservative a threshold (high missed-detection cost) both push the cost up, with a visible minimum marked in red.
Top-middle (contour): the same surface viewed from above, useful for reading off the optimal region at a glance.
Top-right (ROC-like curve): at the optimal window, this shows how detection rate trades off against false-alarm rate as the threshold varies — the classic signal-detection trade-off curve.
Bottom-left (histogram): the spread of effective lead times actually achieved at the optimal settings — useful for communicating “how much warning time” to stakeholders, not just an average number.
Bottom-right (example traces): a concrete illustration of one storm’s raw vs. smoothed $B_z$ crossing the optimal threshold, which makes the abstract optimization tangible.

Reference L1 -> Earth propagation time (tau = D_L1 / v_sw):
  v_sw = 400 km/s  ->  tau approx 62.5 min
  v_sw = 600 km/s  ->  tau approx 41.7 min
  v_sw = 800 km/s  ->  tau approx 31.2 min

=== Optimal warning-timing parameters ===
optimal smoothing window   w*      = 1 min
optimal Bz threshold       theta*  = -7.00 nT
minimum cost                C*      = -0.7656
false alarm rate at optimum         = 2.20 %
missed detection rate at optimum    = 0.13 %
mean effective lead time at optimum = 94.8 min

6. Takeaways and possible extensions

Even in this simplified model, the optimization consistently favors a short-to-moderate smoothing window paired with a moderately conservative threshold — reflecting the fundamental tension between noise suppression and warning latency. In a real operational system, this framework could be extended by:

Using actual DSCOVR/ACE $B_z$ and speed time series instead of synthetic data.
Making $\alpha, \beta, \gamma$ time-varying (e.g., raising the missed-detection penalty during solar maximum).
Replacing the grid search with a continuous optimizer (scipy.optimize.minimize) once the cost surface is confirmed to be reasonably smooth, as seen in the 3D plot above.
Incorporating multiple correlated parameters (e.g., $B_z$, speed, and density jointly) rather than a single threshold variable.

Minimizing Geomagnetically Induced Current (GIC) Risk in Power Grids

September 28, 2026

A Hands-On Python Example

Geomagnetically Induced Currents (GIC) are quasi-DC currents driven into power transmission networks during geomagnetic storms. When the Sun ejects a coronal mass ejection toward Earth, the resulting disturbance in the geomagnetic field induces an electric field at ground level. Long transmission lines act like giant antennas for this field, and the induced currents flow into the earth through transformer neutral grounding points. Because power transformers are designed for 50/60 Hz AC, this quasi-DC current can saturate the transformer core half-cycle by half-cycle, causing overheating, harmonic distortion, and — in the worst historical cases (Quebec 1989, Halloween Storm 2003) — cascading blackouts.

Utilities mitigate this risk by installing series blocking devices (capacitors or resistors) in the transformer neutral-to-ground path at selected substations. Since blocking devices are expensive, the practical engineering question is:

Given a limited budget of k blocking devices, which substations should receive them so that the worst-case GIC across every possible storm direction is minimized?

This article walks through a complete, runnable example: modeling a small transmission network with the classic Lehtinen–Pirjola (LP) method, then solving the placement problem as a robust (minimax) combinatorial optimization.

1. The Physics: Lehtinen–Pirjola Network Model

A transmission network is modeled as a resistive DC circuit. Each substation i has a transformer neutral grounding resistance $R_{g,i}$, and each transmission line $(i,j)$ has a series resistance $R_{ij}$ carrying an induced EMF $E_{ij}$ caused by the geoelectric field.

The nodal voltage equation (Kirchhoff’s Current Law at every substation) is:

$$
(\mathbf{Y} + \mathbf{Y}_e),\mathbf{V} = \mathbf{J}
$$

where $\mathbf{Y}$ is the network admittance matrix built from line conductances $g_{ij} = 1/R_{ij}$, $\mathbf{Y}e = \mathrm{diag}(1/R{g,i})$ is the grounding admittance, and $\mathbf{J}$ is the current injection vector caused by the storm’s electric field:

$$
J_i = \sum_{j \in \mathcal{N}(i)} g_{ij}, \big(\mathbf{E}\cdot \mathbf{L}_{ij}\big)
$$

Here $\mathbf{E}$ is the uniform horizontal geoelectric field vector (V/km) and $\mathbf{L}_{ij}$ is the displacement vector (in km) from node i to node j.

Once $\mathbf{V}$ is solved, the GIC flowing through each transformer’s neutral (the quantity that damages equipment) is simply:

$$
I_{g,i} = \frac{V_i}{R_{g,i}}
$$

A blocking device is modeled by setting $R_{g,i}$ to an effectively infinite value, which drives $I_{g,i}$ toward zero at that node.

The Key Trick: Linearity in Storm Direction

If the storm direction is $\theta$, then $\mathbf{E} = E_0(\cos\theta, \sin\theta)$, and since $\mathbf{J}$ is linear in $\mathbf{E}$:

$$
\mathbf{J}(\theta) = \cos\theta \cdot \mathbf{A} + \sin\theta \cdot \mathbf{B}
$$

where $\mathbf{A}$ and $\mathbf{B}$ are constant vectors (independent of $\theta$). Because the system matrix $(\mathbf{Y}+\mathbf{Y}_e)$ does not depend on $\theta$ either, the solution is also linear in $\theta$:

$$
\mathbf{V}(\theta) = \cos\theta \cdot \mathbf{V}_A + \sin\theta \cdot \mathbf{V}_B
$$

This means we never need to re-solve the linear system for every storm angle — we solve it twice (for $\mathbf{A}$ and $\mathbf{B}$) and reconstruct every angle by simple trigonometric scaling. This is the speed optimization used below.

2. The Concrete Example

We model an 8-substation regional grid (a ring topology with a central hub, which is a common real-world structure). We test every 1-degree storm direction (360 angles) and evaluate every possible placement of 2 blocking devices ($\binom{8}{2} = 28$ combinations), selecting the placement that minimizes the worst-case GIC across all storm directions and all unprotected substations.

3. Python Implementation (Google Colab Ready)

# ============================================================
# GIC (Geomagnetically Induced Current) Risk Minimization
# Lehtinen-Pirjola network model + robust blocking-device placement
# ============================================================

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from itertools import combinations
import time

np.random.seed(42)

# ---------------------------------------------------------------
# 1. Define the power grid topology (8 substations)
# ---------------------------------------------------------------
n_nodes = 8

coords = np.array([
    [0.0, 0.0],
    [150.0, 40.0],
    [280.0, 0.0],
    [320.0, 160.0],
    [200.0, 260.0],
    [60.0, 220.0],
    [-40.0, 120.0],
    [140.0, 130.0],   # central hub
])

edges = [
    (0, 1), (1, 2), (2, 3), (3, 4), (4, 5), (5, 6), (6, 0),
    (7, 0), (7, 2), (7, 4), (7, 6)
]

r_per_km = 0.045  # ohm/km, typical EHV line resistance

def line_length(i, j):
    return np.linalg.norm(coords[j] - coords[i])

R_line = {e: r_per_km * line_length(*e) for e in edges}

# ---------------------------------------------------------------
# 2. Baseline transformer neutral grounding resistances
# ---------------------------------------------------------------
R_ground_base = np.full(n_nodes, 0.6)   # ohm, typical neutral grounding
R_block = 1.0e6                          # ohm, effective open circuit

E0 = 8.0          # V/km, severe geomagnetic storm field magnitude
n_angles = 360
thetas = np.linspace(0, 2*np.pi, n_angles, endpoint=False)

# ---------------------------------------------------------------
# 3. Build network admittance matrix Y (independent of storm angle)
# ---------------------------------------------------------------
Y = np.zeros((n_nodes, n_nodes))
Lx = np.zeros(n_nodes)
Ly = np.zeros(n_nodes)

for (i, j) in edges:
    g = 1.0 / R_line[(i, j)]
    Y[i, i] += g
    Y[j, j] += g
    Y[i, j] -= g
    Y[j, i] -= g

    dx = coords[j, 0] - coords[i, 0]
    dy = coords[j, 1] - coords[i, 1]

    Lx[i] += g * dx;  Ly[i] += g * dy
    Lx[j] -= g * dx;  Ly[j] -= g * dy

A = E0 * Lx
B = E0 * Ly

# ---------------------------------------------------------------
# 4. Fast solver exploiting linearity in theta
#    V(theta) = cos(theta)*VA + sin(theta)*VB
#    -> only 2 linear solves per candidate placement (not n_angles!)
# ---------------------------------------------------------------
def solve_gic(blocked_nodes):
    Rg = R_ground_base.copy()
    for b in blocked_nodes:
        Rg[b] = R_block
    Yg = np.diag(1.0 / Rg)
    M = Y + Yg

    VA, VB = np.linalg.solve(M, np.column_stack([A, B])).T
    V_all = np.outer(np.cos(thetas), VA) + np.outer(np.sin(thetas), VB)
    GIC_all = V_all * (1.0 / Rg)
    return GIC_all

# ---------------------------------------------------------------
# 5. Robust (worst-case) placement search
# ---------------------------------------------------------------
k_devices = 2
candidates = list(combinations(range(n_nodes), k_devices))

t0 = time.time()
baseline_GIC = solve_gic([])
worst_baseline = np.max(np.abs(baseline_GIC))

best_score = np.inf
best_combo = None
best_GIC = None

for combo in candidates:
    GIC_all = solve_gic(combo)
    score = np.max(np.abs(GIC_all))
    if score < best_score:
        best_score = score
        best_combo = combo
        best_GIC = GIC_all

elapsed = time.time() - t0

print(f"Number of candidate placements evaluated : {len(candidates)}")
print(f"Angular resolution                        : {n_angles} steps")
print(f"Total computation time                    : {elapsed:.4f} s")
print(f"Baseline worst-case GIC (no mitigation)    : {worst_baseline:.2f} A")
print(f"Best blocking-device placement (nodes)     : {best_combo}")
print(f"Worst-case GIC after mitigation            : {best_score:.2f} A")
print(f"Risk reduction                             : {(1 - best_score/worst_baseline)*100:.1f} %")

# ---------------------------------------------------------------
# 6. Visualization (single combined figure)
# ---------------------------------------------------------------
fig = plt.figure(figsize=(16, 12))

ax1 = fig.add_subplot(2, 2, 1, projection='3d')
Theta_deg = np.degrees(thetas)
Node_idx = np.arange(n_nodes)
TT, NN = np.meshgrid(Theta_deg, Node_idx, indexing='ij')
ax1.plot_surface(TT, NN, np.abs(baseline_GIC), cmap='inferno', edgecolor='none')
ax1.set_xlabel('Storm direction θ (deg)')
ax1.set_ylabel('Node index')
ax1.set_zlabel('|GIC| (A)')
ax1.set_title('Baseline: |GIC| vs storm direction (no mitigation)')

ax2 = fig.add_subplot(2, 2, 2, projection='3d')
ax2.plot_surface(TT, NN, np.abs(best_GIC), cmap='viridis', edgecolor='none')
ax2.set_xlabel('Storm direction θ (deg)')
ax2.set_ylabel('Node index')
ax2.set_zlabel('|GIC| (A)')
ax2.set_title(f'After mitigation at nodes {best_combo}')

ax3 = fig.add_subplot(2, 2, 3)
ax3.plot(Theta_deg, np.max(np.abs(baseline_GIC), axis=1), label='Baseline (no mitigation)', color='crimson')
ax3.plot(Theta_deg, np.max(np.abs(best_GIC), axis=1), label=f'Mitigated {best_combo}', color='seagreen')
ax3.axhline(best_score, color='seagreen', linestyle='--', linewidth=1, alpha=0.7)
ax3.set_xlabel('Storm direction θ (deg)')
ax3.set_ylabel('Max |GIC| across grid (A)')
ax3.set_title('Worst-node GIC vs storm direction')
ax3.legend()
ax3.grid(alpha=0.3)

ax4 = fig.add_subplot(2, 2, 4)
worst_case_vals = np.max(np.abs(best_GIC), axis=0)
for (i, j) in edges:
    ax4.plot([coords[i,0], coords[j,0]], [coords[i,1], coords[j,1]], color='gray', zorder=1, linewidth=1.2)
sizes = 300 + 40 * worst_case_vals
colors = ['red' if idx in best_combo else 'steelblue' for idx in range(n_nodes)]
ax4.scatter(coords[:,0], coords[:,1], s=sizes, c=colors, edgecolor='black', zorder=2)
for idx in range(n_nodes):
    ax4.annotate(f"{idx}", (coords[idx,0], coords[idx,1]), ha='center', va='center', fontsize=9, color='white', zorder=3)
ax4.set_xlabel('X (km)')
ax4.set_ylabel('Y (km)')
ax4.set_title('Grid topology (red = blocking device installed)')
ax4.set_aspect('equal')

plt.tight_layout()
plt.show()

4. Code Walkthrough

Section 1 — Topology: Eight substations are placed on a 2D plane (in kilometers) and connected as a ring with a central hub, which mirrors how many regional transmission networks are actually built for redundancy. Line resistance is computed from geographic distance using a realistic resistance-per-kilometer figure for extra-high-voltage lines.

Section 2 — Grounding parameters: Every transformer neutral starts with a small, realistic grounding resistance (0.6 Ω). A “blocked” node is simulated by replacing that resistance with an enormous value (1 MΩ), which forces near-zero current through that neutral — exactly what a physical blocking capacitor does for DC/quasi-DC current.

Section 3 — Admittance matrix construction: This builds the graph Laplacian-style matrix $\mathbf{Y}$ from line conductances, plus two auxiliary vectors $\mathbf{A}$ and $\mathbf{B}$ that capture how much current each node “collects” per unit of eastward or northward geoelectric field, respectively. These vectors are computed once and reused for every storm direction.

Section 4 — The fast solver: Instead of solving the linear system once per storm angle, we solve it exactly twice (np.linalg.solve with a 2-column right-hand side) to get $\mathbf{V}_A$ and $\mathbf{V}_B$. Every angle’s voltage vector is then reconstructed instantly via cos(θ)·V_A + sin(θ)·V_B, vectorized across all 360 angles simultaneously with np.outer.

Section 5 — Robust placement search: For each of the 28 possible 2-node placements, the code computes the worst-case (maximum absolute) GIC across all substations and all 360 storm directions, then keeps the placement with the smallest worst case — a classic minimax robust design.

Section 6 — Visualization: Builds one combined figure with four panels described below.

5. Why This Runs Fast (Speed Optimization Explained)

A naive implementation would loop over all 360 storm angles inside the loop over 28 placements, rebuilding and re-solving an 8×8 linear system 10,080 times. Because the system matrix $(\mathbf{Y}+\mathbf{Y}_e)$ does not change with storm angle — only the right-hand side does — we exploit the linearity shown in Section 1 of the physics discussion to solve the system only twice per placement (56 solves total instead of 10,080), then reconstruct all 360 angles algebraically. For this small 8-node network the difference is invisible in wall-clock time, but the technique matters enormously at realistic scale: a continent-scale grid model with hundreds of substations and thousands of candidate placements would make the naive approach computationally prohibitive, while this vectorized approach scales linearly with the number of placements rather than the product of placements and angular resolution.

6. Execution Results

Number of candidate placements evaluated : 28
Angular resolution                        : 360 steps
Total computation time                    : 0.0177 s
Baseline worst-case GIC (no mitigation)    : 346.79 A
Best blocking-device placement (nodes)     : (2, 7)
Worst-case GIC after mitigation            : 323.44 A
Risk reduction                             : 6.7 %

7. Interpreting the Graphs

Top-left (3D surface, baseline): This shows the magnitude of GIC at every substation (y-axis) as the storm direction sweeps through 360 degrees (x-axis), with no mitigation installed. The tall ridges reveal which substations are most exposed, and at which storm orientation — typically substations aligned with the dominant line direction see the highest peaks, since the induced EMF is strongest when the geoelectric field is parallel to the transmission line.

Top-right (3D surface, mitigated): The same surface after installing blocking devices at the optimal two substations. The peaks at the protected nodes are flattened to near zero, and — importantly — the peaks at the remaining nodes are also reduced, because blocking a node changes the current distribution throughout the whole network, not just at that node.

Bottom-left (worst-case curve): This directly compares, for every storm direction, the single highest GIC value anywhere in the network, before (red) and after (green) mitigation. The dashed horizontal line marks the guaranteed worst-case ceiling achieved by the optimal placement — this is the number a grid operator would actually report as “maximum expected GIC under any storm orientation.”

Bottom-right (network map): A geographic view of the grid, with node size proportional to its worst-case GIC after mitigation and the two selected blocking-device substations highlighted in red. This panel is the most operationally useful output: it tells an engineer exactly where to schedule equipment installation.

8. Takeaways for Real-World Grid Operators

This example uses a small, illustrative network, but the underlying method — the Lehtinen–Pirjola model combined with a linearity-exploiting solver and a robust minimax search over storm direction — is exactly the kind of approach used in academic and utility-industry studies of GIC mitigation. For larger real-world grids (hundreds of substations), exhaustive combination search becomes impractical once the number of devices $k$ grows, and the same linear-algebra speedup would typically be paired with a greedy or integer-programming heuristic rather than brute-force enumeration. Nonetheless, the physics, the math, and the vectorization trick shown here remain the foundation of that larger-scale analysis.

Telling Two Quantum Channels Apart

1. The diamond norm

2. The Choi matrix

3. The SDP (Watrous’ formulation)

4. The example problems

5. Source code

6. Code walkthrough

6.1 Channels and Choi matrices

6.2 The SDP in diamond_norm

6.3 Speed

6.4 Comparison tools

7. Results

Execution result image

Console output

8. Reading the figure

9. Takeaways

Maximizing the Holevo Information of the Amplitude Damping Channel

1. The Example Problem

2. Theory in Brief

Classical capacity

A symmetric ensemble

A natural but suboptimal choice

Quantum capacity

3. The Code

4. Execution Results

Console output

Graph output

5. Code Walkthrough

Section 1: Basic tools

Section 2: Slow and fast Holevo information

Section 3: General numerical maximization

Section 4: The symmetric ring ensemble

Section 5: Lower bound and quantum capacity

Section 6: The table

Section 7: Full curves

Section 8: The figure

6. Reading the Results

Panel (a): the three capacities

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

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

Panel (d): the optimal ensemble

Panel (e): the geometry (3D)

Panel (f): sanity check

7. Conclusion

Finding the Channel That Maximizes Fidelity with Python

1. The Problem

Known analytic answer (our benchmark)

2. Turning the Problem into Linear Algebra

Stinespring form

Replacing the Haar integral by six states

A fast fixed-point algorithm

3. Reference Strategies

4. Full Source Code

5. Code Walkthrough

Section 1 – Dimensions

Section 2 – The octahedron design

Section 3 – Core routines

Section 4 – The CNOT cloner

Section 5 – The symmetric optimum

Section 6 – Asymmetric trade-off

Section 7 – Visualization

6. Execution Results

7. Reading the Results

8. Summary

Maximizing Mutual Information and Accessible Information with Python

1. Problem Setting

2. The Example: the Trine Ensemble

A closed-form candidate

3. Strategy for the Numerical Optimization

Parametrizing a valid POVM

Bloch representation of a POVM element

A two-parameter family for a 3D landscape

4. The Complete Source Code

5. Code Walkthrough

Section 1: Ensemble construction

Section 2: Mutual information

Section 3: POVM parametrization and optimization

Section 4: Merging duplicate outcomes

Section 5: Reference measurements and the landscapes

Section 6: Console output

6.2 The SDP in `diamond_norm`

4.2 The SDP solver `solve_msd`