如果给量子计算的应用场景排个名组合优化一定排在最前面。原因很简单日常遇到的排产、路径规划、芯片布线、社交网络划分几乎都能塞进组合优化这个筐里而这些问题一旦规模变大就极其难算。QAOAQuantum Approximate Optimization Algorithm量子近似优化算法就是冲着这个方向来的2014年由Edward Farhi等人提出之后一直是NISQ时代研究最广泛的量子算法之一。这篇文章我会从组合优化为什么难入手讲清楚QAOA的物理原理和数学推导再用手写代码的方式把基于Qiskit的QAOA完整跑通。适合两类人看懂一点量子计算但没写过QAOA的以及熟悉组合优化但想看看量子算法到底怎么落地的。1. 组合优化问题的现实压力QAOA出现的逻辑起点1.1 一个典型的MaxCut问题长什么样先把问题具体化。设想一张无向图顶点代表社交网络里的用户边代表好友关系。你想把用户分成两个社群让被切断的好友关系数量尽可能多。为什么切边多就好因为一个理想社群内部连接紧密、外部连接稀疏切割掉的外部连接越多两个群体的边界越清晰。这个直觉性问题就是最大割问题MaxCut。形式化一点图G(V,E)对每个顶点i标记一个变量x_i∈{1,-1}表示它属于左边还是右边。一条边(i,j)被切割当且仅当x_i≠x_j。那么切割数可以写成C(x) 1/2 * Σ_{(i,j)∈E} (1 - x_i x_j)x_i x_j-1时这一项得1x_i x_j1时这一项得0再乘上1/2正好是切割边数。目标就是最大化C(x)。以4个顶点的环形图为例顶点0-1-2-3-0连成一个环。最优切割一定是按二部图分开即{0,1}在一组、{2,3}在另一组四条边全被切断C4。看起来很简单但换成一般图就不一样了。MaxCut是Karp列出的21个NP完全问题之一言下之意没有人能找到多项式时间的精确算法。暴力枚举所有2^n种划分n50时就超过千万亿次n100时彻底不可行。1.2 经典算法的天花板与量子路线的机会面对NP-hard工程界的标准动作是退而求其次近似算法或启发式搜索。近似算法里最著名的是Goemans-Williamson算法用半定规划松弛在理论保证上能达到0.878的近似比。就是说算法给出的切割数至少是最优解的87.8%。这个数字听起来不错但存在几个麻烦半定规划求解本身开销不小大规模问题上跑不动工业界的组合优化往往带各种复杂约束SDP松弛很难覆盖而且0.878这个常数被UGC猜想唯一游戏猜想认为是不可突破的你要想再提升一点近似比可能就要推翻一个基础猜想。启发式算法更灵活模拟退火、遗传算法、禁忌搜索都是工程常用工具。它们能处理大实例但没有性能保证结果好坏要看参数和运气。这给了我这类做优化的人一种感觉经典路线在一个精度-效率-通用性的三角里很难同时照顾好三个角。QAOA的出现给了另一条路。它不试图设计更好的松弛或搜索策略而是把组合优化问题映射到量子系统的能量景观上用一台可以执行浅层量子线路的设备以变分的方式逼近最优解。它在线路深度上只和问题的规模边数、层数p多项式相关而不是和希尔伯特空间维度指数相关。即使p很小QAOA在某些问题上也能给出非平凡的近似保证。放在NISQ硬件的框架下看这正是当前技术条件下少有的、能直接跑在通用量子计算机上的优化算法。2. 绝热定理到参数化线路QAOA的物理内核2.1 绝热演化的朴素想法QAOA的思想根源是绝热量子计算所以得先从绝热定理讲起。量子系统演化遵循薛定谔方程。如果系统哈密顿量随时间缓慢变化且初态处于某个瞬时本征态那么系统会一直停留在对应的瞬时本征态上不会跳到别的能级。特别是如果从基态出发、缓慢演化就能始终保持基态。把这个问题和优化联系起来设计一个终态哈密顿量H_C它的基态编码了我们要找的最优解。再找一个容易制备基态的初始哈密顿量H_B。然后构建含时哈密顿量H(s) (1-s)H_B sH_C, s从0到1从H_B的基态出发把s缓慢从0推到1系统就会演化到H_C的基态附近。最后测量会以高概率得到最优解。这就是绝热量子计算的框架。问题出在缓慢两个字上。绝热定理要求演化时间T远大于系统最小能隙平方的倒数。很多NP-hard问题在中途会出现能隙几乎闭合的点T随之指数增长退相干早就把量子态毁掉了。所以纯粹绝热计算在真实硬件上很难落地量子退火机比如D-Wave能做的其实也是受限的“模拟绝热”。2.2 QAOA把绝热过程剪短成可优化的变分线路Farhi等人的思路非常聪明既然完整绝热过程太慢那就只做p步离散化并且别把每一步的时间当成固定值而是当成可以训练的参数。具体来说把H(s)的时间演化拆成p层每层先演化H_C一段时间γ_l再演化H_B一段时间β_l。整条线路作用在初态|⟩^n上得到|ψ(γ,β)⟩ e^{-iβ_p H_B} e^{-iγ_p H_C} ··· e^{-iβ_1 H_B} e^{-iγ_1 H_C} |⟩^n这里γ(γ_1...γ_p)和β(β_1...β_p)就是QAOA的变分参数。和绝热过程不同这些参数不需要满足任何慢演化条件完全可以由经典优化器自由调整目的只有一个让关于成本哈密顿量的期望值F(γ,β)⟨ψ|H_C|ψ⟩最大。当p→∞时如果参数选择得均匀QAOA可以逼近绝热演化但QAOA真正给人信心的是另一个方向哪怕p很小比如p1或p2在特定问题上也已经有可证明的性能保证并且线路非常浅能在NISQ设备上跑。p在这里通俗讲就是QAOA的迭代深度p越大表达空间越强但线路越深、经典参数优化越难。2.3 两个幺正算子到底在干什么成本层与混合层的分工理解QAOA关键在理解每一层里两个算子的角色别把它们当成黑盒。成本层U_C(γ)e^{-iγH_C}做的事情是让量子态中每个计算基分量|z⟩按照它的成本值C(z)获得一个相位e^{-iγC(z)}。在这个阶段量子态本身不改变振幅只是给好的分量和坏的分量加了不同的相位。直观理解它像一个向导把目标函数的几何信息写进了量子态的相位里。混合层U_B(β)e^{-iβH_B}其中H_BΣX_i是横场项。这一层会在不同计算基之间产生转移让各个基态的振幅重新分配。它像一个搜索器负责把振幅从相位上读出的信息转化为真正的概率流动。用登山类比成本层告诉你哪个方向是下坡但只告诉你方向至少走几步怎么分配精力得靠混合层来扰动和扩散。两者交替作用量子态就在目标函数的景观里不断演化。经典优化器通过调整γ和β决定每一步听向导的多一点还是多探索一点。这种结构和经典模拟退火中的降温-扰动循环有神似之处但由于量子系统存在相干叠加搜索方式本质上是并行的这也是QAOA在理论上可能超越经典启发式的原因。3. 从MaxCut到成本哈密顿量QAOA的第一份推导作业3.1 用Pauli Z矩阵重写切割数量子线路只能操作量子比特所以必须把MaxCut的目标函数翻译成量子算符。核心做法是把经典的二值变量x_i换成Pauli Z矩阵的本征值。Z|0⟩|0⟩、Z|1⟩-|1⟩正好和x_i∈{1,-1}同构。于是切割数C(x)变成量子力学算符H_C 1/2 Σ_{(i,j)∈E} (I - Z_i Z_j)对任意计算基态|z⟩也就是一个比特串的量子版本有⟨z|H_C|z⟩C(z)。也就是说H_C的本征值就是对应的切割数本征值越大表示切割数越多。那么MaxCut就等价于制备H_C的一个高能本征态——注意是高能态因为我们的目标C取正方向的最大值。很多读QAOA论文的人会在这一步犯迷糊通常量子算法都在找基态最低能量怎么这里要找高能态了其实只是符号约定问题。如果你想把MaxCut写成找基态的问题只需要定义H_C -H_C 1/2 Σ(Z_i Z_j - I)然后求基态。两者的数学本质完全一样。我习惯用H_C -H_C的版本做推导但为了直观起见下面代码和本文统一用最大化⟨H_C⟩来写。3.2 成本层的门实现与RZZ符号陷阱现在要处理U_C(γ)e^{-iγH_C}。由于H_C是各项之和、且各项互相对易因为都是互不重叠或共享比特的对角算符矩阵乘法顺序不影响结果指数可以拆开成乘积e^{-iγH_C} ∏_{(i,j)∈E} e^{-iγ (I - Z_i Z_j)/2}对每个边单独分析这一项e^{-iγ (I-Z_iZ_j)/2} e^{-iγ/2} · e^{iγ Z_iZ_j/2}这里e^{-iγ/2}是一个全局相位对所有计算基分量都相同测量时不影响概率分布可以安全丢掉。真正的作用是e^{iγ Z_iZ_j/2}。回忆Qiskit中RZZ门的定义RZZ(θ) e^{-iθ Z_iZ_j/2}所以e^{iγ Z_iZ_j/2} RZZ(-2γ)。这就是为什么QAOA的MaxCut线路中每条边要加一个RZZ(-2γ)门。我在不少教程和开源代码里看到有人写RZZ(2γ)、RZZ(γ)甚至RZ的变体经常其实是同一个算法用了不同的目标函数符号或不同的参数约定。只要你和自己手上的目标函数方向严格对得上就行。但为了少踩坑我的经验是第一次实现时先用一个最简单的比特串比如全部为0的比特串手工算一遍期望值确认符号正确再跑优化否则你会发现优化器收敛到一个完全无意义的参数上还很难排查。3.3 混合层与初始态的配合混合层H_BΣX_i同样是各自独立项的求和所以e^{-iβH_B} ∏_i e^{-iβX_i} ∏_i RX(2β)因为RX(θ)e^{-iθX/2}。所以在线路实现中每个量子比特上作用RX(2β)即可。初始态为什么选|⟩^n一方面的原因是它容易制备只需对所有量子比特做H门更本质的原因是|⟩是X的本征态本征值为1因而是H_B的基态。这正好和绝热演化的起点一致H_B的基态经过演化最终接近H_C的高能态。当然QAOA毕竟是变分方法初始态并非不可改但|⟩^n是经过验证的可靠默认选择。4. 基于Qiskit的完整代码实现构建、运行与验证4.1 环境选型与版本避坑我用的环境是Python 3.10配Qiskit 1.0以上的版本。这里必须提醒一个常见的坑Qiskit从1.0版本开始把模拟器从主包里拆分出去了Aer模拟器需要用pip install qiskit-aer单独安装导入语句也从from qiskit import Aer改成了from qiskit_aer import AerSimulator。如果你用的是老教程里的Aer.get_backend(qasm_simulator)在新版本里大概率直接报错。完整依赖pip install qiskit qiskit-aer scipy numpy matplotlibscipy不是量子计算必须的但它提供了minimize函数QAOA的经典优化环节直接用它不用自己写优化器。4.2 关键函数逐段拆解下面是最小可运行的完整实现。先写目标函数和线路构建部分。import numpy as np from qiskit import QuantumCircuit, QuantumRegister from qiskit_aer import AerSimulator from scipy.optimize import minimize def maxcut_value(bitstring, edges): 给定Qiskit测量返回的比特串和图的边列表计算切割数。 注意Qiskit counts字典的键最右边的字符对应第0号量子比特 所以这里把bitstring反转后再索引。 x [int(b) for b in bitstring[::-1]] cut 0 for i, j in edges: if x[i] ! x[j]: cut 1 return cut def build_qaoa_circuit(n_qubits, edges, p, gamma, beta): 构建QAOA线路。gamma和beta是长度为p的数组。 初态为|^n然后交替作用成本层和混合层p次。 qr QuantumRegister(n_qubits, q) circ QuantumCircuit(qr) # 初态对所有量子比特做H门得到|^n circ.h(qr) for layer in range(p): # 成本层每条边一个RZZ(-2*gamma) # 符号推导见正文e^{-i*gamma*(I - Z_i Z_j)/2} 的全局相位去掉后 # 等价于 RZZ(-2*gamma) for i, j in edges: circ.rzz(-2.0 * gamma[layer], i, j) # 混合层每个量子比特旋转 RX(2*beta) for i in range(n_qubits): circ.rx(2.0 * beta[layer], i) circ.measure_all() return circ def expectation_from_counts(counts, edges): 从测量统计结果中估计期望切割数 E[C]。 total sum(counts.values()) exp 0.0 for bitstring, cnt in counts.items(): exp maxcut_value(bitstring, edges) * cnt return exp / total然后是主优化流程def run_qaoa(edges, p1, shots8192, seed42): n_qubits max(max(e) for e in edges) 1 backend AerSimulator(seed_simulatorseed) def objective(params): gamma params[:p] beta params[p:] circ build_qaoa_circuit(n_qubits, edges, p, gamma, beta) counts backend.run(circ, shotsshots).result().get_counts() return -expectation_from_counts(counts, edges) # 参数初始化随机在[0, pi]之间取 rng np.random.default_rng(seed) init_params rng.uniform(0.0, np.pi, 2 * p) result minimize( objective, init_params, methodCOBYLA, options{maxiter: 500, tol: 1e-4} ) return result几个实现细节值得展开说。第一maxcut_value的位序处理。Qiskit测量结果例如1010这个字符串里最右边是q0的状态最左边是q_{n-1}的状态。所以必须反转后再比较。这个细节搞错整个期望值计算就错了但程序不会报错因为数字看起来都合理。我调试时浪费过不少时间在类似问题上。第二目标函数取负号。因为scipy.optimize.minimize是求最小值而MaxCut要最大化切割数所以在期望值前面加负号。优化器的返回值result.fun就是负的期望切割数的估计值。第三COBYLA方法是无梯度方法。它每次调用评估只跑一次量子线路对噪声相对鲁棒。在后面参数优化章节会详细展开为什么选它。4.3 跑通第一个QAOA实例并验证最优解测试图用4节点的环形图这是QAOA的hello world。edges [(0, 1), (1, 2), (2, 3), (3, 0)] n 4 res run_qaoa(edges, p1, shots8192, seed42) print(COBYLA返回的最优目标函数值负的期望切割数:, res.fun) # 同时做暴力搜索验证QAOA是否逼近全局最优 best_cut 0 best_strings [] for i in range(2 ** n): bitstring f{i:0{n}b} val maxcut_value(bitstring, edges) if val best_cut: best_cut val best_strings [bitstring] elif val best_cut: best_strings.append(bitstring) print(暴力搜索最优切割数:, best_cut) print(最优比特串:, best_strings)我自己跑出来的典型结果是COBYLA最后得到的目标函数值在-3.85到-4.0之间也就是说期望切割数非常接近4。暴力搜索确认最优切割数是4最优比特串有两个0101和1010。这两个字符串正好是彼此比特翻转的结果对应两个互补的顶点划分。这是MaxCut的天然对称性把两边对调切割数不变。然后拿出最优参数重新用更多采样次数确认最优解概率opt_params res.x gamma opt_params[:1] beta opt_params[1:] final_circ build_qaoa_circuit(n, edges, 1, gamma, beta) final_counts backend.run(final_circ, shots20000).result().get_counts() opt_prob sum(final_counts.get(s, 0) for s in best_strings) / 20000 print(最优解的总概率两个互补比特串:, opt_prob)我实测p1时这个概率通常在0.65到0.85之间。也就是说即使只跑一层QAOA测量最优解的概率也已经显著高于随机猜测的1/16。这看起来很棒但要说明的是4节点环图太小了小到经典算法几乎不需要动脑就能解。QAOA的真正价值还没有在这个规模上体现出来。做到这里重点是验证整条链路——从线路构建、期望值估计到经典优化——是正确的。接下来再进入调参和硬件问题。5. 参数优化中容易踩的坑与调参实战5.1 优化器选择背后的噪声逻辑QAOA是一个量子-经典混合算法量子计算机只负责在给定参数下跑线路、返回期望值估计经典优化器负责根据这些估计更新参数。所以优化器的选择直接影响收敛速度、稳定性和最终解质量。我用过几种常见优化器结论是当目标是模拟器且采样噪声较弱时L-BFGS-B等梯度类方法收敛最快但一旦进入真实硬件或者模拟中加入噪声模拟梯度方法的优势会迅速消失甚至失效。原因是梯度类方法对每次函数评估的精度很敏感采样噪声会让梯度估计产生偏差优化过程容易发散。COBYLA是QAOA社区里最常用的无梯度优化器之一。它不需要计算导数每次只做一次线路采样对噪声容忍度高。缺点是收敛慢尤其是在参数多的时候。如果你用p10、需要优化20个参数COBYLA可能需要上千次迭代。SPSA是另一个不错的选项它用随机方向估计梯度每次只需要两次线路采样特别适合参数多、噪声大的场景但超参数学习率、扰动幅度需要耐心调。我在模拟器上做小规模实验的建议是优先COBYLA把maxiter放开到500以上如果发现目标函数长期不降再尝试SPSA或者换一个随机种子重新初始化。5.2 初始化、局部极值与目标函数形态QAOA目标函数的形态并不简单。对4节点环图p1时如果固定βπ/4、扫描γ从0到π你会看到期望切割数的曲线有多个峰。这个多峰结构意味着经典优化器很可能陷入局部极值。最典型的陷阱是全零初始化。如果γβ0线路退化为只有初态H门期望切割数正好是所有边的一半4条边里期望切2条而且它的梯度正好为0。优化器检测到梯度为0以后会认为已经收敛——实际上你根本没有开始优化。我在早期实验里就吃过这个亏一开始以为算法坏了后来才发现是初始化的问题。所以参数初始化有两个基本经验不要全零初始化尽量在[0,π]内随机撒点并且跑多个种子比较最终目标值取最优的。利用参数的周期性。从RZZ(-2γ)和RX(2β)的结构可以看出γ和β的有效参数空间是紧致的重复多次随机初始化时不必担心跑出区间。我还建议读者在正式优化前写一个扫描脚本观察目标函数形态。对p1的小图固定一个参数、扫描另一个参数画出期望值曲线能直观理解为什么说这个目标函数充满局部极值。5.3 采样次数和期望值估计的统计边界量子线路每次运行都会坍缩到某个计算基态所以期望值C的估计是从有限次采样中得到的统计平均。采样次数shots直接决定了目标函数的信噪比。shots1024时单次评估的统计波动可能达到0.1以上这会干扰COBYLA的收敛判断shots8192时波动明显减小但每个期望值评估的时间也更长。我的习惯是分阶段设置shots优化过程中用4096到8192次采样保证优化器能看到足够平滑的目标函数又不至于太慢。最终确认最优参数时用20000次甚至更多采样重新评估得到一个可靠的解概率。还有一个容易被忽略的问题COBYLA在每次迭代中重新评估同一组参数时结果会有统计波动。这会导致优化器在最优参数附近来回抖而不是精确收敛。如果发现收敛后的目标函数在某个值附近波动不要怀疑代码这是采样噪声的正常表现。缓解办法是适当增加shots或者在最终阶段用固定的大shots重估。6. 在真实量子硬件前必须知道的工程现实6.1 线路深度、门错误率与成功概率的账QAOA在模拟器上跑得再漂亮最终还是要考虑硬件实现。先算一笔最简单的账。线路里成本层使用了RZZ门大部分真实超导硬件不直接支持RZZ需要分解为CNOT加单比特旋转RZZ(θ) CNOT(q_i, q_j) · RZ(θ, q_j) · CNOT(q_i, q_j)也就是说每隔RZZ门要消耗2个CNOT门。如果有m条边、p层QAOA那么总CNOT数大约为2mp。假设硬件的CNOT错误率是1e-2那么线路整体成功概率大约是(1-0.01)^{2mp}。我用一个真实感更强的例子m20、p10总CNOT数400成功概率(0.99)^400约等于0.018。这意味着绝大多数运行都因为门错误而得不到有效信号。所以NISQ时代的QAOA能处理的图规模是非常受限的。如果你要在大规模图m超过50上跑p10的QAOA在真实硬件上基本不可能得到有用结果。这不是某个特定平台的限制而是当前所有通用量子硬件的共性约束。应对思路有两种一是减少p在浅线路上做文章二是选择稀疏图避免完全连接图降低RZZ门总量。6.2 测量错误、拓扑约束与模拟器落差真实硬件的另一个问题是量子比特之间的连接关系。超导芯片通常只能对相邻比特做两比特门如果你的MaxCut图里有不相邻的边需要插入SWAP门把量子比特临时搬到相邻位置。一个SWAP要拆成3个CNOT代价非常高。所以我在选问题图时会先看一眼硬件拓扑尽量让图的边和硬件连接关系匹配。测量错误同样不可忽视。真实硬件的读数不是100%准确的某些比特的读数错误率可能高达百分之几。在QAOA中测量错误主要影响对期望切割数的估计从而干扰最终参数判断。你可以用校准矩阵去修正测量概率分布但实践中我发现测量错误对最优点位置的影响通常不大更需要注意的是它对这个解到底离最优多远的判断。模拟器和真机之间的落差是我在实战中感受最深的。在无噪声模拟器上QAOA的最优参数收敛得非常稳定但把同样参数放到真机上由于门错误、串扰、漂移等问题目标函数会整体变差。我通常采用的做法是先在模拟器上找到一个可靠的参数区域再基于该区域在真机上做局部微调避免从头再优化的痛苦。还有一点经验值得分享真机运行会受排队、校准状态影响同一组参数在不同时间段跑出的结果可能差异很大。所以比较参数好坏时最好集中在一个时间段内跑完避免跨越大半天引起的系统误差。这块儿的工程细节论文里很少写但实操中确实决定结果能不能看。最后我自己跑QAOA最大的体会是这个算法最难的其实不是量子部分而是经典优化和噪声工程。每次在模拟器上调通一组参数换一个图又要重新调一遍每次以为找到最优p时标题里所谓的近似优化几个字总会在硬件结果里重新提醒你一次。如果你也想用QAOA做实际问题我建议从MaxCut这种教科书问题开始先在模拟器上把它彻底跑明白再考虑真机。这个过程中的坑比任何论文里的公式都更值得踩一遍。