变分量子算法的优化与贫瘠高原

1. 参数化量子电路（ansatz）U(θ)
   把初态映射为 |ψ(θ)⟩ = U(θ)|0...0⟩，θ 有 m 个可调参数
2. 代价函数 C(θ)
   通常为期望值 C(θ) = ⟨ψ(θ)| H |ψ(θ)⟩
   或分类损失、保真度等更一般的目标
3. 经典优化器
   迭代 θ ← θ - η∇C(θ)，梯度可由量子设备估计

迭代循环：
  初始化 θ_0
  repeat:
      量子设备：制备 |ψ(θ)⟩，测量若干次 → 估计 C(θ) 或 ∇C(θ)
      经典设备：根据估计更新 θ
  until 收敛
→ 量子设备只做「前向传播」，经典设备做「反向传播」

∂C/∂θ_j = [ C(θ_j + π/2) - C(θ_j - π/2) ] / 2
推导要点：
  U(θ) = cos(θ/2)·I - i·sin(θ/2)·P，∂U/∂θ = -(i/2)·P·U
  代入 C(θ) = ⟨0|U†HU|0⟩ 展开，三角恒等式给出上式
→ 无需额外线路，只需把该参数 ±π/2 各跑一次

# PennyLane：参数移位规则是默认梯度
import pennylane as qml
from pennylane import numpy as np

dev = qml.device("default.qubit", wires=2)

@qml.qnode(dev, diff_method="parameter-shift")
def cost(params):
    qml.RY(params[0], wires=0)
    qml.CNOT(wires=[0, 1])
    qml.RY(params[1], wires=1)
    return qml.expval(qml.PauliZ(0) @ qml.PauliZ(1))

params = np.array([0.3, 0.7], requires_grad=True)
print(qml.grad(cost)(params))   # 内部用参数移位规则

形式化：对随机初始化的 θ，
  E[∂C/∂θ_j] = 0（无偏，代价函数在随机点处近似平坦）
  Var[∂C/∂θ_j] ∝ exp(-c·n)，c > 0 为常数
→ 梯度信号幅度 ~ exp(-c·n/2)
  要分辨它，需采样次数 ~ exp(c·n)
→ 优化器「看到」一片没有信息的平原

关键结果（McClean 等，2018）：
  若 ansatz 在参数分布上形成 2-设计，且 H 的展开不含单位算符项
  则 E_θ[∂C/∂θ] = 0，且
     Var_θ[∂C/∂θ] = f(n) · ||H||^2 / (2^(2n) - 1)
→ 方差以 2^(-2n) 衰减，即梯度幅度 ~ 2^(-n)

直觉：单参数 θ_j 只影响一个小角度旋转，其贡献被线路其余
      部分的随机化稀释；2-设计下任何局部算子的期望都以
      1/dim = 2^(-n) 的因子收缩

全局可观测量（Z^⊗n）：梯度方差 ~ 2^(-2n) → 指数衰减
局部可观测量（单比特 Z_1）：梯度方差 ~ 常数 → 与 n 无关
→ 这是「局部代价函数」缓解贫瘠高原的数学依据

成因一：随机电路的 2-设计性
  深层随机电路在参数分布上接近 Haar 随机（2-设计）
  → 输出态在希尔伯特空间均匀分布，区分性信号被抹平

成因二：纠缠增长与线路深度
  层数 L 越大、纠缠越强 → 越接近 2-设计
  存在临界深度 L*，超过后梯度方差指数衰减

成因三：全局可观测量与全局代价函数
  C = ⟨Z_1 ⊗ ... ⊗ Z_n⟩ 需估计全部比特的关联
  → 方差 ~ 2^(-2n)，最直接的「放大器」

成因四：硬件噪声
  噪声本身也贡献随机化，即使 ansatz 理论上不是 2-设计
  → 噪声诱导的贫瘠高原是 NISQ 上更现实的威胁

