Optimal Water Allocation Over 20 Turns

October 22, 2026

Balancing Agriculture, Households, and Industry

Water is the one resource that every sector needs and no sector can replace. Farms need it to feed people, households need it to keep people alive, and factories turn it into the tax revenue that pays for tomorrow’s infrastructure. When a drought arrives, every cubic meter given to one sector is taken from another.

In this article we formulate a concrete water allocation problem as a 20-turn dynamic optimization, solve it numerically with Python, and visualize the result with a single combined figure that includes two 3D plots.

1. Problem Statement

A region receives a limited amount of water $W_t$ in each turn $t = 0, 1, \dots, 19$. The planner splits it among three sectors:

$$
\alpha_t^{\mathrm{agr}} + \alpha_t^{\mathrm{dom}} + \alpha_t^{\mathrm{ind}} = 1, \qquad \alpha_t^{\bullet} \ge 0 .
$$

The goal is to maximize the sum of the final population and the final asset value:

$$
\max_ ; J = P_{20} + w , A_{20}, \qquad w = 1 .
$$

Here $P$ is the population in thousands of people and $A$ is the asset value in millions of USD.

2. Mathematical Model

2.1 Water supply

The base supply drops sharply in a drought centered on turn 9. Infrastructure built from past tax revenue reduces leakage and raises the effective supply:

$$
W_t^{\mathrm{base}} = W_0\left(1 - 0.55, e^{-\left(\frac{t-9}{3}\right)^2}\right), \qquad
W_t = W_t^{\mathrm{base}}\left(1 + \eta_{\max}\frac{A_t}{A_t + A_{1/2}}\right).
$$

2.2 Satisfaction ratios

The water delivered to each sector is $W_t^{\mathrm{agr}} = \alpha_t^{\mathrm{agr}} W_t$, and similarly for the other two. The demand is proportional to the population, with per-capita coefficients $q_a$ and $q_d$:

$$
r_t^{\mathrm{agr}} = \frac{W_t^{\mathrm{agr}}}{q_a P_t}, \qquad
r_t^{\mathrm{dom}} = \frac{W_t^{\mathrm{dom}}}{q_d P_t}.
$$

The satisfaction level is a smooth version of $\min(1, r)$, which keeps the objective differentiable for gradient-based optimization:

$$
s(r) = 1 - \frac{1}{\beta}\log!\left(1 + e^{\beta (1-r)}\right).
$$

2.3 Population dynamics

Water shortage raises mortality, and a food shortage also suppresses births. The logistic factor represents the carrying capacity $K$:

$$
m_t = m_0 + m_d,(1 - s_t^{\mathrm{dom}})^2 + m_a,(1 - s_t^{\mathrm{agr}})^2,
$$

$$
P_{t+1} = P_t\left[1 + b, s_t^{\mathrm{agr}} s_t^{\mathrm{dom}}\left(1 - \frac{P_t}{K}\right) - m_t\right].
$$

2.4 Industry, tax revenue, and assets

Industrial output saturates with industrial water and scales with the square root of the workforce. A fraction $\kappa$ of the output is collected as tax and reinvested, while assets depreciate at rate $\delta$:

$$
Y_t = Y_{\max}\left(1 - e^{-W_t^{\mathrm{ind}}/k_i}\right)\sqrt{\frac{P_t}{P_0}}, \qquad
A_{t+1} = (1-\delta)A_t + \kappa Y_t .
$$

This creates the central trade-off. Industrial water increases future assets and future water efficiency, while household and agricultural water keeps people alive and growing.

2.5 Concrete example

Symbol	Value	Meaning
$P_0$, $A_0$	100, 200	Initial population (thousand) and asset value (million USD)
$W_0$	85	Base water supply per turn
$q_a$, $q_d$	0.50, 0.30	Per-capita agricultural and domestic demand
$K$, $b$	260, 0.09	Carrying capacity and birth coefficient
$m_0, m_d, m_a$	0.010, 0.25, 0.20	Mortality coefficients
$Y_{\max}$, $k_i$	60, 25	Industrial output scale and saturation constant
$\kappa$, $\delta$	0.5, 0.03	Reinvestment rate and depreciation
$\eta_{\max}$, $A_{1/2}$	0.4, 400	Infrastructure efficiency gain

3. Solution Approach

We compare five policies:

Equal split, with one third of the water to each sector.
Agriculture first, with shares (0.60, 0.30, 0.10).
Industry first, with shares (0.25, 0.20, 0.55).
Best constant, found by a grid search over all fixed allocations.
Optimized (dynamic), where the allocation can change every turn.

The dynamic policy has $20 \times 3 = 60$ logits $z_t$, mapped to valid shares by a softmax:

$$
\alpha_t^{k} = \frac{e^{z_t^{k}}}{\sum_{j} e^{z_t^{j}}} .
$$

The softmax removes the simplex constraint, so we can use L-BFGS-B directly. Because the objective is non-convex, we run six starting points and keep the best result. The simulator uses plain scalar arithmetic inside the 20-step loop, so the full optimization finishes in about a second and no further acceleration is needed.

4. Source Code

import math
import time
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
from scipy.optimize import minimize

plt.style.use("dark_background")

# ----------------------------------------------------------------------
# 1. Model parameters
# ----------------------------------------------------------------------
T = 20
P0, A0 = 100.0, 200.0
W0 = 85.0
Q_A, Q_D = 0.50, 0.30
K_POP = 260.0
BIRTH = 0.09
M0, M_D, M_A = 0.010, 0.25, 0.20
Y_MAX, K_I = 60.0, 25.0
KAPPA, DELTA = 0.5, 0.03
EFF_MAX, EFF_HALF = 0.4, 400.0
WEIGHT = 1.0
BETA = 40.0

def base_supply(t):
    return W0 * (1.0 - 0.55 * math.exp(-((t - 9.0) / 3.0) ** 2))

def smooth_sat(r):
    x = 1.0 - r
    sp = np.logaddexp(0.0, BETA * x) / BETA
    return float(min(1.0, max(0.0, 1.0 - sp)))

def simulate(alloc):
    P, A = P0, A0
    Ps, As, Sa, Sd, Ws = [P], [A], [], [], []
    for t in range(T):
        eff = 1.0 + EFF_MAX * A / (A + EFF_HALF)
        W = base_supply(t) * eff
        wa, wd, wi = alloc[t, 0] * W, alloc[t, 1] * W, alloc[t, 2] * W
        s_a = smooth_sat(wa / (Q_A * P))
        s_d = smooth_sat(wd / (Q_D * P))
        mort = M0 + M_D * (1.0 - s_d) ** 2 + M_A * (1.0 - s_a) ** 2
        growth = BIRTH * s_a * s_d * max(0.0, 1.0 - P / K_POP)
        P_new = max(P * (1.0 + growth - mort), 1.0)
        out = Y_MAX * (1.0 - math.exp(-wi / K_I)) * math.sqrt(P / P0)
        A_new = (1.0 - DELTA) * A + KAPPA * out
        P, A = P_new, A_new
        Ps.append(P); As.append(A); Sa.append(s_a); Sd.append(s_d); Ws.append(W)
    J = P + WEIGHT * A
    return J, np.array(Ps), np.array(As), np.array(Sa), np.array(Sd), np.array(Ws)

def constant_policy(a, d, i):
    return np.tile(np.array([a, d, i], dtype=float), (T, 1))

def softmax_rows(z):
    z = z - z.max(axis=1, keepdims=True)
    e = np.exp(z)
    return e / e.sum(axis=1, keepdims=True)

def neg_objective(zflat):
    alloc = softmax_rows(zflat.reshape(T, 3))
    return -simulate(alloc)[0]

# ----------------------------------------------------------------------
# 2. Baseline policies and best constant policy
# ----------------------------------------------------------------------
policies = {
    "Equal split": constant_policy(1 / 3, 1 / 3, 1 / 3),
    "Agriculture first": constant_policy(0.60, 0.30, 0.10),
    "Industry first": constant_policy(0.25, 0.20, 0.55),
}

step = 0.025
grid = np.arange(0.0, 1.0 + 1e-9, step)
gx, gy, gz = [], [], []
best_const, best_const_J = None, -1e18
for a in grid:
    for i in grid:
        d = 1.0 - a - i
        if d < -1e-9:
            continue
        d = max(d, 0.0)
        J = simulate(constant_policy(a, d, i))[0]
        gx.append(a); gy.append(i); gz.append(J)
        if J > best_const_J:
            best_const_J, best_const = J, (a, d, i)
policies["Best constant"] = constant_policy(*best_const)

# ----------------------------------------------------------------------
# 3. Dynamic optimization (multi-start L-BFGS-B)
# ----------------------------------------------------------------------
t0 = time.time()
rng = np.random.default_rng(42)
starts = [np.zeros(T * 3)]
bc = np.log(np.maximum(np.array(best_const), 1e-3))
starts.append(np.tile(bc, T))
for _ in range(4):
    starts.append(rng.normal(0.0, 0.7, T * 3))

best_res = None
for z0 in starts:
    res = minimize(neg_objective, z0, method="L-BFGS-B",
                   bounds=[(-8.0, 8.0)] * (T * 3),
                   options={"maxiter": 300})
    if best_res is None or res.fun < best_res.fun:
        best_res = res
elapsed = time.time() - t0
policies["Optimized (dynamic)"] = softmax_rows(best_res.x.reshape(T, 3))

# ----------------------------------------------------------------------
# 4. Evaluation
# ----------------------------------------------------------------------
results = {name: simulate(al) for name, al in policies.items()}
print("Optimization time: %.2f s" % elapsed)
print("Best constant allocation (agri, domestic, industry): "
      "(%.3f, %.3f, %.3f)" % best_const)
print()
print("%-22s %10s %12s %10s" % ("Policy", "Population", "Asset value", "J"))
for name, r in results.items():
    print("%-22s %10.2f %12.2f %10.2f" % (name, r[1][-1], r[2][-1], r[0]))
print()
opt_alloc = policies["Optimized (dynamic)"]
print("Optimized allocation by turn (agri, domestic, industry):")
for t in range(T):
    print("turn %2d  W=%6.1f  (%.3f, %.3f, %.3f)" %
          (t + 1, results["Optimized (dynamic)"][5][t],
           opt_alloc[t, 0], opt_alloc[t, 1], opt_alloc[t, 2]))

# ----------------------------------------------------------------------
# 5. Visualization
# ----------------------------------------------------------------------
colors = {
    "Equal split": "#4cc9f0",
    "Agriculture first": "#90be6d",
    "Industry first": "#f9c74f",
    "Best constant": "#b388eb",
    "Optimized (dynamic)": "#ff4d6d",
}
turns = np.arange(1, T + 1)
fig = plt.figure(figsize=(20, 12))
gs = GridSpec(2, 3, figure=fig)

ax1 = fig.add_subplot(gs[0, 0])
ax1.stackplot(turns, opt_alloc[:, 0], opt_alloc[:, 1], opt_alloc[:, 2],
              labels=["Agriculture", "Domestic", "Industry"],
              colors=["#90be6d", "#4cc9f0", "#f9c74f"], alpha=0.9)
ax1.set_title("Optimized water allocation share")
ax1.set_xlabel("Turn"); ax1.set_ylabel("Share of available water")
ax1.set_xlim(1, T); ax1.set_ylim(0, 1)
ax1b = ax1.twinx()
ax1b.plot(turns, results["Optimized (dynamic)"][5], color="white", lw=2, ls="--",
          label="Effective supply")
ax1b.set_ylabel("Effective supply")
h1, l1 = ax1.get_legend_handles_labels()
h2, l2 = ax1b.get_legend_handles_labels()
ax1.legend(h1 + h2, l1 + l2, loc="lower left", fontsize=8)

ax2 = fig.add_subplot(gs[0, 1])
for name, r in results.items():
    ax2.plot(np.arange(0, T + 1), r[1], lw=2.5 if "Optimized" in name else 1.6,
             color=colors[name], label=name)
ax2.set_title("Population trajectory")
ax2.set_xlabel("Turn"); ax2.set_ylabel("Population (thousand)")
ax2.legend(fontsize=8); ax2.grid(alpha=0.2)

ax3 = fig.add_subplot(gs[0, 2])
for name, r in results.items():
    ax3.plot(np.arange(0, T + 1), r[2], lw=2.5 if "Optimized" in name else 1.6,
             color=colors[name], label=name)
ax3.set_title("Asset value trajectory")
ax3.set_xlabel("Turn"); ax3.set_ylabel("Asset value (million USD)")
ax3.legend(fontsize=8); ax3.grid(alpha=0.2)

ax4 = fig.add_subplot(gs[1, 0], projection="3d")
surf = ax4.plot_trisurf(np.array(gx), np.array(gy), np.array(gz),
                        cmap="plasma", linewidth=0.1, antialiased=True, alpha=0.95)
ax4.scatter([best_const[0]], [best_const[2]], [best_const_J], color="white",
            s=60, label="Best constant")
ax4.scatter([best_const[0]], [best_const[2]], [results["Optimized (dynamic)"][0]],
            color="#ff4d6d", s=80, marker="*", label="Optimized (dynamic)")
ax4.set_title("Objective J vs constant allocation")
ax4.set_xlabel("Agriculture share"); ax4.set_ylabel("Industry share")
ax4.set_zlabel("J")
ax4.view_init(elev=28, azim=-125)
ax4.legend(fontsize=8)
fig.colorbar(surf, ax=ax4, shrink=0.55, pad=0.08, label="J")

ax5 = fig.add_subplot(gs[1, 1], projection="3d")
for name, r in results.items():
    ax5.plot(r[1], r[2], np.arange(0, T + 1), lw=2.5 if "Optimized" in name else 1.5,
             color=colors[name], label=name)
    ax5.scatter([r[1][-1]], [r[2][-1]], [T], color=colors[name], s=30)
ax5.set_title("State-space trajectory")
ax5.set_xlabel("Population (thousand)"); ax5.set_ylabel("Asset value")
ax5.set_zlabel("Turn")
ax5.view_init(elev=22, azim=-60)
ax5.legend(fontsize=8)

ax6 = fig.add_subplot(gs[1, 2])
names = list(results.keys())
pops = [results[n][1][-1] for n in names]
assets = [WEIGHT * results[n][2][-1] for n in names]
xpos = np.arange(len(names))
ax6.bar(xpos, pops, color="#4cc9f0", label="Population")
ax6.bar(xpos, assets, bottom=pops, color="#f9c74f", label="Asset value")
for k, n in enumerate(names):
    ax6.text(k, pops[k] + assets[k] + 3, "%.1f" % results[n][0],
             ha="center", fontsize=9)
ax6.set_xticks(xpos)
ax6.set_xticklabels([n.replace(" ", "\n") for n in names], fontsize=8)
ax6.set_ylim(0, max(p + a for p, a in zip(pops, assets)) * 1.12)
ax6.set_title("Final objective J = population + asset value")
ax6.set_ylabel("J"); ax6.legend(fontsize=8); ax6.grid(alpha=0.2, axis="y")

fig.suptitle("Water allocation optimization over 20 turns", fontsize=18)
fig.tight_layout(rect=[0, 0, 1, 0.96])
plt.show()

5. Code Walkthrough

5.1 Parameters and helper functions

The constants in the first block correspond one-to-one with the table in Section 2.5. base_supply implements the Gaussian-shaped drought, which cuts the base supply by up to 55% around turn 9.

smooth_sat implements the smoothed satisfaction function $s(r)$. It uses np.logaddexp(0, x), which computes $\log(1+e^{x})$ without numerical overflow. The result is clipped to $[0, 1]$ to remove tiny floating-point excursions.

5.2 The simulator

simulate runs the 20-turn dynamics for any allocation matrix of shape (T, 3). In each turn it:

Computes the effective supply, including the infrastructure bonus from current assets.
Splits the supply into agricultural, domestic, and industrial water.
Converts each delivery into a satisfaction level.
Updates the population from the mortality and growth formulas.
Updates the assets from the industrial output.

It returns the objective $J$ together with all time series, which are reused later for plotting.

5.3 Softmax parametrization

softmax_rows converts a $20 \times 3$ matrix of logits into valid shares. Every row is non-negative and sums to one, so the optimizer works in an unconstrained space. Subtracting the row maximum before exponentiation is a standard stabilization step.

5.4 Baselines and the constant-policy grid

