Demo 3：QAOA 投资组合

第 14 章 优化与数值量子算法 · 本章演示 · 代码演示总索引

可视化：双图 —— 3 资产 Markowitz QUBO

对应教程：第 14 章 优化与数值量子算法 中相关小节 · 依赖安装：pip install "unified-quantum[visualization,simulation]"

实现要点

核心：Toy Markowitz portfolio over 3 assets A, B, C with expected returns μ = (0.10, 0.15, 0.08), covariance Σ, risk aversion λ = 4.0, and a cardinality penalty that biases toward 2-asset portfolios. The quadratic utility is encoded as a QUBO H_C (negated because uniqc's qaoa_ansatz MAXIMIZES ⟨H_C⟩) and solved by p=1 QAOA. We sweep (γ, β), overlay a classical hill-climb trajectory, then histogram the top bitstrings at the optimum — the dominant sampling peaks correspond to the asset subsets QAOA suggests holding.

完整源码

demos/ch14-optimization-numerics/03_qaoa_portfolio.py

"""
ch14-demo-3: QAOA portfolio optimization (toy 3-asset mean-variance)

对应教程：chenzhaoyun.com/quantum-tutorial/ch14-optimization-numerics/constraint-satisfaction.html
核心：Toy Markowitz portfolio over 3 assets A, B, C with expected returns
     μ = (0.10, 0.15, 0.08), covariance Σ, risk aversion λ = 4.0, and a
     cardinality penalty that biases toward 2-asset portfolios. The
     quadratic utility is encoded as a QUBO H_C (negated because uniqc's
     qaoa_ansatz MAXIMIZES ⟨H_C⟩) and solved by p=1 QAOA. We sweep (γ, β),
     overlay a classical hill-climb trajectory, then histogram the top
     bitstrings at the optimum — the dominant sampling peaks correspond
     to the asset subsets QAOA suggests holding.
依赖：pip install unified-quantum[visualization]
运行：python demos/ch14-optimization-numerics/03_qaoa_portfolio.py
"""
from __future__ import annotations

from pathlib import Path

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np

from uniqc import qaoa_ansatz, submit_task, wait_for_result

OUT = Path(__file__).parent

# === 1. Toy 3-asset portfolio: expected returns & covariance ===
ASSETS = ["A", "B", "C"]
MU = np.array([0.10, 0.15, 0.08])
COV = np.array([
    [0.0100, 0.0015, -0.0008],
    [0.0015, 0.0200, 0.0010],
    [-0.0008, 0.0010, 0.0120],
])
# Risk aversion
LAMBDA = 4.0
# Cardinality constraint penalty (encourages holding ≈ 2 assets)
CARD_PEN = 0.30
CARD_TARGET = 2

# QUBO encoding: w_i ∈ {0, 1} is whether we hold asset i. Bitstring 'qb0
# qb1 qb2' (MSB first) holds asset i iff qb_i = 1.
N_QUBITS = 3

# Build H_C from the Markowitz objective. C(w) = w^T (λ·COV) w - μ^T w +
# c · (Σ w_i - 2)^2. Ising: w_i = (1 - Z_i)/2, so the per-axis Z coefficient is
#   α_i = -linear_i/2 - Σ_{j≠i} quadratic_ij / 4
# QAOA in uniqc MAXIMIZES ⟨H_C⟩, so we negate: H_C = -α_i Z_i - β_ij Z_i Z_j.
linear = np.array([
    LAMBDA * COV[i, i] - MU[i] + CARD_PEN * (1 - 2 * CARD_TARGET)
    for i in range(N_QUBITS)
])
quad = {
    (i, j): 2 * (LAMBDA * COV[i, j] + CARD_PEN)
    for i in range(N_QUBITS) for j in range(i + 1, N_QUBITS)
}
alpha = np.array([
    -linear[i] / 2 - sum(quad.get((min(i, j), max(i, j)), 0) / 4
                         for j in range(N_QUBITS) if j != i)
    for i in range(N_QUBITS)
])
beta = {ij: q / 4 for ij, q in quad.items()}

H_C = [(f"Z{i}", -alpha[i]) for i in range(N_QUBITS)]
H_C += [(f"Z{i}Z{j}", -b) for (i, j), b in beta.items()]


def portfolio_cost(bs: str) -> float:
    w = np.array([int(c) for c in bs.zfill(N_QUBITS)])
    return float(w @ (LAMBDA * COV) @ w - MU @ w
                 + CARD_PEN * (w.sum() - CARD_TARGET) ** 2)


def portfolio_dict(bs: str) -> str:
    w = np.array([int(c) for c in bs.zfill(N_QUBITS)])
    held = [ASSETS[i] for i in range(N_QUBITS) if w[i] == 1]
    return ", ".join(held) if held else "(cash)"


# === 2. QAOA evaluation ===
SHOTS = 600


def run_qaoa(gamma: float, beta: float, shots: int = SHOTS):
    c = qaoa_ansatz(H_C, p=1, gammas=np.array([gamma]), betas=np.array([beta]))
    for q in range(N_QUBITS):
        c.measure(q)
    return wait_for_result(
        submit_task(circuit=c, backend="dummy:local:simulator", shots=shots)
    ).counts


def evaluate_qaoa(gamma: float, beta: float, shots: int = SHOTS) -> float:
    counts = run_qaoa(gamma, beta, shots)
    total = sum(counts.values())
    if total <= 0:
        return 0.0
    return sum(portfolio_cost(bs) * cnt for bs, cnt in counts.items()) / total


# === 3. γ × β landscape scan + hill-climb ===
print("ch14-demo-3: QAOA portfolio optimization (3 assets)")
print(f"  assets={ASSETS}, μ={MU.tolist()}, λ={LAMBDA}, card target={CARD_TARGET}")
print(f"  shots per QAOA run: {SHOTS}")

G = np.linspace(0.10, 1.30, 11)
B = np.linspace(0.10, 0.90, 11)
landscape = np.zeros((len(G), len(B)))
for gi, g in enumerate(G):
    for bi, b in enumerate(B):
        landscape[gi, bi] = evaluate_qaoa(g, b)
print(f"  grid scan: {len(G)} x {len(B)} = {len(G) * len(B)} points")

gi0, bi0 = np.unravel_index(int(np.argmin(landscape)), landscape.shape)
g_star, b_star = float(G[gi0]), float(B[bi0])
c_star = float(landscape[gi0, bi0])
trajectory = [(g_star, b_star, c_star)]
rng = np.random.default_rng(0)
for step in range(12):
    sigma = max(0.04, 0.18 * (0.55 ** step))
    g_try = float(np.clip(g_star + rng.normal(0, sigma), 0.05, np.pi))
    b_try = float(np.clip(b_star + rng.normal(0, sigma), 0.05, np.pi))
    c_try = evaluate_qaoa(g_try, b_try)
    if c_try < c_star:
        g_star, b_star, c_star = g_try, b_try, c_try
    trajectory.append((g_star, b_star, c_star))

counts = run_qaoa(g_star, b_star, shots=4000)
hist = {bs.zfill(N_QUBITS): cnt for bs, cnt in counts.items()}
total_best = sum(hist.values())
hist_sorted = sorted(hist.items(), key=lambda kv: -kv[1])
top_states = hist_sorted[:8]

# === 5. visualisation ===
fig = plt.figure(figsize=(12.5, 5.4))
gs = fig.add_gridspec(1, 2, width_ratios=[1.05, 1.0])

ax1 = fig.add_subplot(gs[0, 0])
im = ax1.imshow(
    landscape, origin="lower", aspect="auto", cmap="magma",
    extent=[B[0] - (B[1] - B[0]) / 2, B[-1] + (B[1] - B[0]) / 2,
            G[0] - (G[1] - G[0]) / 2, G[-1] + (G[1] - G[0]) / 2],
)
cs = ax1.contour(B, G, landscape, levels=4, colors="white", linewidths=0.7, alpha=0.5)
ax1.clabel(cs, inline=True, fontsize=8, fmt="%.2f")
ax1.plot([b_star], [g_star], "*", color="#ffd700", ms=18, mec="black", mew=0.9,
         label=rf"opt $\gamma^*={g_star:.3f},\beta^*={b_star:.3f}$")
ax1.plot([t[1] for t in trajectory], [t[0] for t in trajectory],
         "-o", color="#00d4ff", ms=3, lw=1.0, alpha=0.8,
         label="classical hill-climb")
ax1.set_xlabel(r"mixer angle $\beta$")
ax1.set_ylabel(r"cost angle $\gamma$")
ax1.set_title(r"p=1 QAOA $\langle H_C\rangle$ on Markowitz QUBO (11×11 scan)")
plt.colorbar(im, ax=ax1).set_label(r"$\langle H_C\rangle$ (lower = better)")
ax1.legend(loc="upper right", fontsize=8)

ax2 = fig.add_subplot(gs[0, 1])
labels = [bs for bs, _ in top_states]
probs = [cnt / total_best for _, cnt in top_states]
costs = [portfolio_cost(bs) for bs, _ in top_states]
best_cost = min(portfolio_cost(f"{i:03b}") for i in range(8))
colors = ["#2ca02c" if abs(c - best_cost) < 1e-6 else
          "#1f77b4" if c < 0 else "#d62728" for c in costs]
ax2.bar(range(len(labels)), probs, color=colors, ec="black", lw=0.5)
for k, (lab, p, c) in enumerate(zip(labels, probs, costs)):
    ax2.text(k, p + 0.012, f"{lab}\n[{portfolio_dict(lab)}]\ncost={c:.3f}",
             ha="center", va="bottom", fontsize=7)
ax2.set_xticks(range(len(labels)))
ax2.set_xticklabels(labels, rotation=0, fontsize=8)
ax2.set_ylabel("sampling probability")
ax2.set_title(rf"top bitstrings at $(\gamma^*,\beta^*)$ — $\langle H_C\rangle$={c_star:.3f}")
ax2.set_ylim(0, max(probs) * 1.55)
ax2.grid(alpha=0.3, axis="y")

plt.suptitle("ch14 demo 3: QAOA portfolio optimization (3-asset Markowitz)", y=1.02)
plt.tight_layout()
out_path = OUT / "03_qaoa_portfolio.png"
plt.savefig(out_path, dpi=120, bbox_inches="tight")
plt.close()

# === 6. stdout 关键统计 ===
print(f"\n  total QAOA circuit evaluations: "
      f"{len(G) * len(B) + len(trajectory) - 1 + 1}")
print(f"  grid best      : cost={float(landscape.min()):.4f} "
      f"at γ={G[gi0]:.3f}, β={B[bi0]:.3f}")
print(f"  after hill-climb: γ*={g_star:.4f}, β*={b_star:.4f}")
print(f"  mean cost       : {c_star:.4f}")

print(f"  classical ground-truth costs over all 8 portfolios:")
for i in range(8):
    bs = f"{i:03b}"
    print(f"    |{bs}>  [{portfolio_dict(bs)}]  cost={portfolio_cost(bs):+.4f}")

print(f"  top 4 QAOA samples:")
for bs, cnt in hist_sorted[:4]:
    print(f"    |{bs}>  [{portfolio_dict(bs)}]  p={cnt/total_best:.4f}  "
          f"cost={portfolio_cost(bs):+.4f}")
print("  PNG ->", out_path)

运行结果

运行方式

# 仓库根目录执行；输出 PNG 到 demos/ch14-optimization-numerics/03_qaoa_portfolio.png，关键统计到 stdout
.venv/bin/python3 demos/ch14-optimization-numerics/03_qaoa_portfolio.py

---

← 本章演示 · ↑ 代码演示总索引 · ⌂ 教程门户