COBYLA：约束优化的线性近似，对噪声鲁棒、广泛用于 VQE
        缺点：收敛慢，参数多时表现差
Nelder-Mead：单纯形法，实现简单，但高维退化、易早停
SPSA（同时扰动随机近似）：
  每次迭代只扰动所有参数的一个随机方向
    θ± = θ ± c_k · Δ_k（Δ_k 为随机 ±1 向量）
    梯度估计 = [C(θ+) - C(θ-)] / (2 c_k) · Δ_k^(-1)
  优点：每次迭代仅需 2 次评估，与参数数无关
  缺点：收敛慢、噪声敏感

from scipy.optimize import minimize
import numpy as np

def cost_fn(params):
    return float(cost(params))   # cost 为上面的 PennyLane qnode

res = minimize(cost_fn, np.random.uniform(0, 2 * np.pi, 2),
               method="COBYLA", options={"maxiter": 200, "tol": 1e-6})
print(res.fun, res.x)

贫瘠高原：梯度处处极小（全局平坦）→ 无法优化
局域极小：局部平坦但存在下降方向 → 可能靠重启逃脱
两者症状相似（都表现为优化停滞），但应对不同

import numpy as np
from qiskit import QuantumCircuit
from qiskit.circuit import ParameterVector
from qiskit.quantum_info import SparsePauliOp
from qiskit.primitives import StatevectorEstimator

n = 4
theta = ParameterVector("θ", 2 * n)

def ansatz(params):
    qc = QuantumCircuit(n)
    for i in range(n):
        qc.ry(params[i], i)            # 单比特旋转层
    for i in range(n - 1):
        qc.cx(i, i + 1)                # 纠缠层
    for i in range(n):
        qc.ry(params[n + i], i)
    return qc

H = SparsePauliOp("Z" * n)             # 全局可观测量
estimator = StatevectorEstimator()

def cost(params):
    qc = ansatz(params)
    return estimator.run([(qc, H)]).result()[0].data.evs

def param_shift_grad(cost, params, j, shift=np.pi / 2):
    """参数移位规则：对第 j 个参数做 ±π/2 移位"""
    p_plus = np.array(params, dtype=float)
    p_minus = np.array(params, dtype=float)
    p_plus[j] += shift
    p_minus[j] -= shift
    return 0.5 * (cost(p_plus) - cost(p_minus))

# 观察梯度方差随 n 的衰减（贫瘠高原实验）
for nn in [2, 4, 6, 8]:
    grads = [param_shift_grad(cost, np.random.uniform(0, 2 * np.pi, 2 * nn), 0)
             for _ in range(50)]
    print(f"n={nn}: Var[grad] = {np.var(grads):.3e}")
# 预期：方差随 n 指数下降 → 半对数图上近似直线

import pennylane as qml
from pennylane import numpy as np

def make_cost(n, local=True):
    dev = qml.device("default.qubit", wires=n)

    @qml.qnode(dev, diff_method="parameter-shift")
    def cost(params):
        for i in range(n):
            qml.RY(params[i], wires=i)
        for i in range(n - 1):
            qml.CNOT(wires=[i, i + 1])
        for i in range(n):
            qml.RY(params[n + i], wires=i)
        if local:      # 局部代价：梯度方差与 n 无关
            return sum(qml.expval(qml.PauliZ(i)) for i in range(n)) / n
        return qml.expval(qml.PauliZ(0) @ qml.PauliZ(n - 1))  # 全局：衰减

    return cost

for n in [4, 6, 8, 10]:
    cl, cg = make_cost(n, True), make_cost(n, False)
    vl = np.var([qml.grad(cl)(np.random.uniform(0, 2 * np.pi, 2 * n))[0]
                 for _ in range(40)])
    vg = np.var([qml.grad(cg)(np.random.uniform(0, 2 * np.pi, 2 * n))[0]
                 for _ in range(40)])
    print(f"n={n}: local Var={vl:.3e}  global Var={vg:.3e}")