The three hand-made policies are the reference points. The double loop then scans every fixed allocation $(\alpha^{\mathrm{agr}}, \alpha^{\mathrm{ind}})$ on a grid of step 0.025, with $\alpha^{\mathrm{dom}} = 1 - \alpha^{\mathrm{agr}} - \alpha^{\mathrm{ind}}$. These roughly 800 evaluations serve two purposes. They give the best constant policy, and they provide the data for the 3D surface.

5.5 Dynamic optimization

The decision vector has 60 entries. We run L-BFGS-B from six starting points:

the uniform allocation,
the best constant allocation,
four random perturbations with a fixed seed.

The smooth satisfaction function matters here. A hard min would produce kinks that confuse gradient-based search, while the softplus version has a well-defined gradient everywhere. Bounds of $\pm 8$ on the logits keep the softmax numerically tame.

5.6 Evaluation and printing

Every policy is re-simulated, and the script prints a comparison table followed by the optimized allocation for each turn. The console output is the numerical counterpart of the figure.

5.7 The combined figure

All plots are drawn in a single $2 \times 3$ figure:

Top left: the optimized allocation as a stacked area, with the effective supply as a dashed line.
Top middle and right: population and asset trajectories for all five policies.
Bottom left (3D): the objective $J$ as a surface over constant allocations, with the best constant policy and the dynamic optimum marked.
Bottom middle (3D): the trajectory of each policy in (population, assets, turn) space.
Bottom right: a stacked bar chart that splits the final $J$ into its population and asset components.

6. Results

Combined figure

Console output

Optimization time: 4.39 s
Best constant allocation (agri, domestic, industry): (0.350, 0.250, 0.400)

Policy                 Population  Asset value          J
Equal split                 90.11       390.24     480.36
Agriculture first          115.21       238.17     353.38
Industry first              68.17       415.54     483.71
Best constant               89.24       410.27     499.50
Optimized (dynamic)         90.79       414.54     505.33

Optimized allocation by turn (agri, domestic, industry):
turn  1  W=  96.3  (0.308, 0.254, 0.438)
turn  2  W=  97.0  (0.309, 0.250, 0.441)
turn  3  W=  97.4  (0.309, 0.246, 0.445)
turn  4  W=  97.2  (0.307, 0.244, 0.449)
turn  5  W=  95.3  (0.303, 0.244, 0.454)
turn  6  W=  89.9  (0.295, 0.247, 0.458)
turn  7  W=  79.4  (0.282, 0.258, 0.460)
turn  8  W=  64.7  (0.267, 0.280, 0.453)
turn  9  W=  50.9  (0.272, 0.313, 0.415)
turn 10  W=  45.1  (0.324, 0.333, 0.342)
turn 11  W=  51.0  (0.375, 0.310, 0.314)
turn 12  W=  65.0  (0.379, 0.268, 0.353)
turn 13  W=  80.2  (0.362, 0.235, 0.403)
turn 14  W=  91.3  (0.348, 0.219, 0.433)
turn 15  W=  97.5  (0.344, 0.216, 0.440)
turn 16  W= 100.1  (0.347, 0.219, 0.433)
turn 17  W= 101.2  (0.353, 0.225, 0.422)
turn 18  W= 101.6  (0.358, 0.232, 0.410)
turn 19  W= 101.9  (0.362, 0.239, 0.399)
turn 20  W= 102.1  (0.367, 0.245, 0.388)

7. Discussion

7.1 The one-sided policies fail in opposite ways

The agriculture-first policy keeps the population large, since farmers and households are almost never short of water. But it starves industry, so the asset value barely moves from its starting level and the final $J$ is by far the lowest.

The industry-first policy shows the opposite failure. Assets grow quickly, but during the drought the population falls to its lowest level. Because industrial output scales with $\sqrt{P_t}$, the shrinking workforce eventually drags on the very thing the policy was trying to maximize.

7.2 The 3D surface has a single broad peak

The surface of $J$ over constant allocations is smooth and unimodal. Its maximum lies at roughly 35% agriculture, 25% domestic use, and 40% industry. The surface falls off steeply when the agricultural share drops very low, because that region produces famine-level mortality, and also when the industrial share approaches zero, because the asset engine stalls. The practical lesson is that a moderately balanced allocation already captures most of the achievable value.

7.3 Dynamic reallocation adds a further gain

The optimized policy beats the best constant allocation. The stacked area chart shows why:

Before the drought, the optimizer gives industry its largest share. Early assets compound through reinvestment and also raise the effective water supply, so water spent on industry pays back twice.
During the drought, industry is cut back while the domestic and agricultural shares rise. The effective supply is at its lowest, and protecting the population from the quadratic mortality penalty is worth more than the lost output.
After the drought, the industrial share recovers and the allocation drifts back toward its pre-drought pattern.

In other words, the optimizer treats industry as an investment to be front-loaded and treats household and farm water as insurance to be strengthened exactly when the supply is scarce.

7.4 State-space view

In the 3D trajectory plot, every policy starts at the same point. The agriculture-first path moves toward high population and low assets, and the industry-first path moves toward high assets and low population. The optimized path stays close to the efficient frontier of the other policies. It accepts a temporary population dip during the drought but ends with the highest combined value.

7.5 Final composition of the objective

The stacked bars make the trade-off explicit. Populations and asset values sit on different scales, so the weight $w$ matters: increasing $w$ shifts the optimum toward industry, and decreasing it shifts the optimum toward households and agriculture. Because the code exposes WEIGHT as a single constant, the policy preference can be explored by changing one line.

8. Conclusion

We modeled the allocation of scarce water as a 20-turn dynamic optimization problem with three coupled effects: shortage-driven mortality, industrial tax revenue that compounds into assets, and infrastructure that improves the effective supply. A smooth satisfaction function made the problem suitable for gradient-based optimization, and a softmax parametrization removed the allocation constraints.

The main findings were:

Extreme policies are dominated. Neglecting either industry or agriculture sharply lowers the final objective.
A well-chosen constant split already performs well, and the objective surface around it is smooth and forgiving.
A time-varying policy performs better still. Invest in industry while water is plentiful, then shift water to people and food during the drought.

The same framework extends naturally to stochastic droughts, multiple regions sharing a river, or water-saving technology investments as an additional decision variable.

Optimizing Farmland vs. Housing with Python

October 21, 2026

How Many People Can a Village Feed and House?

A village of 100 people owns a fixed amount of land. Every hectare turned into a field produces food, and every hectare turned into housing raises the population cap. Ten years from now, how many people can the village support, and what split of land gets it there?

This is a small problem, but it shows how two limiting factors, food and housing, pull against each other. In this post I model the village, solve it by brute-force search, and plot the result in 2D and 3D.

1. The Problem

The village has:

Item	Value
Total land $L$	60 ha
Initial population $P_0$	100 people
Planning horizon	10 years
Farming productivity $A$	8.0
Diminishing-returns exponent $\alpha$	0.85
Housing density $D$	10 people/ha
Maximum growth rate $r$	0.35 per year
Starvation severity $m$	0.60
Food needed per person $c$	1.0 per year

The decision variable is $x$, the hectares of farmland. The remaining $L - x$ hectares become residential land. The goal is to maximize the population after 10 years.

2. The Model

Food production shows diminishing returns, since extra fields need more labor, are less fertile and so on:

$$
F(x) = A,x^{\alpha}
$$

Housing capacity is linear in residential land:

$$
H(x) = D,(L - x)
$$

Each year the population grows logistically toward the housing capacity, and it shrinks if food runs short. The food sufficiency ratio is

$$
q_t = \frac{F(x)}{c,P_t}, \qquad s_t = \max\bigl(0,\ 1 - q_t\bigr)
$$

and the population update is

$$
P_{t+1} = P_t + r,P_t\left(1 - \frac{P_t}{H(x)}\right) - m,P_t,s_t
$$

The optimization problem is

$$
\max_{0 \le x \le L}\ P_{10}(x)
$$

Why is there an interior optimum? Too little farmland starves the village. Too little housing caps the population. At a long-run equilibrium with food shortage, the growth term and the starvation term balance:

so the best split balances the two constraints rather than maximizing either one.

3. The Complete Source Code

Everything is in one script. The search evaluates thousands of candidate land splits simultaneously with NumPy array operations instead of a Python loop over candidates.

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

# ------------------------------------------------------------
# 1. Parameters of the village
# ------------------------------------------------------------
L = 60.0        # total usable land [ha]
P0 = 100.0      # initial population
YEARS = 10      # planning horizon
A = 8.0         # farming productivity coefficient
ALPHA = 0.85    # diminishing returns on farmland (0 < ALPHA <= 1)
D = 10.0        # people housed per hectare of residential land
R = 0.35        # maximum annual growth rate
M = 0.60        # starvation severity
C = 1.0         # food needed per person per year


# ------------------------------------------------------------
# 2. Vectorised simulator (x, A, D may be scalars or arrays)
# ------------------------------------------------------------
def simulate(x, A=A, D=D):
    """Return the population history, shape (YEARS + 1, *x.shape)."""
    x, A_, D_ = np.broadcast_arrays(
        np.asarray(x, dtype=float),
        np.asarray(A, dtype=float),
        np.asarray(D, dtype=float),
    )
    food = A_ * x ** ALPHA                       # annual food production
    cap = np.maximum(D_ * (L - x), 1e-9)         # housing capacity
    P = np.full(x.shape, P0)
    hist = [P.copy()]
    for _ in range(YEARS):
        q = food / (C * np.maximum(P, 1e-9))     # food sufficiency ratio
        shortage = np.clip(1.0 - q, 0.0, 1.0)
        P = P + R * P * (1.0 - P / cap) - M * P * shortage
        P = np.maximum(P, 0.0)
        hist.append(P.copy())
    return np.array(hist)


# ------------------------------------------------------------
# 3. Optimisation by exhaustive search on a fine grid
# ------------------------------------------------------------
x_grid = np.linspace(0.0, L, 6001)               # step = 0.01 ha
final_pop = simulate(x_grid)[-1]
i_best = int(np.argmax(final_pop))
x_best, p_best = x_grid[i_best], final_pop[i_best]

print("=" * 60)
print("Village land-allocation problem")
print("=" * 60)
print(f"Total land L = {L:.0f} ha, initial population = {P0:.0f}, horizon = {YEARS} years")
print(f"Optimal farmland   : {x_best:.2f} ha ({x_best / L:.1%} of the land)")
print(f"Optimal housing    : {L - x_best:.2f} ha")
print(f"Population in year {YEARS}: {p_best:.2f}")
print("-" * 60)
print(f"{'Farmland [ha]':>14} {'Housing [ha]':>13} {'Pop. (year 10)':>15}")
for xv in [10, 15, 20, 25, 30, 35, 40, 45, 50]:
    pv = simulate(float(xv))[-1]
    print(f"{xv:>14.0f} {L - xv:>13.0f} {float(pv):>15.2f}")
print("-" * 60)

best_hist = simulate(x_best)
food_best = A * x_best ** ALPHA
cap_best = D * (L - x_best)
print(f"Food production at optimum : {food_best:.2f} (feeds {food_best / C:.1f} people)")
print(f"Housing capacity at optimum: {cap_best:.2f} people")
print("Year-by-year population at the optimum:")
for t, p in enumerate(best_hist):
    print(f"  year {t:>2}: {float(p):7.2f}")

# ------------------------------------------------------------
# 4. Data for the figures
# ------------------------------------------------------------
years = np.arange(YEARS + 1)
sample_x = [10, 20, 30, float(np.round(x_best, 2)), 45, 50]
x_3d = np.linspace(0.0, L - 1.0, 70)
T3, X3 = np.meshgrid(years, x_3d)
Z3 = simulate(x_3d)[:, :].T                       # shape (len(x_3d), YEARS + 1)

a_vals = np.linspace(4.0, 12.0, 60)
XX, AA = np.meshgrid(np.linspace(0.0, L - 1.0, 60), a_vals)
ZZ = simulate(XX, A=AA)[-1]
best_x_for_a = XX[0][np.argmax(ZZ, axis=1)]

# ------------------------------------------------------------
# 5. One combined figure (2 x 3 panels)
# ------------------------------------------------------------
fig = plt.figure(figsize=(21, 12))

# (a) final population vs farmland
ax1 = fig.add_subplot(2, 3, 1)
ax1.plot(x_grid, final_pop, color="tab:green", lw=2)
ax1.axvline(x_best, color="red", ls="--", label=f"optimum: {x_best:.2f} ha")
ax1.scatter([x_best], [p_best], color="red", zorder=5)
ax1.set_xlabel("Farmland [ha]")
ax1.set_ylabel("Population in year 10")
ax1.set_title("(a) Final population vs. farmland")
ax1.grid(alpha=0.3)
ax1.legend()

# (b) trajectories
ax2 = fig.add_subplot(2, 3, 2)
for xv in sample_x:
    ax2.plot(years, simulate(float(xv)), marker="o", ms=4, label=f"farmland = {xv:g} ha")
ax2.set_xlabel("Year")
ax2.set_ylabel("Population")
ax2.set_title("(b) Population trajectories")
ax2.grid(alpha=0.3)
ax2.legend(fontsize=8)

# (c) the two limiting factors at year 0 conditions
ax3 = fig.add_subplot(2, 3, 3)
food_curve = A * x_grid ** ALPHA / C
house_curve = D * (L - x_grid)
ax3.plot(x_grid, food_curve, color="tab:orange", lw=2, label="People fed by farmland")
ax3.plot(x_grid, house_curve, color="tab:blue", lw=2, label="People housed by residential land")
ax3.plot(x_grid, np.minimum(food_curve, house_curve), color="black", ls=":", lw=2,
         label="Bottleneck = min(food, housing)")
ax3.axvline(x_best, color="red", ls="--")
ax3.set_xlabel("Farmland [ha]")
ax3.set_ylabel("People")
ax3.set_title("(c) Food limit vs. housing limit")
ax3.grid(alpha=0.3)
ax3.legend(fontsize=8)

# (d) 3D surface: time x farmland -> population
ax4 = fig.add_subplot(2, 3, 4, projection="3d")
ax4.plot_surface(T3, X3, Z3, cmap=cm.viridis, edgecolor="none", alpha=0.95)
ax4.set_xlabel("Year")
ax4.set_ylabel("Farmland [ha]")
ax4.set_zlabel("Population")
ax4.set_title("(d) 3D: population over time and farmland")
ax4.view_init(elev=28, azim=-130)

# (e) 3D surface: farmland x productivity -> final population
ax5 = fig.add_subplot(2, 3, 5, projection="3d")
ax5.plot_surface(XX, AA, ZZ, cmap=cm.plasma, edgecolor="none", alpha=0.95)
ax5.set_xlabel("Farmland [ha]")
ax5.set_ylabel("Productivity A")
ax5.set_zlabel("Population in year 10")
ax5.set_title("(e) 3D: final population vs. farmland and productivity")
ax5.view_init(elev=28, azim=-125)

# (f) heat map with the optimal ridge
ax6 = fig.add_subplot(2, 3, 6)
cs = ax6.contourf(XX, AA, ZZ, levels=25, cmap="plasma")
fig.colorbar(cs, ax=ax6, label="Population in year 10")
ax6.plot(best_x_for_a, a_vals, color="white", lw=2.5, label="Optimal farmland")
ax6.set_xlabel("Farmland [ha]")
ax6.set_ylabel("Productivity A")
ax6.set_title("(f) Optimal farmland as productivity changes")
ax6.legend(loc="upper right")

plt.suptitle("Farmland vs. Housing: Maximizing the Village Population After 10 Years",
             fontsize=17, y=1.0)
plt.tight_layout()
plt.show()

4. Code Walkthrough

Section 1: Parameters

All constants of the model are defined at the top, so you can experiment by changing a single number. ALPHA < 1 encodes diminishing returns on farmland. M controls how brutally the population shrinks when food runs short.

Section 2: The vectorised simulator

simulate(x, A, D) implements the update rule

$$
P_{t+1} = P_t + r,P_t\left(1 - \frac{P_t}{H}\right) - m,P_t,s_t
$$

Key details:

np.broadcast_arrays gives x, A and D a common shape. You can pass a single farmland value, an array of 6,001 candidates, or a 2D grid of farmland × productivity, and the same function handles all of them.
np.maximum(..., 1e-9) guards against division by zero. When farmland is 0 hectares, food is 0, and when it is 60 hectares, housing is 0. Without the guard you would get NaN values or warnings.
np.clip(1.0 - q, 0.0, 1.0) implements $s_t = \max(0, 1-q_t)$ and also caps the shortage at 100%, so the population cannot drop by more than $m$ in a year.
P = np.maximum(P, 0.0) prevents negative populations.
The loop runs only 10 times (one per year). The “loop” over candidate land splits is handled by NumPy internally.

