简介：这份资源提供模拟退火（SA）、模拟量子退火（SQA）与路径积分蒙特卡罗（PIQMC）的完整实现，面向从事组合优化、统计物理建模及量子退火算法研究的开发者与研究生。代码接口基于 Python 3 编写，底层蒙特卡洛采样由 Cython 加速，可支撑 2D Edwards-Anderson 模型、Sherrington-Kirkpatrick 模型与 Wishart Planted Ensemble 三类典型体系的实验复现，并用于生成 arXiv:2101.10154 论文中的部分结果。资源包共 184 个文件，以 166 个 txt 数据与说明文件、9 个 py 接口脚本、2 个 pyx 与 2 个 c 加速源码为主，另含 md 文档与 yml 配置，压缩后约 3.64MB，结构紧凑。该实现改编自 Hadayat Seddiqi 的 Cython 代码，修复了若干错误，为量子退火蒙特卡罗动力学加入全局移动，并精简了冗余逻辑。已有 700 余人学习，适合需要对照论文复现实验、研究退火动力学或二次开发优化求解器的读者下载研读。
1. 模拟退火与量子退火：为什么经典热涨落搞不定的优化问题，路径积分能救
一个组合优化问题，解空间动辄 2 的几百次方，穷举不可能，梯度下降又因为能量面坑坑洼洼而频繁卡在局部极小。模拟退火（SA）的思路很朴素：给系统一个虚拟温度，允许它偶尔往“更差”的方向跳，温度慢慢降下来，最终大概率停在全局最优附近。这个机制在二维 Edwards-Anderson 自旋玻璃上跑得不错，但一进入 Sherrington-Kirkpatrick 这种全连接模型，或者 Wishart Planted Ensemble 这种带随机矩阵结构的场景，经典热涨落就有点力不从心了——能量势垒又高又窄，热噪声要么跳不过去，要么跳过去之后又跳回来。
模拟量子退火（SQA）换了一条路：不靠热涨落，靠量子涨落。横场项充当那个“允许隧穿”的角色，让系统直接穿过势垒而不是翻过去。问题在于，量子退火在真实硬件上跑代价高，于是路径积分蒙特卡罗（PIQMC）成了经典计算机上模拟量子退火的标准工具——把 d 维量子系统映射成 d+1 维经典系统，多出来的那一维就是虚时间。这份用 Python 3 做接口、Cython 写采样的代码，把 SA、SQA、PIQMC 三条路线放在同一个仓库里，支持三种模型，还修了原始 Cython 实现里的一些 bug，给量子退火动力学加了全局移动。搞懂它，等于搞懂“经典热涨落和量子隧穿在优化问题里到底差在哪”。
2. 模拟退火的底层机制：Metropolis 准则为什么能跳出局部极小
2.1 从能量面到接受概率：Metropolis 到底在算什么
模拟退火的核心不是“降温”这个动作，而是 Metropolis 接受准则。给定当前状态 ( s ) 和候选状态 ( s' )，能量差 ( \Delta E = E(s') - E(s) )，接受概率是：
[ P_{\text{accept}} = \min\left(1, e^{-\Delta E / T}\right) ]
这个式子的物理含义是：温度 ( T ) 高的时候，( e^{-\Delta E / T} ) 接近 1，几乎什么候选都接受，系统在解空间里大范围游走；温度低的时候，只有 ( \Delta E < 0 ) 的候选才容易被接受，系统逐渐收敛到低能区域。关键在于，即使 ( \Delta E > 0 )，仍然有一个非零概率接受，这就是跳出局部极小的“后悔药”。
常见的实现思路是：每次迭代随机选一个自旋翻转，计算能量变化，然后和 ( \log(\text{rand}()) ) 比较。代码里用 Cython 写这部分，是因为 Metropolis 循环是热点，Python 的循环开销会让采样效率掉一个数量级。Cython 把自旋数组和能量计算都编译成 C 级别，Python 层只负责调度和参数配置。
2.2 降温 schedule：线性、指数、还是自适应
降温策略直接决定 SA 能不能收敛到全局最优。线性降温 ( T_k = T_0 - k \cdot \Delta T ) 简单，但后期温度降得太快，容易冻结在局部极小。指数降温 ( T_k = T_0 \cdot \alpha^k )（( \alpha ) 通常取 0.9 到 0.99）更常用，因为它在高温区停留时间短，低温区停留时间长，符合“先粗搜后精修”的直觉。
代码里run_SA_EA.py和run_SA_Wishart.py分别对应 Edwards-Anderson 模型和 Wishart Planted Ensemble 的 SA 运行脚本。EA 模型是稀疏图上的自旋玻璃，Wishart 模型则来自随机矩阵理论，能量面结构完全不同。同一个 SA 核心，换一个能量函数和邻域定义，就能适配不同模型——这也是为什么代码把模型定义放在models.py里，而不是硬编码在采样循环中。
2.3 全局移动：为什么局部翻转不够用
原始 SA 实现通常只做单自旋翻转，这在 EA 模型上没问题，因为每个自旋只和少数邻居耦合。但在 SK 模型这种全连接系统里，单自旋翻转的能量变化涉及所有其他自旋，计算量大，而且每次只改一个自旋，在解空间里移动得太慢。代码里给量子退火蒙特卡罗动力学加了全局移动，本质上是一次性翻转多个自旋或者整块更新，让系统在低温下也能快速跨越多个局部极小之间的势垒。
全局移动的代价是接受率会下降，因为一次改太多，( \Delta E ) 往往很大。常见的做法是引入一个“移动幅度”参数，在高温时用大移动，低温时退回小移动。代码里没有暴露这个参数给 Python 层，而是在 Cython 内部根据温度自适应调整——这也是为什么直接读qmc.c比读 Python 接口更能理解实际行为。
3. 量子退火与路径积分：虚时间维度是怎么把量子问题变成经典问题的
3.1 横场项与量子隧穿：SQA 到底在模拟什么
模拟量子退火的核心是在经典哈密顿量上再加一个横场项：
[ H = H_{\text{classical}} + \Gamma \sum_i \sigma_i^x ]
( \Gamma ) 是横场强度，( \sigma_i^x ) 是 Pauli X 矩阵。横场项的作用是让自旋在 ( z ) 方向的分量不再确定，产生量子涨落。在优化问题里，这个涨落允许系统通过隧穿效应穿过能量势垒，而不是像 SA 那样靠热涨落翻过去。当 ( \Gamma ) 逐渐减到 0，系统退回到经典哈密顿量，完成一次“量子退火”。
问题在于，量子系统不能直接在经典计算机上模拟，因为量子态是叠加态。路径积分蒙特卡罗的做法是：把量子配分函数 ( Z = \text{Tr}(e^{-\beta H}) ) 用 Trotter 公式展开，得到一系列虚时间切片，每个切片是一个经典副本，相邻切片之间通过耦合项连接。这样，d 维量子系统就变成了 d+1 维经典系统，多出来的那一维就是虚时间。
3.2 Trotter 数：虚时间切片多少才够
Trotter 数 ( P ) 是 PIQMC 里最关键的参数。( P ) 越大，虚时间离散化越细，模拟越接近真实量子系统，但计算量也线性增长。代码里run_PIQMC_EA.py会显式传入 ( P )，而qmc.c里则用这个值来分配副本数组和计算相邻副本之间的耦合强度。
常见的经验是：( P ) 至少取到 ( \beta \cdot \Gamma_{\max} ) 的几倍，否则 Trotter 误差会显著改变相变行为。在 EA 模型上，如果 ( P ) 太小，系统会表现出虚假的经典行为，量子涨落被离散化误差淹没。代码里没有自动调 ( P ) 的机制，需要用户根据 ( \beta ) 和 ( \Gamma ) 手动设置——这也是为什么 README 里强调“更多详情请下载后阅读”。
3.3 副本耦合与全局移动：PIQMC 里的并行更新
PIQMC 的每个虚时间切片是一个经典副本，副本之间通过耦合项 ( J_\perp \propto -P/\beta^2 ) 连接。在采样过程中，每个副本独立做 Metropolis 更新，但能量计算要加上与相邻副本的耦合。代码里qmc.c的采样循环同时遍历所有副本，对每个副本尝试单自旋翻转，然后根据总能量变化决定是否接受。
全局移动在这里的作用更明显：如果只做单自旋翻转，副本之间的关联会变得很强，系统在虚时间方向上扩散很慢。代码里加的全局移动允许一次性更新一个副本里的多个自旋，甚至跨副本更新，从而加速虚时间方向上的混合。这个改动是原始 Cython 代码里没有的，也是这份实现值得细看的地方。
4. 代码结构拆解：Python 接口与 Cython 采样怎么分工
4.1 文件布局：谁负责什么
仓库里的文件可以分成三层：
| 文件 | 角色 | 关键内容 |
|---|---|---|
models.py | 模型定义 | EA、SK、Wishart 的能量函数和邻域结构 |
python_interface.py | Python 接口 | 参数解析、调用 Cython 模块、结果输出 |
qmc.c | PIQMC 核心 | 路径积分采样、副本耦合、全局移动 |
sa.c | SA 核心 | Metropolis 循环、降温 schedule |
run_*.py | 运行脚本 | 针对不同模型和算法的入口 |
models.py里定义的能量函数是纯 Python 的，方便调试和验证；qmc.c和sa.c里的能量计算则是 C 级别，为了速度。两者之间通过 Cython 的cdef接口传递自旋数组和耦合矩阵。这种分工的好处是：模型逻辑改起来容易，采样核心不用动；采样核心优化时，模型逻辑不受影响。
4.2 Cython 采样循环：为什么不用纯 Python
Metropolis 循环里每次迭代要做的事情包括：随机选自旋、计算能量变化、生成随机数、比较接受概率。纯 Python 做这些，每次迭代大概几微秒，跑一百万次就是几秒。Cython 编译后，同样的循环能压到几十纳秒级别，整体快两个数量级。
代码里qmc.c的采样循环大致是这样的结构：
// 伪代码示意，基于 qmc.c 的实际逻辑 for (int step = 0; step < n_steps; step++) { for (int replica = 0; replica < P; replica++) { int spin = rand() % N; double dE = compute_energy_change(spins, replica, spin); // 加上与相邻副本的耦合项 dE += compute_coupling_change(spins, replica, spin); if (dE < 0 || rand_uniform() < exp(-dE / T)) { flip_spin(spins, replica, spin); } } // 全局移动：每隔若干步尝试一次 if (step % global_move_interval == 0) { try_global_move(spins, P, N, T); } }compute_energy_change只计算翻转一个自旋带来的局部能量变化，不需要重新算整个系统能量。compute_coupling_change处理虚时间方向上的耦合，只涉及相邻两个副本。try_global_move是后来加的，具体实现取决于模型——在 EA 模型上可能是翻转一个连通块，在 SK 模型上可能是翻转一组随机选中的自旋。
4.3 参数传递：Python 层怎么控制 Cython 行为
python_interface.py里暴露的参数包括：模型类型、系统大小、温度范围、Trotter 数、迭代步数、随机种子。这些参数通过 Cython 的def函数传给 C 层，C 层再分配数组、初始化自旋、跑采样循环。结果返回给 Python 层后，可以存成文本或者直接画图。
一个容易踩坑的地方是：C 层用的随机数生成器和 Python 的random模块不是同一个。如果不在 Python 层固定种子，每次运行结果会不一样。代码里在python_interface.py里显式调用了 C 层的种子设置函数，保证可复现性。这也是为什么直接跑run_*.py比手动调接口更省心——脚本里已经把种子、参数、输出路径都配好了。
5. 避坑与排查：SA/SQA/PIQMC 实现里最容易翻车的几个地方
5.1 现象：SA 跑出来的能量比预期高很多
原因：降温太快，或者初始温度太低。SA 的理论保证是“降温足够慢时收敛到全局最优”，但实际代码里为了速度，降温 schedule 往往偏激进。如果初始温度低于某个势垒高度，系统一开始就卡在局部极小，后面再降温也没用。
解决：先把初始温度设成比典型能量变化大一个数量级，跑一次看接受率。如果接受率低于 0.5，说明温度太低；如果高于 0.9，说明温度太高，浪费计算。代码里run_SA_EA.py的默认参数是调过的，但换模型后需要重新标定。
5.2 现象：PIQMC 结果和 SA 几乎一样，看不出量子效应
原因：Trotter 数 ( P ) 太小，或者横场 ( \Gamma ) 太小。( P ) 太小时，虚时间方向只有几个切片，量子涨落被离散化误差抹平；( \Gamma ) 太小时，横场项相对于经典哈密顿量可以忽略，系统退化成经典 SA。
解决：先把 ( \Gamma ) 设成和经典耦合强度同一量级，然后把 ( P ) 翻倍，看结果有没有变化。如果翻倍后结果不变，说明 ( P ) 够了；如果变了，继续加 ( P )。代码里run_PIQMC_EA.py的 ( P ) 默认值偏保守，跑小系统够用，大系统需要手动调大。
5.3 现象：全局移动开启后，接受率骤降
原因：全局移动一次改太多自旋，( \Delta E ) 太大，Metropolis 接受概率 ( e^{-\Delta E / T} ) 接近 0。在低温下尤其明显，因为 ( T ) 小，指数衰减快。
解决：把全局移动的幅度和温度挂钩。高温时允许大移动，低温时退回小移动。代码里qmc.c的try_global_move内部有温度判断，但阈值是硬编码的，换模型后可能需要调整。如果接受率低于 0.01，说明移动幅度太大，需要减小。
5.4 现象：Cython 编译报错，找不到qmc.h
原因：Cython 模块的编译依赖.c和.h文件，如果setup.py或者编译命令里没有正确指定 include 路径，就会报找不到头文件。
解决：确认qmc.c和qmc.h在同一个目录下，并且编译命令里加了-I.。如果用的是pyximport，需要把.c文件放在和.pyx相同的目录。代码仓库里没有提供setup.py，需要自己写一个或者用cythonize手动编译。
5.5 现象：Wishart 模型跑出来的结果和论文对不上
原因：Wishart Planted Ensemble 的生成参数（信号强度、矩阵维度、稀疏度）和论文里不一致。代码里models.py的 Wishart 实现依赖几个超参数，如果这些参数和论文里的设置不同，能量面结构会完全不一样。
解决：先确认models.py里 Wishart 类的初始化参数，然后对照论文里的实验设置逐项检查。如果论文里用了特定的随机矩阵系综，代码里也要用同样的系综。这一步没有捷径，只能逐参数核对。
6. 设计取舍：为什么用 Cython 而不是纯 C 或纯 Python
这份代码最值得琢磨的设计取舍，是“Python 接口 + Cython 采样”这个组合。纯 Python 写采样循环，开发快但跑得慢，跑一次 EA 模型的全温区扫描可能要几小时；纯 C 写全部逻辑，跑得快但改模型、调参数、画图都要重新编译，迭代成本高。Cython 卡在中间：采样循环用 C 级别速度跑，模型定义和参数配置留在 Python 层，改起来不用重新编译整个项目。
另一个取舍是“不自动调参”。代码里没有自动调 Trotter 数、自动降温、自动调全局移动幅度的机制，所有关键参数都要用户手动指定。这看起来不友好，但符合这类代码的定位：它是用来复现论文结果的，不是给生产环境用的优化器。论文里的实验设置是固定的，自动调参会引入额外变量，反而让结果难以对比。我一般会这么理解：这份代码的价值不在于“开箱即用”，而在于“每个参数都能追溯到物理含义”。
还有一个细节是全局移动的实现。原始 Cython 代码只有单自旋翻转，这份实现加了全局移动，但全局移动的策略是硬编码的，没有暴露给 Python 层。这意味着换一个模型，可能需要改qmc.c里的全局移动逻辑，而不是在 Python 层配置。这个取舍的代价是灵活性下降，好处是 C 层可以针对特定模型做优化，比如 EA 模型上翻转连通块比随机翻转一组自旋更有效。
验证这份代码是否跑对了，最直接的方法是：先用 SA 跑一个小系统，确认能量收敛曲线符合预期；再用 PIQMC 跑同样的系统，把 ( \Gamma ) 设成 0，看结果是否退化成 SA。如果两者一致，说明 PIQMC 的经典极限是对的；然后逐渐增大 ( \Gamma )，观察能量和接受率的变化。如果 ( \Gamma ) 增大后能量反而升高，说明量子涨落太强，系统被推离了低能区域，需要调整 ( P ) 或者降温 schedule。
从那以后我每次跑这类蒙特卡洛代码，都强制走一遍“经典极限验证”：先把量子项关掉，确认经典部分没问题，再逐步打开量子项。这个习惯帮我省了很多调试时间，也让我对“量子退火到底在做什么”有了更具体的感知。希望帮到你。
本文还有配套的精品资源，点击获取
