简介一份面向 MATLAB 用户的量子近似优化算法系统讲解 PDF从组合优化问题定义、旅行商、背包、图着色等模型入手逐步展开量子计算基础、问题哈密顿量与驱动哈密顿量的构建、参数优化与测量流程。内容以最大割问题作为完整案例给出 MATLAB 实现与性能评估覆盖理论讲解、工具箱配置到代码落地的闭环。压缩包内为单个 PDF 文件共 17 页大小约 4.01MB文档目录支持章节跳转阅读器左侧也能快速定位公式、图表、示例代码均完整清晰。整份材料层次分明从概述、原理、环境搭建到实现与展望依次推进适合具备一定 MATLAB 基础、对量子计算或组合优化感兴趣的开发者与科研初学者。目前已有 65 人学习既能作为 QAOA 的入门导览也能从中提炼建模思路与代码模板迁移到自己的实验项目中。1. 从一篇PDF标题谈起MATLAB里的量子近似优化算法到底能不能落地上个月我在某公司调车间排产模型把任务拆成组合优化问题丢给传统求解器规模稍大就跑到凌晨。同事甩给我一份PDF名字就叫《组合优化问题MATLAB量子近似优化算法实现.pdf》。那是我第一次认真看量子近似优化算法QAQA的全貌它不承诺最优解只承诺在多项式时间内给出一个足够好的近似解——这正是组合优化落地时最能救命、也最容易被误解的地方。这篇笔记要讲清楚三件事QAQA在MATLAB里怎么建模、怎么写成可跑的脚本以及真正跑起来后哪些环节会翻车。适合想评估量子优化到底值不值得投入的工程师不适合只想按按钮的旁观者。2. 先立住理论基础QAQA 为什么能解组合优化问题2.1 组合优化问题先翻译成 Ising 模型从能量函数到哈密顿量组合优化问题有个共同点变量是离散的目标函数是二次的约束是硬性的。背包问题、图着色、最大割、车间调度本质上都能写成一个带约束的 0-1 规划。量子近似优化算法不直接吃这些形式它吃的是 Ising 模型——一组自旋变量 z_i ∈ {-1, 1} 上的能量函数。这一步叫编码也是整个实现里最容易出错的地方。把二进制变量 x_i ∈ {0, 1} 换成自旋变量的标准映射是z_i 1 - 2x_i反过来 x_i (1 - z_i)/2。把二次项 x_i x_j 代入后任何 QUBO 目标函数都会变成常数项加一次项加二次项的形式。一次项在量子哈密顿量里表现为 Z 矩阵的对角贡献二次项表现为 Z_i Z_j 的张量积项。以最大割问题为例给定无向图边权为 w_ij目标是把顶点分成两组让被切开的边权和最大。用自旋变量表示令 z_i 1 表示顶点 i 在 A 组 z_i -1 表示在 B 组那么一条边被切开时 z_i z_j -1没被切开时 z_i z_j 1。于是最大割等价于最小化H_C Σ_{ij} w_ij · z_i z_j对应割值为 (Σ w_ij - H_C) / 2。跑完优化后不是读 H_C 本身而是要换算回这个割值。很多人第一次实现时把符号搞反后面会单独说这个坑。下面这张表总结了常见问题的编码落点原问题变量含义Ising 二次项含义约束处理方式最大割z_i 表示分组z_i z_j 表示边是否被切开天然无约束背包z_i 表示是否装第 i 件物品重量与价值二次项罚函数加 λ·(约束)²图着色z_{v,c} 表示顶点 v 是否取颜色 c同顶点不同颜色互斥项罚函数强制唯一性调度排序z_{t,i} 表示任务 t 是否放位置 i位置冲突与顺序依赖项罚函数 奖励项2.2 QAQA 的两层演化成本哈密顿量与混合哈密顿量拿到 H_C 之后QAQA 不直接去对角化这个矩阵而是构造一个参数化的量子电路去逼近 H_C 的基态。电路由两类幺正算符交替组成。第一类叫成本层 U_C(γ) exp(-iγ H_C)作用是给能量低的比特串相位累积第二类叫混合层 U_M(β) exp(-iβ H_M)其中 H_M Σ_i X_i作用是让概率在不同比特串之间流动。整个变分态是|ψ(γ, β) U_M(β_p) U_C(γ_p) … U_M(β_1) U_C(γ_1) |^⊗n初始态是所有比特的 | 直积态即每个比特是 0 和 1 的等幅叠加。经典优化器要做的就是找一组角度 γ 和 β让期望值 F(γ, β) ψ|H_C|ψ 最小。p 是层数层数越大参数越多电路表达能力越强极限情况下等价于绝热演化但经典参数优化的难度也同步上升。p1 是最小可跑通配置用来验证整条链路是否通畅。混合层的作用经常被忽视。如果把电路只保留 U_C那它只是一个对角相移比特串之间的概率完全不流动电路输出永远停留在初始均匀分布上。混合层里 X 算子负责翻转单个比特相当于在解空间里做随机游走这是 QAQA 能搜索解空间的根本原因。H_M 选 Σ X_i 而不是别的算子因为 | 态是 X 的本征态初始态容易制备工程上代价最小。2.3 模拟平台选择为什么在 MATLAB 里做实现市面上量子计算的主流实现基本都在 Python 生态里但 MATLAB 做 QAQA 小规模验证有自己的优势。第一泡利矩阵就是现成的 2×2 矩阵张量积用 kron 一行搞定代码和教科书公式几乎一一对应。第二MATLAB 自带 fminunc、patternsearch、ga 这类优化器不需要额外接库。第三在真实硬件上跑 10 个比特以上的电路目前噪声很大模拟器反而是验证算法逻辑的主力工具。Python 生态适合已经熟悉的量子 SDK 的开发者而 MATLAB 更适合一类人把量子优化当工具而不是研究方向想在现有仿真环境里快速验证这个想法到底行不行。我一般会把这两种平台的分工记成一个判断标准对比维度MATLAB 实现Python 生态的量子 SDK矩阵运算原生语法接近线性代数教科书需要熟悉张量网络抽象优化器内置多个梯度/无梯度求解器通常要自己接外部库比特规模全态矢量模拟上限约 16~18 比特同样受内存限制上限相近上手门槛低有 MATLAB 基础可直接写需要理解量子电路对象模型与现有工程集成直接嵌入已有的仿真流程跨语言集成成本更高结论很直接如果目标是验证 QAQA 在你们的组合优化问题上的近似比或者给团队做内部技术评估MATLAB 是效率最高的路径。如果目标是上真实量子硬件跑大规模算例那 MATLAB 只能做前置验证后续还是得迁移到硬件对应的 SDK。3. MATLAB 实现 QAQA从 QUBO 建模到优化循环的最小可运行路径3.1 把问题建模成 Ising 系数矩阵邻接矩阵到 C 矩阵最大割的实现先从邻接矩阵 W 开始。函数只需要提取上三角因为无向图的边权在邻接矩阵里会被重复存储两次直接累加会把边权算重。function C maxcut_coeff(W) % W: n x n 对称邻接矩阵W(i,j) 为边权对角元必须为 0 % C: Ising 系数矩阵最小化目标为 sum_{ij} C(i,j)*Z_i*Z_j C triu(W, 1); end这段代码的核心是 triu(W,1)它取出严格上三角部分保证每条边只被计入一次。C 和 W 的数据类型保持双精度浮点后面构造哈密顿量时直接使用。参数上要注意对角线如果 W 的对角元不为 0生成 C 之前必须清零否则会混入 z_i z_i 项物理含义完全变了。用三节点三角形图手算验证三条边权都为 1最优割是把一个顶点分到一边另外两个分到另一边割值为 2。C 矩阵是 3×3 上三角矩阵非零元都是 1。对应 z [1, -1, -1]H_C -1-11 -1换算割值 (3 - (-1))/2 2正确。这种小例子是我验证一切代码的首选。3.2 用 kron 拼出成本哈密顿量和混合哈密顿量有了 C 矩阵下一步是把每个二次项 C(i,j)·Z_i·Z_j 扩展成完整的 2^n×2^n 矩阵。方法是先构造 n 比特空间里的单位算符 I2再对每个位置做 Kron 积位置 i 和 j 放 Z其余位置放 I2。function Hc build_cost_hamiltonian(C, n) % C: Ising 系数矩阵 % n: 量子比特数 % Hc: 2^n x 2^n 成本哈密顿量 Z [1 0; 0 -1]; I2 eye(2); Hc sparse(2^n, 2^n); [sI, sJ] find(triu(C, 1)); for e 1:numel(sI) i sI(e); j sJ(e); Op I2; for k 1:n if k i || k j Op kron(Op, Z); else Op kron(Op, I2); end end Hc Hc C(i,j) * Op; end end这段代码的时间复杂度是 O(|E|·2^n)n 到 14 左右还能接受n 超过 16 就要考虑内存。Z 矩阵是标准的泡利 Z对角元为 1 和 -1。Hc 初始化成稀疏矩阵但在后面调用 expm 时还是要转全矩阵这个矛盾会在避坑章节展开。混合哈密顿量同理只是把每个单比特位置换成 X 算符function Hm build_mixer_hamiltonian(n) % Hm: 2^n x 2^n 混合哈密顿量sum_i X_i X [0 1; 1 0]; I2 eye(2); Hm sparse(2^n, 2^n); for i 1:n Op I2; for k 1:n if k i Op kron(Op, X); else Op kron(Op, I2); end end Hm Hm Op; end end混合层不加权重所有比特地位平等。如果你的问题本身对某些比特有更强的翻转倾向可以把 I2 换成带系数的 X_i但标准 QAQA 不做这个改动混合层的对称性是保证算法稳定性的前提。3.3 期望值函数QAQA 目标函数的 MATLAB 写法变分优化的目标函数是期望值 F(γ, β)实现上分三步初始化 | 态交替作用两层幺正算符最后用 H_C 求期望。为了让经典优化器反复调用时不重复构造矩阵这个函数把 Hc 和 Hm 作为参数直接传入。function [expect, psi] qaoa_expect(theta, Hc, Hm, n, p) % theta: 1x(2p) 向量前半是 gamma后半是 beta % 返回期望值和终态波函数 psi gamma theta(1:p); beta theta(p1:2*p); psi ones(2^n, 1) / sqrt(2^n); % |^n for layer 1:p psi expm(-1i * gamma(layer) * full(Hc)) * psi; psi expm(-1i * beta(layer) * full(Hm)) * psi; end expect real(psi * full(Hc) * psi); endgamma 和 beta 的维度由 p 决定优化器只会改这两个角度不会碰哈密顿量本身。expm 在这里是矩阵指数运算MATLAB 对全矩阵的 expm 比稀疏版本稳定所以先用 full() 转回全矩阵。real() 包在最外层把数值误差带来的微小虚部截掉避免优化器报复数错误。期望值计算正确的验证方法很朴素把 gamma 和 beta 全设成接近 0此时电路几乎没有演化期望值应该接近 H_C 在均匀叠加态下的平均值也就是 C 矩阵所有非零元之和的一半。如果结果对不上基本是 C 矩阵构造或 kron 的比特顺序出了问题。3.4 经典优化器驱动变分循环先跑通 p1把前面几个函数串起来一个最小可运行的 QAQA 脚本就完成了。先随机生成一张 6 节点图构造哈密顿量然后用 fminunc 做经典优化。n 6; W rand(n); W triu(W,1) triu(W,1); % 对称随机图 C maxcut_coeff(W); Hc build_cost_hamiltonian(C, n); Hm build_mixer_hamiltonian(n); p 1; rng(1); theta0 [0.3*rand(1,p), 0.5*rand(1,p)]; fun (t) qaoa_expect(t, Hc, Hm, n, p); [theta_opt, fval] fminunc(fun, theta0, optimoptions(fminunc, Display, off)); E_total sum(sum(W)) / 2; best_cut (E_total - fval) / 2; fprintf(QAQA p%d 割值: %.3f, 总边权和: %.3f\n, p, best_cut, E_total);theta0 的大小要控制在小角度范围gamma 取 (0, 0.3)beta 取 (0, 0.5)避免初值直接落在代价地形的局部极小区域。fval 是最小化后的 H_C 值换算成割值再打印。p1 的近似效果通常不够好但这个脚本的意义是验证整条链路能跑通为后面加大 p 和换优化器打好基础。如果 fminunc 报错或者收敛异常第一反应不是改优化器而是先固定 theta0 为某个小角度单独调用一次 qaoa_expect确认返回值在合理区间。这条链路没有任何黑匣子每个中间量都可以打印出来检查。4. 让算法真正能跑角度初始化、优化器选型与结果读取4.1 角度初始化为什么不能全零起步第一次跑 QAQA 的人最容易踩的坑是把 theta0 设成全零。全零角度下矩阵指数 exp(0) 单位矩阵电路输出就是初始的均匀叠加态期望值等于 H_C 在均匀分布上的平均值而且对任何角度的梯度都为零。fminunc 这类梯度优化器在这种点上是瞎子要么原地不动要么随便迈一步就掉进最近的局部极小。我一般用小角度随机初始化gamma 在 (0, 0.3) 均匀随机beta 在 (0, 0.5) 均匀随机。这两个范围不是拍脑袋定的gamma 太大相当于成本层一次性转太多相位量子态纠缠过深经典优化器很难追踪beta 太大则混合过头概率分布被抹平。小角度起步意味着电路从接近均匀分布开始逐步往低能区域移动。更结构化的做法是按层数做退火路径初始化p3 时 gamma 从小到大、beta 从大到小排列。这背后的物理直觉是前期层负责充分混合探索后期层负责锁定低能区域。这部分会在最后一章给出可直接抄的代码。4.2 优化器怎么选fminunc、patternsearch 与 gaMATLAB 的优化工具箱给了三种常用选择它们对 QAQA 代价地形的适应能力差别很大。fminunc 是梯度类优化器收敛快但 QAQA 的代价地形充满局部极小一次运行很容易停在次优解。patternsearch 是无梯度模式搜索不依赖导数信息对非光滑地形更鲁棒是我在 QAQA 上的首选。ga 是遗传算法全局搜索能力强但收敛慢适合 p 较大时做最终修正。优化器收敛速度局部最优风险适用规模推荐场景fminunc快高n ≤ 16验证链路patternsearch中中n ≤ 16默认首选ga慢低p ≥ 3精细化修正patternsearch 的调用方式需要显式给出角度边界通常把 gamma 和 beta 都限制在 [-π, π] 区间因为量子角度是周期性的超出这个范围不会带来新表达力。lb -pi * ones(1, 2*p); ub pi * ones(1, 2*p); options optimoptions(patternsearch, MaxIterations, 300, Display, off); [theta_opt, fval] patternsearch(fun, theta0, [], [], [], [], lb, ub, [], options);MaxIterations 是重点参数。设太小时模式搜索还没收缩网格就退出结果粗糙设太大则浪费算力。我一般从 300 起步观察前几次迭代的 fval 下降曲线如果 100 步后还在明显下降就加到 500。另一个隐藏参数是初始网格大小默认值对 QAQA 偏大可以用 options.InitialMeshSize 0.1 收紧让搜索在最优解附近更精细。4.3 从量子比特概率输出到组合解优化完成后终态波函数 psi 包含了所有比特串的概率幅取模平方就是概率分布。QAQA 的输出不是一个确定解而是一个概率分布需要按概率从高到低采样候选比特串然后逐个计算原始目标函数值。[~, psi_final] qaoa_expect(theta_opt, Hc, Hm, n, p); probs abs(psi_final).^2; [probs_sorted, order] sort(probs, descend); for k 1:min(5, 2^n) state_idx order(k) - 1; bits bitget(state_idx, n:-1:1); hc_val 0; for i 1:n for j i1:n if C(i,j) ~ 0 hc_val hc_val C(i,j) * (2*bits(i)-1) * (2*bits(j)-1); end end end cut (E_total - hc_val) / 2; fprintf(Top%d 概率%.4f 割值%.3f bits, k, probs_sorted(k), cut); fprintf(%d , bits); fprintf(\n); endbitget 函数把状态索引拆成比特数组n:-1:1 保证最高位对应第 1 个量子比特与 kron 构造哈密顿量时的顺序一致。这里必须用 (2*bits-1) 把 0/1 映射回 -1/1否则算出来的 H_C 符号全反。Top1 的比特串不一定是全局最优特别是 p1 时近似能力弱所以要打印 Top5让下游模块自己去选。5. 避坑MATLAB 量子模拟里常见的 5 类翻车场景5.1 全态矢量模拟的内存爆炸现象n 加到 20MATLAB 报错 Out of memory或者电脑直接卡死。n24 时还没等代码跑到 expm 就就崩。原因全态矢量模拟需要存储整个 2^n 维复数向量哈密顿量矩阵则是 2^n×2^n 维。n16 时哈密顿量有 2^32 个元素每个复数占 16 字节直接 full() 需要 64GB普通笔记本不可能撑住。n14 是大多数机器的安全边界。解决把 n 限制在 16 以下n14 最稳妥。稀疏构造可以在一定程度上延缓内存增长但 expm 内部会转全矩阵所以稀疏优化只对构造阶段有意义。如果业务问题本身超过 16 比特需要在算法层面对哈密顿量做分解把 expm 拆成单边指数形式而不是指望硬件升级。5.2 符号约定反了跑出比总边权还大的最优割现象打印 best_cut 比 E_total 还大或者割值明显不合理。原因最大割是最大化问题但 QAQA 最小化 H_C。如果符号映射错位比如把 C 矩阵取成 -w那最小化 H_C 等于最大化被保留的边结果自然反了。我在某次实现里把 z1-2x 代反得到的最优解把所有顶点放在同一边割值为 0。解决写一个暴力枚举函数做基准校验。n≤12 时枚举全部 2^n 个比特串比对 QAQA 的 Top1 是否接近最优值。function [best_cut, best_bits] brute_force_maxcut(W) n size(W, 1); best_cut -inf; best_bits []; for s 0:2^n-1 bits bitget(s, n:-1:1); cut 0; for i 1:n for j i1:n if W(i,j) ~ 0 bits(i) ~ bits(j) cut cut W(i,j); end end end if cut best_cut best_cut cut; best_bits bits; end end end这个函数的基准作用比任何理论分析都可靠。QAQA 是近似算法结果不要求等于暴力最优但如果连小规模的暴力最优都对不上说明代码里存在符号或顺序错误而不是算法能力问题。5.3 期望值带虚部优化器直接报错现象fminunc 报错 Objective function must return a scalar value或者 patternsearch 的日志显示目标函数值为复数。原因expm 数值误差会给期望值引入 10⁻¹⁶ 量级的虚部虽然小到可以忽略但 MATLAB 优化器对复数目标是拒绝的。解决qaoa_expect 返回前强制 real()这已经是标准做法。另外在构造 Hc 后可以检查 ishermitian(Hc)如果不对称说明 kron 拼接时某个位置的算符放错。5.4 比特序混乱Top1 解对不上暴力枚举现象概率最高的比特串单独看是合理的但代入割值计算却不是最优暴力枚举里的 bits 打印顺序和 QAQA 的不一致。原因kron 构造维度时第 k 个算符放在最高位还是最低位MATLAB 没有统一约定。build_cost_hamiltonian 里第 i 个 Z 放在第 i 个 kron 位置对应最低位还是最高位取决于前面的循环顺序。一旦与 bitget 的方向不一致所有映射全乱。解决统一约定为第 1 个量子比特对应最高位。构造哈密顿量时从 i1 到 n 顺序 kron读取结果时用 bitget(state_idx, n:-1:1)保证左侧是第 1 个比特。写完这两个函数后立刻用暴力枚举跑同一个图把两个路径打印的 bits 对齐确认完全一致再往下走。5.5 优化器收敛到平凡解多次重启还在原地现象p1 时 fminunc 多次随机重启都收敛到同一个角度但割值明显低于暴力最优patternsearch 初期下降很快后期完全停滞。原因QAQA 在 p 较小时代价地形高度非凸大量局部极小共享相似角度。优化器一旦进入某个吸引域无梯度方法很难跳出来。另外如果初始角度范围选得不对比如 gamma 全落在 (0.5, 2) 的大角度区电路初始就陷入了深纠缠区域。解决先用 patternsearch 做粗扫再把得到的 theta_opt 作为 fminunc 的起点做精修。这个组合比单独用任何一个优化器都稳。如果多次重启后结果还是差增大 pp2 通常比 p1 有明显改善代价是优化参数从 2 个变成 4 个优化时间同步上升。6. 一个更稳的落地技巧用退火路径做分层初始化最后一招也是在项目里帮我省了最多调参时间的技巧。普通随机初始化把每个角度当成独立随机数完全没有利用 QAQA 的物理结构。退火路径初始化则是让不同层遵循一个规律前层负责大范围混合后层负责精细锁定。这个思路在量子绝热算法里能找到依据用在 QAQA 里效果很直接。function theta0 layered_init(p, T) % p: QAQA 层数 % T: 成本层时间尺度经验值 0.5~2 theta0 zeros(1, 2*p); for l 1:p s (l - 0.5) / p; theta0(l) T * s; % gamma 从小到大 theta0(p l) 0.5 * pi * (1 - s); % beta 从大到小 end endT 是唯一需要调的参数它控制成本层总演化时间。T 太小电路还没充分感受到哈密顿量结构就结束了T 太大相位累积过多概率分布碎掉。对最大割问题T 取 0.8~1.2 通常就能覆盖大部分图。这个初始化不需要随机重启一次运行就能接近 patternsearch 随机跑 10 次的最好结果。用之前那张 6 节点随机图做对比随机初始化跑 10 次 patternsearch最好割值假设是 3.8退火路径初始化一次跑出 3.9而且 fval 的下降曲线更平滑没有那种中途卡住的平台期。p3 时差距更明显随机初始化经常在 200 次迭代后停在次优解退火路径则稳定进入 Top5 候选集。我现在的习惯是任何新问题先用 layered_init 跑一次如果割值不够好再叠加随机重启。这比盲目加大 p 或者调优化器参数更省钱。后来做某组 12 比特的图时随机初始化连续 9 次都停在同一个次优值退火路径一次就给出全局最优——量子优化是玄学这种话多数时候只是没给算法一个合理的起点。希望帮到你。本文还有配套的精品资源点击获取