Section 3: Optimisation

The search space is one-dimensional, so exhaustive search is both simple and reliable. np.linspace(0, L, 6001) creates 6,001 candidates with 0.01 ha resolution. A single call to simulate evaluates them all at once, and np.argmax picks the best. A scalar optimizer could be fooled by a flat plateau or an unusual shape, while a dense grid cannot.

The code then prints a summary, a table of selected land splits, and the year-by-year trajectory at the optimum. It also prints the food production and housing capacity at the optimum, which shows which constraint is binding.

Section 4: Data for the figures

sample_x contains the farmland values whose trajectories will be compared in panel (b).
T3, X3, Z3 build the grid for the first 3D surface. simulate(x_3d) returns an array of shape (11, 70) (years × candidates), so it is transposed to (70, 11) to match the meshgrid.
XX, AA, ZZ build the second surface by evaluating the final population over a grid of farmland and productivity A. This works only because simulate accepts array-valued A.
best_x_for_a finds, for each productivity level, the farmland that maximizes the final population, which gives the “optimal ridge” for the heat map.

Section 5: One combined figure

All six panels are drawn into a single fig with a 2 × 3 layout, so one plt.show() displays everything:

(a) the objective function $P_{10}(x)$ with the optimum marked,
(b) population trajectories for different land splits,
(c) the food limit versus the housing limit,
(d) a 3D surface of population over time and farmland,
(e) a 3D surface of final population over farmland and productivity,
(f) a heat map of (e) with the optimal ridge overlaid.

5. Console Output

============================================================
Village land-allocation problem
============================================================
Total land L = 60 ha, initial population = 100, horizon = 10 years
Optimal farmland   : 36.81 ha (61.4% of the land)
Optimal housing    : 23.19 ha
Population in year 10: 191.08
------------------------------------------------------------
 Farmland [ha]  Housing [ha]  Pop. (year 10)
            10            50          105.00
            15            45          135.00
            20            40          157.78
            25            35          174.43
            30            30          185.40
            35            25          190.66
            40            20          189.70
            45            15          148.77
            50            10          100.00
------------------------------------------------------------
Food production at optimum : 171.46 (feeds 171.5 people)
Housing capacity at optimum: 231.90 people
Year-by-year population at the optimum:
  year  0:  100.00
  year  1:  119.91
  year  2:  140.17
  year  3:  159.58
  year  4:  177.00
  year  5:  188.34
  year  6:  190.59
  year  7:  191.00
  year  8:  191.06
  year  9:  191.08
  year 10:  191.08

Reading the numbers

The search finds that the best allocation is about 36.8 hectares of farmland (61.4% of the land) and 23.2 hectares of housing, giving a population of about 191 people after 10 years, nearly double the starting 100.

The comparison table shows the shape of the trade-off:

With 10 ha of farmland, food is the bottleneck and the population barely moves (about 105).
With 50 ha of farmland, only 10 ha remain for housing. The cap of $10 \times 10 = 100$ people is exactly the starting population, so the village never grows (100.00).
The best results sit in a narrow band around 35–40 ha, where the final population is about 190.

At the optimum the village produces food for about 171 people but has housing for about 232. The binding constraint is food, but housing must be large enough that it does not choke growth. Substituting the equilibrium population $P^{*} \approx 191$ into the balance condition from Section 2 gives about $0.0616$ on both sides, so the simulation has settled at the equilibrium predicted by the formula.

The year-by-year trajectory rises quickly at first (100 → 120 → 140 → 160 → 177 → 188) and then flattens out around 191 from year 6 on. In other words, the village reaches its carrying capacity well before the ten years are up, and the optimization is really about choosing which carrying capacity it reaches.

6. Graph Output

Reading the graphs

(a) Final population vs. farmland. The curve rises steeply while food is scarce, peaks near 37 ha, and then falls sharply as housing shrinks. The drop on the right is much steeper than the climb on the left, which shows that over-farming is more dangerous than under-farming: housing vanishes linearly, and a population cannot live in fields it has no room to inhabit.

(b) Population trajectories. Each line is a different land split. Splits near the optimum (30–45 ha) all climb toward about 190. Small farmland shares flatten at low levels because of famine. The 50 ha line stays flat at 100 because the housing cap equals the starting population.

(c) Food limit vs. housing limit. The orange curve (people the farmland can feed) rises with diminishing returns. The blue line (people the housing can hold) falls linearly. The dotted black curve is the minimum of the two, which is the theoretical ceiling. The two curves cross at roughly 41 ha, but the red dashed optimum sits to the left of that crossing, at about 37 ha. The reason is that a population does not suddenly stop at its ceiling. Near the housing cap, logistic growth slows down, so giving up a little food to gain more room is worth it.

(d) 3D surface: time × farmland. The ridge of high population runs along the farmland axis near 35–40 ha and grows with time, then flattens by year 6 or so. The cliff on the high-farmland side marks housing exhaustion, and the low plateau on the low-farmland side marks famine. This is the whole dynamic of the problem in one picture.

(e) 3D surface: farmland × productivity. As productivity $A$ rises, the whole surface lifts and the peak moves toward less farmland. The land needed to feed a given number of people shrinks when each hectare is more productive, so the freed land can become housing.

(f) Heat map with optimal ridge. The white line is the optimal farmland as a function of productivity. It bends toward smaller farmland as $A$ grows, confirming the 3D picture. Better farming technology is not a reason to farm more land. It is a reason to build more houses.

7. Conclusion

With the 60 ha village in this example, the best strategy is to devote roughly 61% of the land to farming and 39% to housing, which brings the population from 100 to about 191 in ten years. The central lesson is that the answer is neither “maximize food” nor “maximize housing”. Growth is capped by the scarcer of the two resources, and the optimum balances them, tilted slightly toward food because starvation is more punishing than crowding. As farming becomes more productive, the best allocation shifts toward housing.

Designing Optimal Probes in Quantum Metrology

October 20, 2026

Maximizing Quantum Fisher Information with Python

How precisely can you measure a phase with $N$ qubits? The answer depends on two design choices: which input state you prepare and which measurement you perform. In this post we work through a concrete example, an $N$-qubit phase estimator, and solve the probe design problem numerically, both without noise and under local dephasing.

1. The Problem

A parameter $\varphi$ is imprinted on a probe state through the unitary

$$
U_\varphi = e^{-i\varphi G}, \qquad G=\frac12\sum_{i=1}^{N}\sigma_z^{(i)} .
$$

The quantum Cramér–Rao bound limits the variance of any unbiased estimator after $\nu$ repetitions:

$$
\mathrm{Var}(\hat\varphi);\ge;\frac{1}{\nu,F_Q[\rho_\varphi]} .
$$

The quantum Fisher information (QFI) $F_Q$ depends only on the probe state. For a pure state it is

$$
F_Q = 4\left(\langle G^2\rangle-\langle G\rangle^2\right).
$$

For a mixed state $\rho=\sum_i\lambda_i|i\rangle\langle i|$ it becomes

$$
F_Q = 2\sum_{i,j}\frac{(\lambda_i-\lambda_j)^2}{\lambda_i+\lambda_j},\big|\langle i|G|j\rangle\big|^2 .
$$

Step 1: Choosing the input state

The eigenvalues of $G$ range from $-N/2$ to $+N/2$, so the variance is bounded by $\mathrm{Var}(G)\le N^2/4$ and therefore

$$
F_Q\le N^2 \quad(\text{Heisenberg limit}).
$$

A product state such as $|+\rangle^{\otimes N}$ only reaches $F_Q=N$ (the standard quantum limit). The state that saturates the bound is the GHZ state

$$
|\mathrm{GHZ}\rangle=\frac{|0\rangle^{\otimes N}+|1\rangle^{\otimes N}}{\sqrt2}.
$$

Step 2: Choosing the measurement

The QFI is an upper bound. A concrete measurement with outcome probabilities $p_k(\varphi)$ yields the classical Fisher information

$$
F_C(\varphi)=\sum_k\frac{\big(\partial_\varphi p_k\big)^2}{p_k};\le;F_Q .
$$

We measure the same tilted Pauli observable on every qubit,

$$
\sigma_\vartheta=\sin\vartheta,\sigma_x+\cos\vartheta,\sigma_z ,
$$

where $\vartheta=0$ is the computational basis (blind to the phase of a GHZ state) and $\vartheta=\pi/2$ is the $x$-basis.

Step 3: Adding noise

Independent dephasing with probability $p$ on each qubit multiplies every density-matrix element by a factor that depends on the Hamming distance $d_H(i,j)$ between the two bit strings:

$$
\rho_{ij};\longrightarrow;(1-2p)^{,d_H(i,j)},\rho_{ij}.
$$

Because dephasing commutes with $U_\varphi$, we can compute the QFI of the noisy state directly. For the two reference probes the closed forms are

$$
F_Q^{\mathrm{GHZ}}=N^2(1-2p)^{2N},\qquad F_Q^{\mathrm{prod}}=N(1-2p)^2 .
$$

GHZ is better only while $(1-2p)^{2N-2}>1/N$. For $N=3$ the crossover is at $p=\tfrac12\left(1-3^{-1/4}\right)\approx 0.120$. Beyond that point the optimal probe must be found numerically.

2. Full Python Code

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

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

# ------------------------------------------------------------------
# 1. Basic building blocks
# ------------------------------------------------------------------
def popcount_table(N):
    return np.array([bin(k).count("1") for k in range(2**N)])

def generator_diag(N):
    # G = (1/2) * sum_i Z_i is diagonal; |0> -> +1/2, |1> -> -1/2
    return 0.5 * (N - 2 * popcount_table(N))

def dephasing_kernel(N, p):
    # Local dephasing: rho_ij -> rho_ij * (1-2p)^{Hamming distance(i, j)}
    pc = popcount_table(N)
    idx = np.arange(2**N)
    ham = pc[idx[:, None] ^ idx[None, :]]
    return (1.0 - 2.0 * p) ** ham

def pure_qfi(psi, g):
    prob = np.abs(psi) ** 2
    return 4.0 * (np.sum(prob * g**2) - np.sum(prob * g) ** 2)

def mixed_qfi(rho, g):
    lam, U = np.linalg.eigh(rho)
    Gm = U.conj().T @ (g[:, None] * U)
    S = lam[:, None] + lam[None, :]
    D = lam[:, None] - lam[None, :]
    W = np.zeros_like(S)
    mask = S > 1e-12
    W[mask] = D[mask] ** 2 / S[mask]
    return 2.0 * np.sum(W * np.abs(Gm) ** 2)

def ghz_state(N):
    psi = np.zeros(2**N, dtype=complex)
    psi[0] = psi[-1] = 1 / np.sqrt(2)
    return psi

def product_state(N):
    return np.ones(2**N, dtype=complex) / np.sqrt(2**N)

def to_state(x):
    d = len(x) // 2
    psi = x[:d] + 1j * x[d:]
    return psi / np.linalg.norm(psi)

# ------------------------------------------------------------------
# 2. Numerical search for the optimal probe state
# ------------------------------------------------------------------
def optimize_pure(N, n_restarts=4):
    g = generator_diag(N)
    best_val, best_psi = -1.0, None
    for _ in range(n_restarts):
        x0 = rng.normal(size=2 * 2**N)
        res = minimize(lambda x: -pure_qfi(to_state(x), g), x0, method="L-BFGS-B")
        if -res.fun > best_val:
            best_val, best_psi = -res.fun, to_state(res.x)
    return best_val, best_psi

def optimize_noisy(N, p, n_restarts=3, x_warm=None):
    g = generator_diag(N)
    K = dephasing_kernel(N, p)
    def objective(x):
        psi = to_state(x)
        rho = np.outer(psi, psi.conj()) * K
        return -mixed_qfi(rho, g)
    starts = [rng.normal(size=2 * 2**N) for _ in range(n_restarts)]
    if x_warm is not None:
        starts.append(x_warm)
    best_val, best_x = -1.0, None
    for x0 in starts:
        res = minimize(objective, x0, method="L-BFGS-B")
        if -res.fun > best_val:
            best_val, best_x = -res.fun, res.x
    return best_val, best_x

# ------------------------------------------------------------------
# 3. Classical Fisher information for a local projective measurement
# ------------------------------------------------------------------
def classical_fisher(psi, g, Mfull, phi):
    psi_phi = np.exp(-1j * phi * g) * psi
    dpsi = -1j * g * psi_phi
    a = Mfull @ psi_phi
    da = Mfull @ dpsi
    prob = np.abs(a) ** 2
    dprob = 2.0 * np.real(np.conj(a) * da)
    mask = prob > 1e-12
    return np.sum(dprob[mask] ** 2 / prob[mask])

def measurement_matrix(N, theta):
    c, s = np.cos(theta / 2), np.sin(theta / 2)
    Mt = np.array([[c, s], [-s, c]])
    M = Mt
    for _ in range(N - 1):
        M = np.kron(M, Mt)
    return M

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

# (A) Noiseless: optimal pure probe for N = 1..5
Ns = np.arange(1, 6)
qfi_opt_A = []
for N in Ns:
    val, _ = optimize_pure(int(N))
    qfi_opt_A.append(val)
qfi_opt_A = np.array(qfi_opt_A)

print("=== (A) Noiseless case ===")
print(f"{'N':>3} {'product':>9} {'GHZ':>9} {'optimized':>11} {'N^2':>6}")
for N, v in zip(Ns, qfi_opt_A):
    g = generator_diag(int(N))
    print(f"{N:>3} {pure_qfi(product_state(int(N)), g):>9.4f} "
          f"{pure_qfi(ghz_state(int(N)), g):>9.4f} {v:>11.4f} {N**2:>6d}")

# (B) Dephasing: N = 3, optimize the probe for each p (warm start)
N3 = 3
g3 = generator_diag(N3)
p_grid = np.linspace(0.0, 0.45, 10)
qfi_opt_B, qfi_ghz_B, qfi_prod_B, x_store = [], [], [], []
x_warm = None
for p in p_grid:
    K = dephasing_kernel(N3, p)
    val, x_warm = optimize_noisy(N3, p, n_restarts=3, x_warm=x_warm)
    qfi_opt_B.append(val)
    x_store.append(x_warm.copy())
    qfi_ghz_B.append(mixed_qfi(np.outer(ghz_state(N3), ghz_state(N3).conj()) * K, g3))
    qfi_prod_B.append(mixed_qfi(np.outer(product_state(N3), product_state(N3).conj()) * K, g3))
qfi_opt_B, qfi_ghz_B, qfi_prod_B = map(np.array, (qfi_opt_B, qfi_ghz_B, qfi_prod_B))

print("\n=== (B) Local dephasing, N = 3 ===")
print(f"{'p':>6} {'product':>9} {'GHZ':>9} {'optimized':>11}")
for p, a, b, c in zip(p_grid, qfi_prod_B, qfi_ghz_B, qfi_opt_B):
    print(f"{p:>6.3f} {a:>9.4f} {b:>9.4f} {c:>11.4f}")

# (C) Measurement: tilted local measurement for GHZ and product probes
theta_grid = np.linspace(0.05, np.pi - 0.05, 60)
phi_grid = np.linspace(0.011, np.pi / N3 - 0.011, 80)
psi_ghz, psi_prod = ghz_state(N3), product_state(N3)
FC_ghz = np.zeros((len(theta_grid), len(phi_grid)))
FC_prod = np.zeros_like(FC_ghz)
for i, th in enumerate(theta_grid):
    M = measurement_matrix(N3, th)
    for j, ph in enumerate(phi_grid):
        FC_ghz[i, j] = classical_fisher(psi_ghz, g3, M, ph)
        FC_prod[i, j] = classical_fisher(psi_prod, g3, M, ph)
ratio_ghz = FC_ghz.mean(axis=1) / N3**2
ratio_prod = FC_prod.mean(axis=1) / N3

i_mid = len(theta_grid) // 2
print("\n=== (C) Measurement ===")
print(f"GHZ, computational-basis-like tilt (theta ~ 0): mean F_C = {FC_ghz[0].mean():.4f}")
print(f"GHZ, near-equatorial tilt (theta ~ pi/2):       mean F_C = {FC_ghz[i_mid].mean():.4f}  (F_Q = {N3**2})")
print(f"Product, near-equatorial tilt (theta ~ pi/2):   mean F_C = {FC_prod[i_mid].mean():.4f}  (F_Q = {N3})")

# (D) Entanglement advantage over (N, p)
N_list = np.arange(1, 7)
p_list = np.linspace(0.0, 0.45, 31)
ratio_surface = np.zeros((len(N_list), len(p_list)))
max_dev = 0.0
for a, N in enumerate(N_list):
    N = int(N)
    g = generator_diag(N)
    pg, pp = ghz_state(N), product_state(N)
    for b, p in enumerate(p_list):
        K = dephasing_kernel(N, p)
        fg = mixed_qfi(np.outer(pg, pg.conj()) * K, g)
        fp = mixed_qfi(np.outer(pp, pp.conj()) * K, g)
        ratio_surface[a, b] = fg / fp
        max_dev = max(max_dev, abs(fg - N**2 * (1 - 2*p)**(2*N)),
                      abs(fp - N * (1 - 2*p)**2))
p_star = 0.5 * (1 - 3 ** (-1 / 4))
print("\n=== (D) Closed-form check ===")
print(f"max |numerical - analytic| over the (N, p) grid = {max_dev:.2e}")
print(f"N = 3: GHZ stops beating the product state near p = {p_star:.4f}")
print(f"Total computation time: {time.time() - t0:.1f} s")

# ------------------------------------------------------------------
# 5. One figure with all results
# ------------------------------------------------------------------
fig = plt.figure(figsize=(20, 12))
gs = GridSpec(2, 3, figure=fig, hspace=0.32, wspace=0.28)

ax1 = fig.add_subplot(gs[0, 0])
Nf = np.linspace(1, 5, 100)
ax1.plot(Nf, Nf, "b--", label=r"Standard quantum limit $N$")
ax1.plot(Nf, Nf**2, "r-", label=r"Heisenberg limit $N^2$")
ax1.plot(Ns, qfi_opt_A, "ko", ms=9, mfc="none", mew=2, label="Numerically optimized probe")
ax1.set_xlabel("Number of qubits $N$")
ax1.set_ylabel(r"Quantum Fisher information $F_Q$")
ax1.set_title("(a) Noiseless optimal probe")
ax1.legend()
ax1.grid(alpha=0.3)

ax2 = fig.add_subplot(gs[0, 1])
ax2.plot(p_grid, qfi_prod_B, "b--o", label=r"Product state $|+\rangle^{\otimes 3}$")
ax2.plot(p_grid, qfi_ghz_B, "r-s", label="GHZ state")
ax2.plot(p_grid, qfi_opt_B, "k-", lw=2.5, label="Optimized probe")
ax2.set_xlabel("Dephasing probability $p$")
ax2.set_ylabel(r"$F_Q$")
ax2.set_title("(b) Probe optimization under dephasing ($N=3$)")
ax2.legend()
ax2.grid(alpha=0.3)

ax3 = fig.add_subplot(gs[0, 2])
labels = [format(k, "03b") for k in range(2**N3)]
width = 0.27
for off, ip in zip((-width, 0.0, width), (0, 3, 6)):
    probs = np.abs(to_state(x_store[ip])) ** 2
    ax3.bar(np.arange(2**N3) + off, probs, width=width, label=f"$p={p_grid[ip]:.2f}$")
ax3.set_xticks(np.arange(2**N3))
ax3.set_xticklabels([f"|{s}>" for s in labels])
ax3.set_ylabel("Population of the optimal probe")
ax3.set_title("(c) Structure of the optimal input state")
ax3.legend()
ax3.grid(alpha=0.3, axis="y")

ax4 = fig.add_subplot(gs[1, 0], projection="3d")
PH, TH = np.meshgrid(phi_grid, theta_grid)
surf = ax4.plot_surface(PH, TH, FC_ghz, cmap="viridis", edgecolor="none", alpha=0.95)
ax4.set_xlabel(r"Phase $\phi$")
ax4.set_ylabel(r"Measurement tilt $\vartheta$")
ax4.set_zlabel(r"$F_C$", labelpad=6)
ax4.set_title("(d) GHZ probe: classical Fisher information (3D)")
ax4.view_init(elev=28, azim=-60)
fig.colorbar(surf, ax=ax4, shrink=0.55, pad=0.08)

ax5 = fig.add_subplot(gs[1, 1])
ax5.plot(theta_grid, ratio_ghz, "r-", lw=2.5, label=r"GHZ: $\bar F_C / N^2$")
ax5.plot(theta_grid, ratio_prod, "b--", lw=2.5, label=r"Product: $\bar F_C / N$")
ax5.axvline(np.pi / 2, color="gray", ls=":")
ax5.set_xlabel(r"Measurement tilt $\vartheta$ (rad)")
ax5.set_ylabel(r"$\bar F_C / F_Q$ (averaged over $\phi$)")
ax5.set_title("(e) How much of $F_Q$ does the measurement extract?")
ax5.legend()
ax5.grid(alpha=0.3)

ax6 = fig.add_subplot(gs[1, 2], projection="3d")
PP, NN = np.meshgrid(p_list, N_list)
surf2 = ax6.plot_surface(PP, NN, ratio_surface, cmap="plasma", edgecolor="none", alpha=0.95)
ax6.plot_surface(PP, NN, np.ones_like(ratio_surface), color="gray", alpha=0.25)
ax6.set_xlabel("Dephasing $p$")
ax6.set_ylabel("Qubits $N$")
ax6.set_zlabel("GHZ / product", labelpad=6)
ax6.set_title("(f) Entanglement advantage (gray plane = 1)")
ax6.view_init(elev=26, azim=-125)
fig.colorbar(surf2, ax=ax6, shrink=0.55, pad=0.08)

fig.suptitle("Optimal probe design in quantum metrology", fontsize=18, y=0.98)
plt.show()

3. Code Walkthrough

Building blocks

Since $G=\tfrac12\sum_i\sigma_z^{(i)}$ is diagonal in the computational basis, it is stored as a vector g. The basis state with $n_1$ ones has eigenvalue $(N-2n_1)/2$, and popcount_table counts those ones.

dephasing_kernel is the main speed-up. Instead of applying $3N$ Kraus operators one at a time, it builds the matrix $K_{ij}=(1-2p)^{d_H(i,j)}$ once. The Hamming distance comes from the bitwise XOR of the indices, idx[:,None] ^ idx[None,:]. The noisy state is then a single element-wise product, rho * K.

pure_qfi implements $4,\mathrm{Var}(G)$ using only the populations $|\psi_k|^2$. mixed_qfi implements the eigendecomposition formula. The mask S > 1e-12 removes the $0/0$ terms that appear when two eigenvalues vanish, which would otherwise produce NaN.

Parametrizing and optimizing the probe

to_state maps a real vector of length $2\cdot 2^N$ to a normalized complex state. This turns the constrained search over unit vectors into an unconstrained problem that L-BFGS-B can handle.

optimize_pure minimizes $-F_Q$ from several random starting points and keeps the best result. It uses the cheap pure-state formula, which is why $N=5$ (a 32-dimensional state) finishes almost instantly.

optimize_noisy does the same for the noisy QFI, but with two additions that speed things up. It evaluates $F_Q$ on rho * K, and it accepts a warm start x_warm. As $p$ increases in small steps, the previous optimum is an excellent initial guess, so only three random restarts per $p$ are needed.

Measurement model

measurement_matrix builds $M^{\dagger\otimes N}$, where $M=R_y(\vartheta)$ has the eigenvectors of $\sigma_\vartheta$ as columns. Measuring $\sigma_\vartheta$ on every qubit means the outcome probabilities are $p_k=|(M^{\dagger\otimes N}\psi_\varphi)_k|^2$.

classical_fisher computes the derivative analytically,

$$
\partial_\varphi p_k = 2,\mathrm{Re}!\left[a_k^*,\partial_\varphi a_k\right],\qquad \partial_\varphi\psi_\varphi=-iG\psi_\varphi ,
$$

so no finite-difference error enters. Outcomes with $p_k\le10^{-12}$ are masked out to avoid division by zero.

Experiments (A)–(D)

(A) For $N=1,\dots,5$ the code optimizes over all pure states and compares the result with $N$ and $N^2$.
(B) For $N=3$ and ten values of $p$ it compares the product state, GHZ and the numerically optimized probe.
(C) It sweeps the measurement tilt $\vartheta$ and the phase $\varphi$ on a $60\times80$ grid. The phase range $0<\varphi<\pi/3$ covers one period of the GHZ signal $\cos 3\varphi$. The grids are offset slightly so they avoid the exact points where outcome probabilities vanish.
(D) It evaluates the ratio $F_Q^{\rm GHZ}/F_Q^{\rm prod}$ for $N=1,\dots,6$ and 31 values of $p$. It also reports the largest deviation from the closed-form expressions, which serves as a built-in correctness check of mixed_qfi.

The figure

All six panels are drawn in a single GridSpec figure. Panels (d) and (f) are 3D surfaces.

4. Results

=== (A) Noiseless case ===
  N   product       GHZ   optimized    N^2
  1    1.0000    1.0000      1.0000      1
  2    2.0000    4.0000      4.0000      4
  3    3.0000    9.0000      9.0000      9
  4    4.0000   16.0000     16.0000     16
  5    5.0000   25.0000     25.0000     25

=== (B) Local dephasing, N = 3 ===
     p   product       GHZ   optimized
 0.000    3.0000    9.0000      9.0000
 0.050    2.4300    4.7830      4.7830
 0.100    1.9200    2.3593      2.6425
 0.150    1.4700    1.0588      1.7240
 0.200    1.0800    0.4199      1.1678
 0.250    0.7500    0.1406      0.7770
 0.300    0.4800    0.0369      0.4867
 0.350    0.2700    0.0066      0.2711
 0.400    0.1200    0.0006      0.1201
 0.450    0.0300    0.0000      0.0300

=== (C) Measurement ===
GHZ, computational-basis-like tilt (theta ~ 0): mean F_C = 0.0000
GHZ, near-equatorial tilt (theta ~ pi/2):       mean F_C = 8.7309  (F_Q = 9)
Product, near-equatorial tilt (theta ~ pi/2):   mean F_C = 2.8973  (F_Q = 3)

=== (D) Closed-form check ===
max |numerical - analytic| over the (N, p) grid = 2.13e-14
N = 3: GHZ stops beating the product state near p = 0.1201
Total computation time: 2.1 s

5. Discussion of the Results

(a) Noiseless optimum: the Heisenberg limit emerges from the optimizer

The black circles are the best QFI values found by the optimizer for $N=1,\dots,5$. They sit on the red curve $N^2$ and well above the blue line $N$. The optimizer started from random states and was never told about GHZ, yet it reproduced the bound $F_Q\le(g_{\max}-g_{\min})^2=N^2$. For $N=3$ the values are $3$ for the product state and $9$ for both GHZ and the optimized probe.

(b) Under dephasing, GHZ is not always best

At $p=0$ the optimized probe equals GHZ ($F_Q=9$). GHZ coherence decays as $(1-2p)^{N}$, so its QFI collapses like $(1-2p)^{2N}$ while the product state’s QFI decays only like $(1-2p)^2$. The curves cross near $p\approx0.12$. In the intermediate regime ($p\approx0.1$–$0.3$) the optimized probe is strictly better than both reference states. At $p=0.10$, for example, the optimized value is about $2.64$, against roughly $2.36$ for GHZ and $1.92$ for the product state. For very large $p$ the optimum merges with the product state. No fixed textbook state is optimal across the whole range, which is the point of probe optimization.

(c) Structure of the optimal state

At $p=0$ the optimizer produces a GHZ state, with weight $0.5$ on $|000\rangle$ and $|111\rangle$ and nothing elsewhere. As noise grows, population leaks into the other six basis states. The weight on the two extreme states drops while the rest rises toward a uniform distribution, which is the signature of a product state $|+\rangle^{\otimes3}$ (every population equal to $1/8$). The optimal probe therefore interpolates between GHZ and a product state, shifting toward the more robust end as noise increases.

(d) 3D view: the measurement matters as much as the state

This surface shows the classical Fisher information $F_C(\varphi,\vartheta)$ of the GHZ probe. Near $\vartheta=\pi/2$, the $x$-basis, there is a high plateau close to $F_Q=9$. As $\vartheta\to0$ or $\pi$ the surface falls to zero. In the computational basis the outcome statistics of a GHZ state do not depend on $\varphi$ at all, so a perfect probe combined with the wrong measurement carries no information. The dips at the edges of the plateau, near $\varphi\to0$ and $\varphi\to\pi/3$, appear because some outcome probabilities approach zero there.

(e) How much of the QFI the measurement extracts

Averaged over $\varphi$, both the GHZ and product curves peak near $\vartheta=\pi/2$ at roughly $0.9$–$1.0$ of their respective $F_Q$. The GHZ curve is narrower, so entangled probes demand more precise measurement alignment. Entanglement raises the achievable sensitivity but makes it more fragile with respect to measurement errors.

(f) 3D view: where entanglement pays off

The surface is the ratio $F_Q^{\rm GHZ}/F_Q^{\rm prod}=N(1-2p)^{2N-2}$, and the gray plane marks ratio $=1$. In the noiseless limit the ratio equals $N$, so it grows linearly with the number of qubits. For any $p>0$ the exponent $2N-2$ makes larger GHZ states lose their advantage faster. The surface drops below the plane in a region that expands as $N$ increases, so the best probe size depends on the noise level. The console check confirms that the numerical QFI agrees with the closed forms to machine precision.

6. Takeaways

The QFI bounds the best possible precision, so maximizing it over input states is the first step of probe design. In the noiseless case this gives the Heisenberg-limited GHZ state.
The measurement must be matched to the probe. A GHZ state read out in the computational basis yields $F_C=0$, while an $x$-basis readout approaches $F_Q$.
With noise, the optimal state is neither GHZ nor a product state. It is a numerically found compromise, and a short gradient-based optimization over $2^N$ amplitudes is enough to locate it.
The same pipeline carries over to other generators, noise channels and qubit numbers. Only generator_diag and dephasing_kernel need to change.

Variational Quantum Algorithms in Practice

October 19, 2026

Solving a Hydrogen-Molecule VQE and a MaxCut QAOA from Scratch in Python

Variational quantum algorithms (VQE and QAOA) split the work between two machines. A parameterized quantum circuit prepares a state $|\psi(\boldsymbol{\theta})\rangle$ and measures a cost function, and a classical optimizer updates $\boldsymbol{\theta}$ to minimize it. This hybrid loop is the standard strategy for today’s noisy, small-scale quantum devices.

In this article we solve two concrete problems with nothing but NumPy and SciPy:

VQE finds the ground-state energy of a two-qubit hydrogen-molecule Hamiltonian.
QAOA solves a MaxCut problem on a 6-node graph, with circuit depths $p = 1, 2, 3, 4$.

Everything is simulated with an exact state-vector, so no quantum SDK is required.

1. Theory

1.1 VQE (Variational Quantum Eigensolver)

Given a Hamiltonian $H$, VQE minimizes the expectation value

$$
E(\boldsymbol{\theta}) = \langle \psi(\boldsymbol{\theta}) | H | \psi(\boldsymbol{\theta}) \rangle ;\ge; E_0 ,
$$

where $E_0$ is the true ground-state energy. The inequality is the variational principle: any trial state gives an upper bound on $E_0$, so minimizing over $\boldsymbol{\theta}$ pushes the energy down toward the ground state.

The Hamiltonian is written as a sum of Pauli strings,

$$
H = \sum_k c_k P_k, \qquad P_k \in {I, X, Y, Z}^{\otimes n},
$$

and on hardware each $\langle P_k\rangle$ is estimated separately from measurements. Our simulator computes the full expectation exactly.

1.2 QAOA (Quantum Approximate Optimization Algorithm)

For MaxCut on a graph $G=(V,E)$, the cost operator is

$$
C = \sum_{(u,v)\in E} \frac{1 - Z_u Z_v}{2},
$$

whose eigenvalue on a computational basis state $|z\rangle$ equals the number of cut edges. A depth-$p$ QAOA state is

$$
|\boldsymbol{\gamma},\boldsymbol{\beta}\rangle = \prod_{l=1}^{p} e^{-i\beta_l B}, e^{-i\gamma_l C}, |+\rangle^{\otimes n}, \qquad B = \sum_{j} X_j .
$$

The classical optimizer maximizes

$$
F_p(\boldsymbol{\gamma},\boldsymbol{\beta}) = \langle \boldsymbol{\gamma},\boldsymbol{\beta} | C | \boldsymbol{\gamma},\boldsymbol{\beta} \rangle ,
$$

and the quality of the result is measured by the approximation ratio

$$
r_p = \frac{F_p^{*}}{C_{\max}} \in [0,1].
$$

As $p\to\infty$, QAOA is guaranteed to approach $r=1$, so a deeper circuit should never do worse.

2. Problem Setup

Problem A: the H$_2$ molecule (VQE)

At a bond length of $0.735,\text{Å}$, after parity mapping and two-qubit reduction, the electronic Hamiltonian of H$_2$ becomes

$$
H = g_0, II + g_1, IZ + g_2, ZI + g_3, ZZ + g_4, XX,
$$

with

$$
g_0=-1.0524,\quad g_1=0.3979,\quad g_2=-0.3979,\quad g_3=-0.0113,\quad g_4=0.1809 \ \text{(Hartree)}.
$$

Our ansatz uses four parameters:

$$
|\psi(\boldsymbol{\theta})\rangle = \big(R_y(\theta_2)\otimes R_y(\theta_3)\big);\mathrm{CNOT}_{0\to1};\big(R_y(\theta_0)\otimes R_y(\theta_1)\big),|00\rangle ,
$$

where

$$
R_y(\theta)=\begin{pmatrix}\cos\frac{\theta}{2} & -\sin\frac{\theta}{2}\[2pt] \sin\frac{\theta}{2} & \cos\frac{\theta}{2}\end{pmatrix}.
$$

Because $H$ is a real matrix, its ground state is real. Every real two-qubit state can be written as a local real rotation applied to $\cos a,|00\rangle+\sin a,|11\rangle$, so this ansatz is expressive enough to reach the exact ground state.

Problem B: MaxCut on a 6-node graph (QAOA)

The graph has 6 nodes and 8 edges: a ring $0!-!1!-!2!-!3!-!4!-!5!-!0$ plus the two chords $(0,3)$ and $(1,4)$. We find the best parameters for $p=1,2,3,4$ and compare the result with the brute-force optimum.

3. Full Source Code

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

plt.style.use('dark_background')
rng = np.random.default_rng(7)
t_start = time.time()

# =====================================================================
# Part A : VQE for the two-qubit H2 Hamiltonian
# =====================================================================
I2 = np.eye(2, dtype=complex)
X = np.array([[0, 1], [1, 0]], dtype=complex)
Z = np.array([[1, 0], [0, -1]], dtype=complex)


def kron(*mats):
    out = np.array([[1.0 + 0j]])
    for m in mats:
        out = np.kron(out, m)
    return out


def ry(t):
    c, s = np.cos(t / 2.0), np.sin(t / 2.0)
    return np.array([[c, -s], [s, c]], dtype=complex)


CNOT = np.array([[1, 0, 0, 0],
                 [0, 1, 0, 0],
                 [0, 0, 0, 1],
                 [0, 0, 1, 0]], dtype=complex)

g0, g1, g2, g3, g4 = (-1.052373245772859, 0.39793742484318045,
                      -0.39793742484318045, -0.01128010425623538,
                      0.18093119978423156)

H = (g0 * kron(I2, I2) + g1 * kron(I2, Z) + g2 * kron(Z, I2)
     + g3 * kron(Z, Z) + g4 * kron(X, X))

E_exact = float(np.linalg.eigvalsh(H)[0])
psi_init = np.array([1, 0, 0, 0], dtype=complex)


def vqe_state(th):
    u1 = kron(ry(th[0]), ry(th[1]))
    u2 = kron(ry(th[2]), ry(th[3]))
    return u2 @ (CNOT @ (u1 @ psi_init))


def vqe_energy(th):
    psi = vqe_state(th)
    return float(np.real(np.vdot(psi, H @ psi)))


histories, results = [], []
for run in range(5):
    hist = []

    def objective(th, hist=hist):
        e = vqe_energy(th)
        hist.append(e)
        return e

    x0 = rng.uniform(-np.pi, np.pi, 4)
    res = minimize(objective, x0, method='COBYLA', tol=1e-10,
                   options={'maxiter': 400, 'rhobeg': 0.5})
    histories.append(np.array(hist))
    results.append(res)

best_run = int(np.argmin([r.fun for r in results]))
th_opt = results[best_run].x
E_vqe = results[best_run].fun

print("=" * 60)
print("Part A : VQE (H2, bond length 0.735 A)")
print("=" * 60)
print(f"Exact ground energy (diagonalization) : {E_exact:.10f} Ha")
print(f"VQE energy (best of 5 restarts)       : {E_vqe:.10f} Ha")
print(f"Absolute error                        : {abs(E_vqe - E_exact):.3e} Ha")
print("Optimal parameters (rad)              :", np.round(th_opt, 5))

# 3D landscape: vary theta_0, theta_1 with theta_2, theta_3 fixed at optimum
grid_n = 61
t0_grid = np.linspace(-np.pi, np.pi, grid_n)
t1_grid = np.linspace(-np.pi, np.pi, grid_n)
T0, T1 = np.meshgrid(t0_grid, t1_grid)
E_land = np.zeros_like(T0)
for i in range(grid_n):
    for j in range(grid_n):
        E_land[i, j] = vqe_energy(np.array([T0[i, j], T1[i, j],
                                            th_opt[2], th_opt[3]]))

w0 = float(np.arctan2(np.sin(th_opt[0]), np.cos(th_opt[0])))
w1 = float(np.arctan2(np.sin(th_opt[1]), np.cos(th_opt[1])))
E_mark = vqe_energy(np.array([w0, w1, th_opt[2], th_opt[3]]))

# =====================================================================
# Part B : QAOA for MaxCut on a 6-node graph
# =====================================================================
n = 6
dim = 2 ** n
edges = [(0, 1), (1, 2), (2, 3), (3, 4), (4, 5), (5, 0), (0, 3), (1, 4)]

idx_all = np.arange(dim)
cut = np.zeros(dim)
for (u, v) in edges:
    cut += (((idx_all >> u) & 1) != ((idx_all >> v) & 1)).astype(float)
max_cut = float(cut.max())
n_opt = int(np.sum(cut == max_cut))


def apply_rx_all(psi, beta):
    c, s = np.cos(beta), np.sin(beta)
    U = np.array([[c, -1j * s], [-1j * s, c]], dtype=complex)
    t = psi.reshape([2] * n)
    for ax in range(n):
        t = np.moveaxis(np.tensordot(U, t, axes=([1], [ax])), 0, ax)
    return t.reshape(-1)


def qaoa_state(params, p):
    gammas, betas = params[:p], params[p:]
    psi = np.full(dim, 1.0 / np.sqrt(dim), dtype=complex)
    for g, b in zip(gammas, betas):
        psi = np.exp(-1j * g * cut) * psi
        psi = apply_rx_all(psi, b)
    return psi


def qaoa_expect(params, p):
    psi = qaoa_state(params, p)
    return float(np.real(np.sum(np.abs(psi) ** 2 * cut)))


def neg_expect(params, p):
    return -qaoa_expect(params, p)


P_MAX = 4
N_RESTART = 10
best_params, best_val = {}, {}
for p in range(1, P_MAX + 1):
    best = None
    for r in range(N_RESTART):
        x0 = np.concatenate([rng.uniform(0, np.pi, p),
                             rng.uniform(0, np.pi / 2, p)])
        res = minimize(neg_expect, x0, args=(p,), method='L-BFGS-B')
        if best is None or res.fun < best.fun:
            best = res
    best_params[p] = best.x
    best_val[p] = -best.fun

ps = np.arange(1, P_MAX + 1)
ratios = np.array([best_val[p] / max_cut for p in ps])

print()
print("=" * 60)
print("Part B : QAOA (MaxCut, 6 nodes, 8 edges)")
print("=" * 60)
print(f"Brute-force maximum cut : {max_cut:.0f}  ({n_opt} optimal bitstrings)")
print(f"{'p':>3} | {'<C>':>10} | {'ratio':>8}")
print("-" * 30)
for p in ps:
    print(f"{p:>3} | {best_val[p]:>10.5f} | {best_val[p] / max_cut:>8.5f}")

psi_final = qaoa_state(best_params[P_MAX], P_MAX)
probs = np.abs(psi_final) ** 2
order = np.argsort(probs)[::-1][:4]
print()
print(f"Top bitstrings at p = {P_MAX} (qubit 0 is the rightmost bit):")
for k in order:
    print(f"  |{k:0{n}b}>  prob = {probs[k]:.4f}   cut = {cut[k]:.0f}")
print(f"Probability of sampling an optimal cut : {probs[cut == max_cut].sum():.4f}")

# 3D landscape for p = 1
gam_grid = np.linspace(0, 2 * np.pi, 61)
bet_grid = np.linspace(0, np.pi, 61)
GG, BB = np.meshgrid(gam_grid, bet_grid)
F_land = np.zeros_like(GG)
for i in range(GG.shape[0]):
    for j in range(GG.shape[1]):
        F_land[i, j] = qaoa_expect(np.array([GG[i, j], BB[i, j]]), 1)

g1_opt = float(best_params[1][0] % (2 * np.pi))
b1_opt = float(best_params[1][1] % np.pi)
F_mark = qaoa_expect(np.array([g1_opt, b1_opt]), 1)

# =====================================================================
# Combined figure
# =====================================================================
fig = plt.figure(figsize=(21, 12))
fig.suptitle("Variational Quantum Algorithms: VQE (H$_2$) and QAOA (MaxCut)",
             fontsize=20, y=0.995)

# (1) VQE convergence
ax1 = fig.add_subplot(2, 3, 1)
for k, h in enumerate(histories):
    err = np.maximum(np.minimum.accumulate(h) - E_exact, 1e-12)
    if k == best_run:
        ax1.semilogy(err, color='cyan', lw=2.5, label='best run')
    else:
        ax1.semilogy(err, color='gray', lw=1.0, alpha=0.7)
ax1.axhline(1.6e-3, color='orange', ls='--', lw=1.5,
            label='chemical accuracy (1.6 mHa)')
ax1.set_xlabel("Cost-function evaluations")
ax1.set_ylabel(r"$E(\theta) - E_{\mathrm{exact}}$  [Ha]")
ax1.set_title("(1) VQE convergence (5 restarts)")
ax1.legend(loc='upper right')
ax1.grid(alpha=0.25)

# (2) VQE 3D landscape
ax2 = fig.add_subplot(2, 3, 2, projection='3d')
ax2.plot_surface(T0, T1, E_land, cmap='viridis', alpha=0.92,
                 linewidth=0, antialiased=True)
ax2.scatter([w0], [w1], [E_mark], color='red', s=90, depthshade=False,
            label='VQE optimum')
ax2.set_xlabel(r"$\theta_0$")
ax2.set_ylabel(r"$\theta_1$")
ax2.set_zlabel("Energy [Ha]")
ax2.set_title(r"(2) VQE energy landscape ($\theta_2,\theta_3$ fixed)")
ax2.legend(loc='upper left')
ax2.view_init(elev=28, azim=-60)

# (3) Graph with best cut
ax3 = fig.add_subplot(2, 3, 3)
best_idx = int(np.argmax(probs))
bits = [(best_idx >> k) & 1 for k in range(n)]
ang = 2 * np.pi * np.arange(n) / n
pos = np.column_stack([np.cos(ang), np.sin(ang)])
for (u, v) in edges:
    is_cut = bits[u] != bits[v]
    ax3.plot([pos[u, 0], pos[v, 0]], [pos[u, 1], pos[v, 1]],
             color='lime' if is_cut else 'gray',
             ls='-' if is_cut else '--', lw=3 if is_cut else 1.5, zorder=1)
for k in range(n):
    ax3.scatter(pos[k, 0], pos[k, 1], s=900,
                color='deepskyblue' if bits[k] == 0 else 'orange',
                edgecolor='white', zorder=2)
    ax3.text(pos[k, 0], pos[k, 1], str(k), ha='center', va='center',
             color='black', fontsize=14, fontweight='bold', zorder=3)
ax3.set_aspect('equal')
ax3.axis('off')
ax3.set_title(f"(3) Most probable cut at p={P_MAX}: {cut[best_idx]:.0f} of "
              f"{max_cut:.0f} edges")

# (4) QAOA p=1 3D landscape
ax4 = fig.add_subplot(2, 3, 4, projection='3d')
ax4.plot_surface(GG, BB, F_land, cmap='magma', alpha=0.92,
                 linewidth=0, antialiased=True)
ax4.scatter([g1_opt], [b1_opt], [F_mark], color='cyan', s=90,
            depthshade=False, label='p=1 optimum')
ax4.set_xlabel(r"$\gamma$")
ax4.set_ylabel(r"$\beta$")
ax4.set_zlabel(r"$\langle C \rangle$")
ax4.set_title(r"(4) QAOA landscape at $p=1$")
ax4.legend(loc='upper left')
ax4.view_init(elev=30, azim=-55)

# (5) Approximation ratio vs depth
ax5 = fig.add_subplot(2, 3, 5)
ax5.plot(ps, ratios, 'o-', color='cyan', lw=2.5, ms=10)
ax5.axhline(1.0, color='orange', ls='--', lw=1.5, label='optimum')
for p, r in zip(ps, ratios):
    ax5.text(p, r - 0.012, f"{r:.3f}", ha='center', va='top', fontsize=11)
ax5.set_xticks(ps)
ax5.set_ylim(min(ratios) - 0.06, 1.03)
ax5.set_xlabel("QAOA depth $p$")
ax5.set_ylabel("Approximation ratio")
ax5.set_title("(5) Approximation ratio vs depth")
ax5.legend(loc='lower right')
ax5.grid(alpha=0.25)

# (6) Output distribution at p = P_MAX
ax6 = fig.add_subplot(2, 3, 6)
colors = plt.cm.plasma(cut / max_cut)
ax6.bar(np.arange(dim), probs, color=colors, width=1.0)
sm = plt.cm.ScalarMappable(cmap='plasma', norm=plt.Normalize(0, max_cut))
sm.set_array([])
fig.colorbar(sm, ax=ax6, label="Cut value")
ax6.set_xlabel("Bitstring index $z$")
ax6.set_ylabel("Probability")
ax6.set_title(f"(6) Output distribution at $p={P_MAX}$")
ax6.grid(alpha=0.25)

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

print()
print(f"Total elapsed time: {time.time() - t_start:.1f} s")

4. Code Walkthrough

4.1 Part A: VQE

Hamiltonian construction. The helper kron chains Kronecker products, so kron(I2, Z) is $I\otimes Z$. The first factor is qubit 0 (the most significant bit). The five Pauli terms are assembled into a $4\times4$ matrix H. We also diagonalize H exactly with eigvalsh, giving E_exact as a reference to judge the VQE result.

Ansatz. vqe_state applies the four-parameter circuit directly to the vector $|00\rangle$: a layer of $R_y$ gates, a CNOT, and a second layer of $R_y$ gates. vqe_energy evaluates $\langle\psi|H|\psi\rangle$ with np.vdot, which conjugates its first argument. The result is real because $H$ is Hermitian.

Optimization loop. We run COBYLA from five random starting points, because variational landscapes can have local minima. COBYLA is gradient-free, which mirrors how shot-noise-limited hardware experiments are often run. Each run stores its full history of energy evaluations, so we can draw convergence curves later. The run with the lowest final energy is kept.

Landscape scan. To visualize the cost function, we vary $\theta_0$ and $\theta_1$ over $[-\pi,\pi]$ while holding $\theta_2,\theta_3$ at their optimal values. The energy has period $2\pi$ in every parameter, so the optimum is wrapped back into the plotted range with arctan2(sin, cos) before it is marked on the surface.

4.2 Part B: QAOA

Cost vector. The cost operator $C$ is diagonal in the computational basis. We therefore precompute cut, a vector of length $2^6=64$ holding the number of cut edges for every bitstring, using bit masks. The phase-separation step $e^{-i\gamma C}$ then becomes one element-wise multiplication, np.exp(-1j*g*cut) * psi.

Mixer. The mixer $e^{-i\beta B}=\bigotimes_j e^{-i\beta X_j}$ acts on each qubit independently. We reshape the 64-element vector into a tensor of shape $(2,2,2,2,2,2)$ and contract a $2\times2$ matrix into one axis at a time with tensordot. Building the full $64\times64$ mixer matrix is never necessary.

Optimization. For each depth $p$ we maximize $F_p$ by minimizing its negative with L-BFGS-B from 10 random starting points. Random restarts matter here because the QAOA landscape has many equivalent and inequivalent local optima.

Post-processing. The code evaluates the optimized state at $p=4$ and extracts the probability distribution $|\psi_z|^2$. It also computes the total probability of sampling a maximum cut, which is the quantity that matters when you actually measure the circuit.

4.3 Why this code is fast

The slow way to write this simulation is to build every gate as a full $2^n\times2^n$ matrix and multiply, which costs $O(4^n)$ per gate. Three choices avoid that:

the diagonal cost is stored as a vector, so the phase step costs $O(2^n)$;
the mixer is applied qubit by qubit through tensordot, costing $O(n,2^n)$ instead of $O(4^n)$;
the Pauli matrices for VQE are tiny ($4\times4$), so exact evaluation is essentially free.

With these choices both landscape scans (3,721 evaluations each) and all optimizations finish quickly on a standard Colab CPU.

5. Execution Results

============================================================
Part A : VQE (H2, bond length 0.735 A)
============================================================
Exact ground energy (diagonalization) : -1.8572750302 Ha
VQE energy (best of 5 restarts)       : -1.8572750302 Ha
Absolute error                        : 1.332e-15 Ha
Optimal parameters (rad)              : [-2.30741  1.87472 -3.95238  1.77846]

============================================================
Part B : QAOA (MaxCut, 6 nodes, 8 edges)
============================================================
Brute-force maximum cut : 8  (2 optimal bitstrings)
  p |         |    ratio
------------------------------
  1 |    5.63245 |  0.70406
  2 |    7.02977 |  0.87872
  3 |    7.04608 |  0.88076
  4 |    7.97840 |  0.99730

Top bitstrings at p = 4 (qubit 0 is the rightmost bit):
  |101010>  prob = 0.4968   cut = 8
  |010101>  prob = 0.4968   cut = 8
  |010001>  prob = 0.0007   cut = 6
  |101110>  prob = 0.0007   cut = 6
Probability of sampling an optimal cut : 0.9935


Total elapsed time: 17.3 s

6. Interpreting the Figure

(1) VQE convergence. Each curve shows the running minimum of the energy error $E(\boldsymbol{\theta})-E_{\text{exact}}$ on a logarithmic axis. By the variational principle this error can never be negative, so the curves decrease monotonically toward zero. The orange dashed line marks chemical accuracy, $1.6\times10^{-3}$ Ha, the threshold at which a quantum-chemistry result becomes practically useful. A run that crosses this line and keeps falling has found the correct basin. Runs that stall higher have been trapped in a local minimum, which is exactly why several restarts are used.

(2) VQE energy landscape (3D). This surface is the energy as a function of $\theta_0$ and $\theta_1$ with the other two angles frozen. Its smooth, periodic shape reflects the trigonometric dependence of the energy on rotation angles. The red marker sits at a minimum of the surface. The picture shows why gradient-free and gradient-based optimizers work well on small, shallow ansätze: the valleys are broad and well-defined. It also hints at the difficulty of larger circuits, where many such valleys coexist.

(3) Best cut found by QAOA. The graph is drawn with its most probable bitstring from the $p=4$ state. Node colors (blue or orange) indicate which side of the partition each node belongs to. Green solid edges are cut, and gray dashed edges are not. A good partition alternates colors around the ring and uses the two chords effectively. Because the cost is symmetric under flipping all bits, the complementary coloring is exactly as probable.

(4) QAOA landscape at $p=1$ (3D). The surface of $\langle C\rangle$ over $(\gamma,\beta)$ shows a periodic ridge-and-valley structure. The cyan marker is the best $p=1$ point found by the optimizer, and it sits on one of the highest ridges. The several equivalent peaks come from the symmetries of the problem, namely the periodicity of the mixer in $\beta$ and the integer-valued cost in $\gamma$. This explains why random restarts consistently find the same maximum value, even from different starting angles.

(5) Approximation ratio versus depth. The ratio $r_p$ should rise as $p$ increases and approach the dashed line at 1. Even $p=1$ already gives a non-trivial guarantee, and each additional layer adds two variational parameters and therefore strictly more expressive power. The gain per layer typically shrinks, which is the usual trade-off between circuit depth, which costs noise on real hardware, and solution quality.

(6) Output distribution. Each bar is the probability of measuring one of the 64 bitstrings, colored by its cut value. After optimization at $p=4$, the probability mass concentrates on bright, high-cut bitstrings, while low-cut strings are suppressed. Because measurements sample from this distribution, the key practical quantity is the total probability on optimal cuts: repeating the circuit a handful of times is usually enough to see a maximum cut at least once.

7. Conclusion

The two examples show the same pattern. A short parameterized circuit prepares a state, a cost function is evaluated, and a classical optimizer improves the parameters.

VQE turned the variational principle into a practical minimization and reproduced the exact ground-state energy of H$_2$ with a four-parameter circuit.
QAOA encoded a combinatorial problem into a phase operator and a mixer, and its solution quality improved with circuit depth $p$.

The same structure scales to larger molecules and graphs. The hard parts there are noise, barren plateaus, and the cost of measurement, but the hybrid quantum-classical loop is unchanged.

Squeezing Quantum Circuits

October 18, 2026

Exact CNOT and Depth Minimization with Breadth-First Search

Every gate on a real quantum device adds noise, and every extra layer of depth gives qubits more time to decohere. Quantum circuit compilation and synthesis is the craft of rewriting a circuit into an equivalent one that is as small and shallow as possible. In its exact form it is a discrete, combinatorial optimization problem.

In this post we solve a complete, concrete instance. We synthesize every linear reversible circuit on 4 qubits (all 20,160 of them) with the provably minimum number of CNOT gates and the provably minimum depth. We then compare the optimum against a standard heuristic, Gauss–Jordan elimination.

1. The problem

A circuit built only from CNOT gates acts linearly on computational basis states over the field $\mathbb{F}_2$:

$$
|x_1, \dots, x_n\rangle ;\mapsto; |M\mathbf{x}\rangle, \qquad M \in GL(n, \mathbb{F}_2).
$$

A $\mathrm{CNOT}(c \to t)$ gate is the elementary row operation

$$
x_t \leftarrow x_t \oplus x_c .
$$

Applying the gates $g_1, \dots, g_k$ in sequence multiplies the matrix by the corresponding elementary matrices:

$$
M(C) = E_{g_k} \cdots E_{g_2} E_{g_1}.
$$

For $n=4$ the number of distinct circuits is the order of the group

Problem A (gate count). Given a target matrix $M$, find a gate sequence $C$ solving

$$
\min_{C} ; |C| \quad \text{subject to} \quad M(C) = M .
$$

Problem B (depth). Let $L$ be a layer of CNOTs acting on pairwise disjoint qubits, and let $\mathcal{L}$ be the set of all such layers. Find

$$
\min_{(L_1,\dots,L_d)} ; d \quad \text{subject to} \quad M(L_d)\cdots M(L_1) = M, \quad L_i \in \mathcal{L}.
$$

Both are shortest-path problems on the Cayley graph of $GL(4,\mathbb{F}_2)$. The vertices are the 20,160 matrices, and the edges are CNOT gates (Problem A) or parallel CNOT layers (Problem B). Because the graph is small, breadth-first search (BFS) from the identity gives the exact optimum for every target matrix at once.

As a baseline, we implement the classical Gauss–Jordan heuristic, which reduces $M$ to the identity row by row. Its worst-case gate count is $O(n^2)$, but it is not optimal.

2. Source code

import itertools
import time
from collections import deque, Counter

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

N = 4
IDENT = tuple(1 << i for i in range(N))
CNOTS = [(c, t) for c in range(N) for t in range(N) if c != t]


def apply_cnot(M, c, t):
    rows = list(M)
    rows[t] ^= rows[c]
    return tuple(rows)


def apply_layer(M, layer):
    for c, t in layer:
        M = apply_cnot(M, c, t)
    return M


def build_layers():
    layers = [(g,) for g in CNOTS]
    for g1, g2 in itertools.combinations(CNOTS, 2):
        if len({*g1, *g2}) == 4:
            layers.append((g1, g2))
    return layers


GATE_MOVES = [(g,) for g in CNOTS]
LAYER_MOVES = build_layers()


def bfs(moves):
    dist = {IDENT: 0}
    parent = {IDENT: None}
    queue = deque([IDENT])
    while queue:
        M = queue.popleft()
        d = dist[M] + 1
        for mv in moves:
            M2 = apply_layer(M, mv)
            if M2 not in dist:
                dist[M2] = d
                parent[M2] = (M, mv)
                queue.append(M2)
    return dist, parent


def reconstruct(parent, M):
    path = []
    while parent[M] is not None:
        M, mv = parent[M]
        path.append(mv)
    return path[::-1]


def greedy_synth(M):
    rows = list(M)
    ops = []
    for i in range(N):
        if not (rows[i] >> i) & 1:
            for j in range(i + 1, N):
                if (rows[j] >> i) & 1:
                    rows[i] ^= rows[j]
                    ops.append((j, i))
                    break
        for r in range(N):
            if r != i and (rows[r] >> i) & 1:
                rows[r] ^= rows[i]
                ops.append((i, r))
    assert tuple(rows) == IDENT
    return ops[::-1]


def simulate(gates):
    M = IDENT
    for c, t in gates:
        M = apply_cnot(M, c, t)
    return M


def circuit_depth(gates):
    level = [0] * N
    for c, t in gates:
        l = max(level[c], level[t]) + 1
        level[c] = level[t] = l
    return max(level) if gates else 0


def draw_circuit(ax, gates, title):
    n_g = len(gates)
    for q in range(N):
        ax.hlines(q, -0.7, n_g - 0.3, color="black", lw=1)
    for k, (c, t) in enumerate(gates):
        ax.vlines(k, min(c, t), max(c, t), color="tab:blue", lw=2)
        ax.plot(k, c, "o", color="tab:blue", ms=8)
        ax.plot(k, t, marker="o", mfc="white", mec="tab:blue", ms=16, mew=2)
        ax.plot(k, t, marker="+", color="tab:blue", ms=12, mew=2)
    ax.set_ylim(N - 0.5, -0.5)
    ax.set_xlim(-0.8, n_g - 0.2)
    ax.set_yticks(range(N))
    ax.set_yticklabels([f"q{q}" for q in range(N)])
    ax.set_xticks([])
    ax.set_title(title, fontsize=11)
    for s in ("top", "right", "bottom"):
        ax.spines[s].set_visible(False)


# ---------------------------------------------------------------
# 1) Exhaustive BFS: optimal gate count and optimal depth
# ---------------------------------------------------------------
t0 = time.time()
dist_gate, parent_gate = bfs(GATE_MOVES)
dist_depth, parent_depth = bfs(LAYER_MOVES)
t_bfs = time.time() - t0

all_mats = list(dist_gate.keys())
assert len(all_mats) == 20160

# ---------------------------------------------------------------
# 2) Compare with the greedy (Gauss-Jordan) heuristic
# ---------------------------------------------------------------
opt_g, opt_d, gr_g, gr_d = [], [], [], []
for M in all_mats:
    g = greedy_synth(M)
    assert simulate(g) == M
    opt_g.append(dist_gate[M])
    opt_d.append(dist_depth[M])
    gr_g.append(len(g))
    gr_d.append(circuit_depth(g))
opt_g, opt_d, gr_g, gr_d = map(np.array, (opt_g, opt_d, gr_g, gr_d))

# Verify that reconstructed optimal circuits really implement the targets
for M in all_mats[:200]:
    assert simulate([g for mv in reconstruct(parent_gate, M) for g in mv]) == M
    assert simulate([g for mv in reconstruct(parent_depth, M) for g in mv]) == M

print(f"BFS finished in {t_bfs:.2f} s  |  |GL(4,2)| = {len(all_mats)}")
print(f"Optimal CNOT count : max = {opt_g.max()}, mean = {opt_g.mean():.3f}")
print(f"Greedy  CNOT count : max = {gr_g.max()}, mean = {gr_g.mean():.3f}")
print(f"Optimal depth      : max = {opt_d.max()}, mean = {opt_d.mean():.3f}")
print(f"Greedy  depth      : max = {gr_d.max()}, mean = {gr_d.mean():.3f}")
print(f"Greedy is already optimal for {np.mean(gr_g == opt_g) * 100:.1f}% of all circuits")
print()
print("Optimal gate count | #circuits | mean greedy count")
for k in range(opt_g.max() + 1):
    mask = opt_g == k
    print(f"{k:>18} | {mask.sum():>9} | {gr_g[mask].mean():.3f}")

# Worst-case example: largest gap between greedy and optimal
gap = gr_g - opt_g
idx = int(np.argmax(gap))
M_star = all_mats[idx]
opt_circ = [g for mv in reconstruct(parent_gate, M_star) for g in mv]
gr_circ = greedy_synth(M_star)
print()
print("Worst-case example (rows of M as bit masks):", M_star)
print(f"Optimal circuit ({len(opt_circ)} CNOTs, depth {circuit_depth(opt_circ)}):", opt_circ)
print(f"Greedy  circuit ({len(gr_circ)} CNOTs, depth {circuit_depth(gr_circ)}):", gr_circ)

# ---------------------------------------------------------------
# 3) All figures in a single output
# ---------------------------------------------------------------
fig = plt.figure(figsize=(21, 12))
gs = GridSpec(2, 3, figure=fig, hspace=0.35, wspace=0.28)

ax1 = fig.add_subplot(gs[0, 0])
bins = np.arange(0, max(opt_g.max(), gr_g.max()) + 2) - 0.5
ax1.hist([opt_g, gr_g], bins=bins, label=["Optimal (BFS)", "Greedy (Gauss-Jordan)"],
         color=["tab:green", "tab:red"], rwidth=0.85)
ax1.set_xlabel("CNOT count")
ax1.set_ylabel("Number of circuits")
ax1.set_title("Gate count distribution over all 20160 circuits")
ax1.legend()

ax2 = fig.add_subplot(gs[0, 1])
H = np.zeros((opt_g.max() + 1, gr_g.max() + 1), dtype=int)
for a, b in zip(opt_g, gr_g):
    H[a, b] += 1
im = ax2.imshow(np.log10(H + 1), origin="lower", cmap="viridis")
for i in range(H.shape[0]):
    for j in range(H.shape[1]):
        if H[i, j] > 0:
            ax2.text(j, i, str(H[i, j]), ha="center", va="center", color="white", fontsize=7)
ax2.set_xlabel("Greedy CNOT count")
ax2.set_ylabel("Optimal CNOT count")
ax2.set_title("Optimal vs greedy (cell = #circuits)")
fig.colorbar(im, ax=ax2, label=r"$\log_{10}(1+\mathrm{count})$")

ax3 = fig.add_subplot(gs[0, 2], projection="3d")
pair_counts = Counter(zip(opt_g.tolist(), opt_d.tolist()))
xs = np.array([k[0] for k in pair_counts], dtype=float)
ys = np.array([k[1] for k in pair_counts], dtype=float)
zs = np.array(list(pair_counts.values()), dtype=float)
colors = plt.cm.plasma(zs / zs.max())
ax3.bar3d(xs - 0.4, ys - 0.4, np.zeros_like(zs), 0.8, 0.8, zs, color=colors, shade=True)
ax3.set_xlabel("Optimal CNOT count")
ax3.set_ylabel("Optimal depth")
ax3.set_zlabel("Number of circuits")
ax3.set_title("Gate count vs depth (3D)")
ax3.view_init(elev=28, azim=-60)

ax4 = fig.add_subplot(gs[1, 0])
ks = np.arange(opt_g.max() + 1)
mean_gr = np.array([gr_g[opt_g == k].mean() for k in ks])
ax4.plot(ks, ks, "g--", label="Optimal")
ax4.plot(ks, mean_gr, "ro-", label="Greedy (mean)")
ax4.fill_between(ks, ks, mean_gr, color="red", alpha=0.15)
ax4.set_xlabel("Optimal CNOT count")
ax4.set_ylabel("CNOT count")
ax4.set_title("Average overhead of the greedy heuristic")
ax4.legend()
ax4.grid(alpha=0.3)

ax5 = fig.add_subplot(gs[1, 1])
draw_circuit(ax5, opt_circ, f"Optimal: {len(opt_circ)} CNOTs, depth {circuit_depth(opt_circ)}")
ax6 = fig.add_subplot(gs[1, 2])
draw_circuit(ax6, gr_circ, f"Greedy: {len(gr_circ)} CNOTs, depth {circuit_depth(gr_circ)}")

plt.show()

3. Code walkthrough

Representing a circuit as a tuple of integers

Each row of $M$ is stored as a 4-bit integer, and the whole matrix is a tuple of four integers. This makes the matrix hashable, so it can be a dictionary key, and a CNOT becomes a single XOR:

1	rows[t] ^= rows[c]

The identity matrix is (1, 2, 4, 8). Representing matrices as bit masks instead of NumPy arrays is the main speed optimization. Each BFS expansion costs a few integer operations rather than a matrix multiplication plus conversion to a hashable key, so both exhaustive searches finish in a fraction of a second. No further acceleration is needed.

Two search spaces, one BFS

The function bfs(moves) is a plain breadth-first search from the identity. The only difference between the two optimization problems is the move set:

GATE_MOVES contains the 12 single CNOT gates $(c,t)$ with $c \neq t$. BFS distance equals the minimum gate count.
LAYER_MOVES contains the 12 single gates plus all pairs of CNOTs acting on four distinct qubits, which can run simultaneously. There are 12 such pairs: 3 ways to split four qubits into two pairs, times 2 orientations for each pair. BFS distance equals the minimum depth.

Because every edge has unit cost, the first time BFS reaches a matrix is guaranteed to be via a shortest path. The parent dictionary stores a back-pointer, so reconstruct recovers an explicit optimal circuit for any target.

Greedy baseline

greedy_synth is the Gauss–Jordan procedure. For each column $i$ it makes sure the pivot $M_{ii}$ is 1, swapping in a row below through an XOR if necessary, and then clears every other 1 in that column. Each row operation is a CNOT. The recorded operations reduce $M$ to the identity, so reversing the list gives a circuit that builds $M$. The function asserts the reduction succeeded.

Self-verification

The code is written to be self-checking:

assert len(all_mats) == 20160 confirms BFS reached exactly the order of the group.
simulate replays any gate list, and the assertions confirm each synthesized circuit reproduces its target matrix.
The reconstructed optimal circuits for the first 200 matrices are replayed in the same way, for both the gate-count and depth searches.

Depth of a gate list

circuit_depth schedules gates as early as possible. A gate placed on qubits $c$ and $t$ lands one level after the later of the two qubits’ current levels, and the circuit depth is the maximum level reached.

4. Console output

BFS finished in 1.16 s  |  |GL(4,2)| = 20160
Optimal CNOT count : max = 9, mean = 5.685
Greedy  CNOT count : max = 15, mean = 7.638
Optimal depth      : max = 6, mean = 4.173
Greedy  depth      : max = 14, mean = 7.011
Greedy is already optimal for 16.5% of all circuits

Optimal gate count | #circuits | mean greedy count
                 0 |         1 | 0.000
                 1 |        12 | 1.000
                 2 |        96 | 2.250
                 3 |       542 | 3.738
                 4 |      2058 | 5.407
                 5 |      5316 | 6.970
                 6 |      7530 | 8.197
                 7 |      4058 | 9.017
                 8 |       541 | 9.571
                 9 |         6 | 10.333

Worst-case example (rows of M as bit masks): (8, 3, 5, 9)
Optimal circuit (4 CNOTs, depth 4): [(0, 1), (0, 2), (0, 3), (3, 0)]
Greedy  circuit (14 CNOTs, depth 13): [(3, 2), (3, 1), (2, 3), (2, 1), (2, 0), (3, 2), (1, 3), (1, 2), (1, 0), (2, 1), (0, 3), (0, 2), (0, 1), (1, 0)]

5. Output figure

6. Results and discussion

The optimum is dramatically smaller than the heuristic

Over all 20,160 circuits, exhaustive search shows that no 4-qubit CNOT circuit ever needs more than 9 gates, and the average optimal circuit needs only about 5.7 gates. Gauss–Jordan elimination, in contrast, reaches 15 gates in the worst case and averages about 7.6. It matches the true optimum for only about 16.5% of the circuits.

The population by optimal gate count is a bell-shaped distribution:

Optimal CNOTs	0	1	2	3	4	5	6	7	8	9
Circuits	1	12	96	542	2058	5316	7530	4058	541	6

Almost all of the group sits at 5–7 gates, and only 6 matrices are as hard as 9 gates. The counts grow by roughly a factor of 8 per level at the start, then saturate once most of the group has been visited. This is the expected shape for BFS on a finite Cayley graph.

Reading each panel

Top left: gate-count histogram. The green bars (optimal) are concentrated between 4 and 7 gates, while the red bars (greedy) are shifted to the right and have a long tail up to 15. Every gate in that tail is wasted noise on real hardware.

Top middle: optimal vs. greedy heat map. Each cell counts the circuits with a given (optimal, greedy) pair. Any circuit on the diagonal is solved optimally by the heuristic, and everything above the diagonal is overhead. The mass lies well above the diagonal, and the gap widens for the easiest circuits. A circuit that needs only 4 CNOTs is often compiled into 5, 6 or many more.

Top right: 3D bar chart. The $x$-axis is the optimal gate count, the $y$-axis is the optimal depth, and the $z$-axis is the number of circuits. Depth can never exceed gate count, because a circuit with $k$ gates can always be scheduled in at most $k$ layers. The bars therefore lie on and below the diagonal of the $(x,y)$ plane. The tall ridge shows that the typical circuit needs about 6 gates but fewer layers, because parallelism lets disjoint CNOTs share a layer. The optimal depth never exceeds 6, and its mean is about 4.2. The greedy method reaches a depth of 14.

Bottom left: average overhead. For each optimal gate count $k$, the red curve shows the mean greedy gate count. The shaded band between the red and dashed green curves is the average overhead, and it is positive at every level above zero. For circuits of optimal size 2, the heuristic already averages 2.25 gates. For optimal size 4 it averages about 5.4 gates.

Bottom middle and right: the worst case. The program picks the circuit with the largest gap between the two methods. For the matrix with rows (8, 3, 5, 9), the exact search finds a circuit with just 4 CNOTs (depth 4), whereas Gauss–Jordan returns 14 CNOTs (depth 13). That is $3.5\times$ more gates and more than $3\times$ the depth, for the very same unitary. A blind row-by-row elimination ignores shared structure, such as one CNOT that clears several entries at once, which an exhaustive search finds immediately.

Scaling and what comes next

The approach scales poorly in the number of qubits because $|GL(n,\mathbb{F}_2)|$ grows as roughly $2^{n^2}$:

$$
|GL(4,\mathbb{F}_2)| = 20160, \qquad |GL(5,\mathbb{F}_2)| = 9{,}999{,}360, \qquad |GL(6,\mathbb{F}_2)| \approx 2\times 10^{10}.
$$

Four qubits is trivial for BFS, five qubits is feasible with a compact encoding, and six or more calls for heuristics such as meet-in-the-middle search, SAT/ILP encodings, or the Patel–Markov–Hayes algorithm, which achieves $O(n^2/\log n)$ gates. Even then, exact tables for small blocks remain useful. Compilers can partition a large circuit into 3- or 4-qubit blocks and replace each one with its provably optimal implementation from a precomputed table like the one we just built. This pattern is called peephole optimization.

7. Conclusion

We cast gate-count and depth minimization as shortest-path problems on the group $GL(4,\mathbb{F}_2)$ and solved them exactly for all 20,160 linear reversible circuits with two breadth-first searches. The key findings are:

Every 4-qubit CNOT circuit can be implemented with at most 9 CNOTs and at most 6 layers.
A standard Gauss–Jordan compiler needs up to 15 CNOTs and is optimal for only about one circuit in six.
In the worst case the heuristic is $3.5\times$ larger than the optimum (14 vs. 4 gates).

Small as it is, this example shows why circuit synthesis matters. An optimization that happens entirely at compile time can remove most of the noisy operations from a circuit before it ever runs on hardware.

Teaching a Qubit to Dance

October 17, 2026

Optimal Quantum Control with GRAPE

A quantum gate is not something you simply “apply.” On real hardware, a gate is the end result of a carefully shaped control pulse acting on a qubit that also has its own drift dynamics. Optimal quantum control (GRAPE: GRadient Ascent Pulse Engineering) tunes the control pulse so that the resulting gate fidelity is maximized.

In this article we synthesize a Hadamard gate on a qubit with a nonzero detuning, using GRAPE with exact analytic gradients, and compare it with a naive textbook pulse sequence.

The Problem

We consider a single qubit with the control Hamiltonian

$$
H(t) = \frac{\Delta}{2},\sigma_z + \frac{u_x(t)}{2},\sigma_x + \frac{u_y(t)}{2},\sigma_y ,
$$

where $\Delta = 0.8$ is a fixed detuning (the drift term we cannot switch off) and $u_x(t), u_y(t)$ are the two control amplitudes we are free to choose, with $|u_{x,y}| \le 2$.

The goal is to realize the Hadamard gate

$$
U_{\text{target}} = \frac{1}{\sqrt{2}}\begin{pmatrix} 1 & 1 \ 1 & -1 \end{pmatrix}
$$

in a total time $T = 6$.

Discretization

We split $[0, T]$ into $N = 60$ slices of width $\delta t = T/N$, with piecewise-constant controls $u_{x,k}, u_{y,k}$. Each slice has an exact propagator. Define

Then

$$
U_k = e^{-iH_k\delta t} = \cos\theta_k, I - i,\frac{\sin\theta_k}{\omega_k},\mathbf{h}_k\cdot\boldsymbol{\sigma}.
$$

The total gate is

$$
V = U_N U_{N-1}\cdots U_1 .
$$

Fidelity and Gradient

The (phase-insensitive) gate fidelity is

$$
F = \left|\frac{1}{2},\mathrm{Tr}!\left(U_{\text{target}}^{\dagger} V\right)\right|^2 = |g|^2 .
$$

Define the forward product $X_k = U_k\cdots U_1$ and the backward product $B_k = U_N\cdots U_{k+1}$. The derivative with respect to a control $u_{j,k}$ is

$$
\frac{\partial F}{\partial u_{j,k}} = 2,\mathrm{Re}!\left[,g^{*}\cdot\frac{1}{2},\mathrm{Tr}!\left(U_{\text{target}}^{\dagger}, B_k,\frac{\partial U_k}{\partial u_{j,k}}, X_{k-1}\right)\right].
$$

Every derivative costs only a few $2\times 2$ matrix products, so the full gradient is obtained at roughly the cost of one forward and one backward sweep. We hand it to L-BFGS-B together with box constraints on the amplitudes.

The Baseline

For comparison we use a pulse that ignores the detuning. It is a $\pi/2$ rotation about $y$ followed by a $\pi$ rotation about $x$, which is a Hadamard gate up to a global phase when $\Delta = 0$.

Source Code

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

np.random.seed(7)

# ---------------- Pauli matrices ----------------
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])

# ---------------- Problem settings ----------------
N = 60            # number of time slices
T = 6.0           # total gate time
DT = T / N
DELTA = 0.8       # qubit detuning (drift term)
UMAX = 2.0        # amplitude bound for each control
U_TARGET = (SX + SZ) / np.sqrt(2)   # Hadamard gate


def slice_unitaries(u, delta=DELTA, scale=1.0, dt=DT):
    """Exact 2x2 exponentials of H = (scale*ux/2) X + (scale*uy/2) Y + (delta/2) Z."""
    h = np.stack([scale * u[0] / 2, scale * u[1] / 2, np.full(u.shape[1], delta / 2)], axis=1)
    w = np.linalg.norm(h, axis=1) + 1e-14
    th = w * dt
    hs = np.einsum('ki,ijl->kjl', h, PAULI)
    U = (np.cos(th)[:, None, None] * I2
         - 1j * (np.sin(th) / w)[:, None, None] * hs)
    return U, h, w, th


def total_unitary(u, delta=DELTA, scale=1.0, dt=DT):
    U, *_ = slice_unitaries(u, delta, scale, dt)
    V = I2.copy()
    for k in range(U.shape[0]):
        V = U[k] @ V
    return V


def fidelity(u, delta=DELTA, scale=1.0):
    V = total_unitary(u, delta, scale)
    g = np.trace(U_TARGET.conj().T @ V) / 2
    return float(np.abs(g) ** 2)


def cost_and_grad(x):
    """Return 1-F and its analytic gradient (exact derivative of each slice propagator)."""
    u = x.reshape(2, N)
    U, h, w, th = slice_unitaries(u)

    # forward products X_k = U_k ... U_1
    X = [I2]
    for k in range(N):
        X.append(U[k] @ X[-1])
    # backward products B_k = U_N ... U_{k+1}
    B = [None] * (N + 1)
    B[N] = I2
    for k in range(N - 1, -1, -1):
        B[k] = B[k + 1] @ U[k]
    V = X[N]
    g = np.trace(U_TARGET.conj().T @ V) / 2
    F = np.abs(g) ** 2

    grad = np.zeros((2, N))
    for k in range(N):
        wk, thk, hk = w[k], th[k], h[k]
        f = np.sin(thk) / wk
        df = (DT * np.cos(thk) * wk - np.sin(thk)) / wk ** 2
        hs = np.einsum('i,ijl->jl', hk, PAULI)
        for j in range(2):
            dU = (-DT * np.sin(thk) * hk[j] / wk * I2
                  - 1j * (f * PAULI[j] + df * hk[j] / wk * hs))
            dU = dU * 0.5   # chain rule: h_j = u_j / 2
            dg = np.trace(U_TARGET.conj().T @ B[k + 1] @ dU @ X[k]) / 2
            grad[j, k] = 2 * np.real(np.conj(g) * dg)
    return 1 - F, -grad.ravel()


# ---------------- GRAPE optimisation (multi-start L-BFGS-B) ----------------
bounds = [(-UMAX, UMAX)] * (2 * N)
best, best_hist = None, None
print("=== GRAPE optimisation (Hadamard gate, detuning = %.2f) ===" % DELTA)
for trial in range(5):
    x0 = 0.3 * np.random.randn(2 * N)
    hist = []
    res = minimize(cost_and_grad, x0, jac=True, method='L-BFGS-B', bounds=bounds,
                   callback=lambda xk: hist.append(1 - cost_and_grad(xk)[0]),
                   options={'maxiter': 500, 'ftol': 1e-14, 'gtol': 1e-10})
    Ft = 1 - res.fun
    print("restart %d: fidelity = %.8f  (iterations: %d)" % (trial + 1, Ft, res.nit))
    if best is None or Ft > 1 - best.fun:
        best, best_hist = res, hist

u_opt = best.x.reshape(2, N)
F_opt = fidelity(u_opt)

# ---------------- Baseline: naive square pulses (ignores detuning) ----------------
u_base = np.zeros((2, N))
u_base[1, : N // 3] = (np.pi / 2) / (T / 3)       # Ry(pi/2)
u_base[0, N // 3:] = np.pi / (2 * T / 3)          # Rx(pi)
F_base = fidelity(u_base)

print("\n=== Result summary ===")
print("Naive square pulses : F = %.6f" % F_base)
print("GRAPE optimised     : F = %.8f" % F_opt)
print("Infidelity (GRAPE)  : %.2e" % (1 - F_opt))

# gradient check
xt = np.random.randn(2 * N) * 0.5
c0, g0 = cost_and_grad(xt)
idx, eps = 17, 1e-6
xp = xt.copy(); xp[idx] += eps
xm = xt.copy(); xm[idx] -= eps
num = (cost_and_grad(xp)[0] - cost_and_grad(xm)[0]) / (2 * eps)
print("Gradient check      : analytic = %.8f, numeric = %.8f" % (g0[idx], num))

# ---------------- Robustness maps ----------------
d_err = np.linspace(-0.4, 0.4, 25)
s_err = np.linspace(0.85, 1.15, 25)
DD, SS = np.meshgrid(d_err, s_err)
FG = np.zeros_like(DD)
FB = np.zeros_like(DD)
for a in range(DD.shape[0]):
    for b in range(DD.shape[1]):
        FG[a, b] = fidelity(u_opt, DELTA + DD[a, b], SS[a, b])
        FB[a, b] = fidelity(u_base, DELTA + DD[a, b], SS[a, b])
print("Fidelity at +-0.2 detuning error (GRAPE / naive): %.5f / %.5f" % (
    fidelity(u_opt, DELTA + 0.2, 1.0), fidelity(u_base, DELTA + 0.2, 1.0)))


# ---------------- State trajectories on the Bloch sphere ----------------
def bloch_path(u, sub=6):
    psi = np.array([1, 0], dtype=complex)
    pts = [psi]
    for k in range(N):
        Uk, *_ = slice_unitaries(u[:, k:k + 1], dt=DT / sub)
        for _ in range(sub):
            psi = Uk[0] @ psi
            pts.append(psi)
    pts = np.array(pts)
    bx = [np.real(p.conj() @ SX @ p) for p in pts]
    by = [np.real(p.conj() @ SY @ p) for p in pts]
    bz = [np.real(p.conj() @ SZ @ p) for p in pts]
    return np.array(bx), np.array(by), np.array(bz)


pg = bloch_path(u_opt)
pb = bloch_path(u_base)
t_edges = np.linspace(0, T, N + 1)

# ---------------- Plot (single figure) ----------------
fig = plt.figure(figsize=(20, 12))

ax1 = fig.add_subplot(2, 3, 1)
ax1.semilogy(1 - np.array(best_hist), lw=2, color='tab:blue')
ax1.set_xlabel('Iteration')
ax1.set_ylabel('Infidelity  $1-F$')
ax1.set_title('GRAPE convergence')
ax1.grid(True, which='both', alpha=0.3)

ax2 = fig.add_subplot(2, 3, 2)
ax2.step(t_edges, np.append(u_opt[0], u_opt[0, -1]), where='post', lw=2, label='$u_x$ (GRAPE)')
ax2.step(t_edges, np.append(u_opt[1], u_opt[1, -1]), where='post', lw=2, label='$u_y$ (GRAPE)')
ax2.step(t_edges, np.append(u_base[0], u_base[0, -1]), where='post', ls='--', color='gray', label='$u_x$ (naive)')
ax2.step(t_edges, np.append(u_base[1], u_base[1, -1]), where='post', ls=':', color='black', label='$u_y$ (naive)')
ax2.set_xlabel('Time')
ax2.set_ylabel('Control amplitude')
ax2.set_title('Optimised control pulses')
ax2.legend(fontsize=9)
ax2.grid(True, alpha=0.3)

ax3 = fig.add_subplot(2, 3, 3, projection='3d')
uu, vv = np.mgrid[0:2 * np.pi:40j, 0:np.pi:20j]
ax3.plot_wireframe(np.cos(uu) * np.sin(vv), np.sin(uu) * np.sin(vv), np.cos(vv),
                   color='lightgray', alpha=0.4, linewidth=0.5)
ax3.plot(*pg, color='tab:red', lw=2.5, label='GRAPE')
ax3.plot(*pb, color='tab:blue', lw=2, ls='--', label='Naive')
ax3.scatter([0], [0], [1], color='black', s=60, label='Start $|0\\rangle$')
ax3.scatter([1], [0], [0], color='green', s=80, marker='*', label='Target $|+\\rangle$')
ax3.set_xlabel('X'); ax3.set_ylabel('Y'); ax3.set_zlabel('Z')
ax3.set_title('Bloch-sphere trajectory')
ax3.legend(fontsize=9, loc='upper left')

ax4 = fig.add_subplot(2, 3, 4, projection='3d')
ax4.plot_surface(DD, SS, FG, cmap=cm.viridis, edgecolor='none', alpha=0.95)
ax4.set_xlabel('Detuning error'); ax4.set_ylabel('Amplitude scale'); ax4.set_zlabel('Fidelity')
ax4.set_zlim(0, 1)
ax4.set_title('Robustness map: GRAPE pulse')

ax5 = fig.add_subplot(2, 3, 5, projection='3d')
ax5.plot_surface(DD, SS, FB, cmap=cm.magma, edgecolor='none', alpha=0.95)
ax5.set_xlabel('Detuning error'); ax5.set_ylabel('Amplitude scale'); ax5.set_zlabel('Fidelity')
ax5.set_zlim(0, 1)
ax5.set_title('Robustness map: naive pulse')

ax6 = fig.add_subplot(2, 3, 6)
tt = np.linspace(0, T, len(pg[2]))
ax6.plot(tt, pg[0], label='$\\langle X\\rangle$ (GRAPE)', color='tab:red')
ax6.plot(tt, pg[1], label='$\\langle Y\\rangle$ (GRAPE)', color='tab:orange')
ax6.plot(tt, pg[2], label='$\\langle Z\\rangle$ (GRAPE)', color='tab:purple')
ax6.plot(tt, pb[0], '--', label='$\\langle X\\rangle$ (naive)', color='tab:blue', alpha=0.7)
ax6.axhline(1.0, color='green', lw=0.8, ls=':')
ax6.set_xlabel('Time')
ax6.set_ylabel('Bloch components')
ax6.set_title('State evolution from $|0\\rangle$')
ax6.legend(fontsize=8, ncol=2)
ax6.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

Code Walkthrough

Setup and constants

The Pauli matrices are stacked into a single array PAULI of shape $(3, 2, 2)$. This lets us build $\mathbf{h}\cdot\boldsymbol{\sigma}$ for all time slices with one einsum call. N, T, DELTA, and UMAX define the discretization, the gate time, the drift, and the amplitude bound. The Hadamard gate is written as $(\sigma_x+\sigma_z)/\sqrt{2}$.

`slice_unitaries`: the fast exponentials

Calling a generic matrix exponential 60 times per cost evaluation would be wasteful. For a $2\times2$ Hamiltonian of the form $\mathbf{h}\cdot\boldsymbol{\sigma}$, the exponential has the closed form shown earlier, so the function computes the vector $\mathbf{h}_k$, its norm $\omega_k$, and the angle $\theta_k$ for all slices at once with NumPy broadcasting. The small constant 1e-14 added to the norm prevents division by zero. The scale argument multiplies the control amplitudes, and delta replaces the detuning. We use both later to probe robustness against pulse-amplitude and detuning errors.

`total_unitary` and `fidelity`

total_unitary multiplies the slice propagators in time order, so later slices multiply from the left, which gives $V = U_N\cdots U_1$. fidelity evaluates $|\mathrm{Tr}(U_{\text{target}}^\dagger V)/2|^2$. Using the absolute value makes the measure insensitive to the global phase, which is physically irrelevant.

`cost_and_grad`: the heart of GRAPE

First we build the lists of forward products X and backward products B. After this, every slice’s contribution to the gradient needs only the local derivative of $U_k$ sandwiched between $B_k$ and $X_{k-1}$.

The derivative of $U_k$ with respect to $h_j$ follows from differentiating $\cos\theta_k$ and $\sin\theta_k/\omega_k$ using $\partial\omega/\partial h_j = h_j/\omega$. The result is the expression for dU. The factor 0.5 accounts for $h_j = u_j/2$. The function returns $1-F$ and the negated gradient, because the optimizer minimizes.

This analytic gradient is what makes the optimization fast. Finite differences would need $2\cdot 2N = 240$ extra cost evaluations per gradient, whereas the analytic version needs one sweep. The script also verifies the gradient numerically at the end, so you can see that it is correct.

Optimization loop

We run five random restarts. Each starts from small random amplitudes and runs L-BFGS-B with box bounds $[-2, 2]$. The callback records the infidelity at every iteration so we can plot the convergence curve for the best run. The tight tolerances (ftol, gtol) let the optimizer keep going until the infidelity is essentially at machine precision.

Baseline pulse

The naive sequence is built directly in the same discretized format. The first third of the time carries a $y$-drive of strength $(\pi/2)/(T/3)$ and the remaining two thirds carry an $x$-drive of strength $\pi/(2T/3)$. Since the rotation angle is amplitude times duration, these correspond to $R_y(\pi/2)$ and $R_x(\pi)$. This is a perfect Hadamard gate when $\Delta = 0$, but the detuning ruins it.

Robustness maps

We evaluate both pulses on a $25\times25$ grid in which the actual detuning is $\Delta + \epsilon_\Delta$ and the actual amplitudes are multiplied by a factor $s$. This costs 625 cheap evaluations per pulse, which takes a fraction of a second thanks to the closed-form propagators.

Bloch-sphere trajectories

bloch_path subdivides each slice into 6 sub-steps. Within a slice the Hamiltonian is constant, so the sub-steps use exactly the same propagator formula with $\delta t/6$. This gives a smooth curve of $(\langle\sigma_x\rangle,\langle\sigma_y\rangle,\langle\sigma_z\rangle)$ starting from $|0\rangle$.

One figure for everything

All six plots are placed in a single figure, so one plt.show() produces everything: convergence, pulse shapes, the Bloch trajectory in 3D, two 3D robustness surfaces, and the time evolution of the Bloch components.

Console Output

=== GRAPE optimisation (Hadamard gate, detuning = 0.80) ===
restart 1: fidelity = 1.00000000  (iterations: 12)
restart 2: fidelity = 1.00000000  (iterations: 11)
restart 3: fidelity = 1.00000000  (iterations: 12)
restart 4: fidelity = 1.00000000  (iterations: 12)
restart 5: fidelity = 1.00000000  (iterations: 12)

=== Result summary ===
Naive square pulses : F = 0.002820
GRAPE optimised     : F = 1.00000000
Infidelity (GRAPE)  : 5.33e-15
Gradient check      : analytic = -0.02026719, numeric = -0.02026719
Fidelity at +-0.2 detuning error (GRAPE / naive): 0.81664 / 0.09848

Reading the Results

Console output

All five random restarts converge to a fidelity of essentially $1$. The infidelity of the best pulse is at the level of numerical round-off, far below $10^{-10}$. This tells us that, with $T=6$ and the given amplitude bound, the Hadamard gate is exactly reachable despite the detuning, and that the optimization landscape is benign enough for the method to find a solution from every start.

The naive pulse reaches a fidelity of only about $0.003$. The detuning is large compared with the gate duration, so the uncorrected rotations are tilted far from the ideal axes, and the gate is almost orthogonal to the target. The gradient check shows the analytic and numerical derivatives agreeing to all printed digits, which confirms the derivation.

Panel 1: convergence

The infidelity drops by many orders of magnitude within roughly a dozen iterations. The semilogarithmic scale makes the final plunge to machine precision visible. This quick convergence is the payoff of the exact gradient combined with a quasi-Newton optimizer.

Panel 2: control pulses

The GRAPE pulses are smooth, irregular-looking waveforms that stay within the $\pm 2$ bounds, while the naive pulse consists of two flat blocks. The optimized shape is not intuitive. It continuously steers the qubit while pre-compensating the rotation about $z$ that the detuning would otherwise cause.

Panel 3: Bloch-sphere trajectory (3D)

The Hadamard gate maps $|0\rangle$ to $|+\rangle$, the point $(1,0,0)$ on the sphere. The GRAPE trajectory winds around the sphere and lands on the green star. The naive trajectory starts the same way but is pushed off course by the drift, so it ends far from the target. This picture shows directly why the naive gate fails.

Panels 4 and 5: robustness surfaces (3D)

The vertical axis is the fidelity, and the horizontal axes are the detuning error and the amplitude scale. The GRAPE surface has a high plateau around the nominal point $(0,1)$ with fidelity $1$. The naive surface is low everywhere. At a detuning error of $+0.2$ the GRAPE pulse still keeps a fidelity of about $0.82$, whereas the naive pulse drops to about $0.10$.

An important point is that the GRAPE pulse was optimized only for the nominal parameters. Its fidelity therefore degrades away from the nominal point, and its robustness here is a side effect of a smooth, well-conditioned solution. When high robustness is the goal, the cost function can be averaged over an ensemble of $(\epsilon_\Delta, s)$ values, which is known as robust GRAPE.

Panel 6: Bloch components over time

The GRAPE curves show $\langle\sigma_x\rangle$ rising to $1$ at the final time while $\langle\sigma_y\rangle$ and $\langle\sigma_z\rangle$ return to $0$. The dashed naive curve for $\langle\sigma_x\rangle$ never gets close to the dotted line at $1$.

Conclusion

With about 100 lines of NumPy, we implemented GRAPE with exact analytic gradients, reaching essentially perfect fidelity for a Hadamard gate on a detuned qubit in just a few iterations. The key ingredients were the closed-form $2\times2$ exponentials, the forward and backward propagator products that turn the gradient into one cheap sweep, and a quasi-Newton optimizer with amplitude bounds. The same structure extends directly to multi-qubit systems by swapping the closed-form slice exponentials for matrix exponentials and the gradient formulas for their exact or approximate derivatives.

Designing Quantum Error-Correcting Codes with Linear Programming

October 16, 2026

The Quantum Delsarte Bound and a Stabilizer Code Search

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

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

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

1. Mathematical Background

1.1 Weight enumerators of a stabilizer code

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

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

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

1.2 Quantum MacWilliams identity

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

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

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

1.3 The shadow enumerator

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

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

1.4 The linear program

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

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

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

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

We compare it with the quantum Singleton bound:

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

1.5 The code search

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

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

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

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

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

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

2. The Concrete Problem

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

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

and compare it with the Singleton bound.

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

3. Full Source Code

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

plt.style.use('dark_background')

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4. Code Walkthrough

4.1 Krawtchouk matrices

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

4.2 The LP feasibility test

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

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

4.3 Speed

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

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

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

4.4 Random stabilizer code

random_stabilizer_code adds generators one at a time:

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

No rejection sampling is needed, so every trial succeeds.

4.5 Consistency checks

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

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

4.6 The figure

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

5. Execution Results

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

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

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

6. Interpreting the Results

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

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

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

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

6.2 The gap to the Singleton bound (top middle)

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

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

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

6.4 The best code found (bottom middle)

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

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

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

6.5 The distance histogram (bottom right)

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

7. Conclusion

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

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

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

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

October 15, 2026

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

1. The problem

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

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

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

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

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

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

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

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

2. Reducing the problem to one variable

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

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

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

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

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

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

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

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

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

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

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

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

3. A concrete example

We use the following setup.

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

The program does four things.

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

4. Source code

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

fig.tight_layout()
plt.show()

5. Code walkthrough

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

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

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

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

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

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

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

6. Execution results

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

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

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

7. Interpreting the results

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

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

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

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

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

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

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

8. Conclusion

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

Quantum Values of Nonlocal Games

October 14, 2026

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

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

That number is the quantum value of the nonlocal game:

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

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

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

1. The CHSH game

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

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

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

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

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

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

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

where

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

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

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

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

An explicit optimal strategy

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

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

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

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

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

Without assuming the strategy: the see-saw method

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

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

The iteration repeats three exactly solvable steps:

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

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

2. The complete source code

Everything is in a single cell.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

plt.tight_layout()
plt.show()

3. Code walkthrough

Section 1–2: the game and the classical value

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

Section 3: the explicit quantum strategy

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

alpha_opt and beta_opt store the optimal angles derived above.

Section 4: the see-saw

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

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

Section 5: Monte Carlo simulation

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

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

Section 6: running everything

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

Section 7: the figure

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

4. Results as graphs

The six panels are in one figure.

How to read each panel

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

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

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

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

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

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

5. Console output

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

How to read the numbers

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

6. Wrap-up

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

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

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

Verifying Quantum Steering with a Semidefinite Program

October 13, 2026

Werner States, Steerable Weight, and Critical Visibility

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

1. The problem

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

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

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

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

Steerable weight

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

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

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

Critical visibility in a single SDP

Take the Werner state

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

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

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

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

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

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

2. The worked example

We study three things:

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

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

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

3. Source code

import itertools
import time

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

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

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


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


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


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


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


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


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


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


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

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

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

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

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

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


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

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

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

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

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

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

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

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

plt.tight_layout()
plt.show()

4. Walking through the code

Building blocks (sections 1 and 2)

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

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

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

SDP (A): steerable weight (section 3)

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

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

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

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

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

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

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

Settings (section 6)

The four sets of measurement directions are:

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

Experiments (sections 7 to 10)

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

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

Plotting (section 11)

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

5. Console output

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

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

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

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

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

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

Reading the console output

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

6. Graphs

Panel (a): steerable weight versus noise

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

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

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

Panel (b): critical visibility versus number of settings

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

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

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

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

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

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

7. Why the fast version matters

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

8. Conclusion

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