1. 项目概述当统计物理遇见机器学习在统计物理的研究工具箱里对偶性一直是一个既强大又神秘的存在。简单来说它就像是为一个复杂的物理系统找到了另一套完全不同的“语言”来描述它但这套新语言描述的物理本质却完全一样。最经典的例子莫过于二维伊辛模型的Kramers-Wannier对偶一个在高温下无序、低温下有序的自旋系统其高温展开竟然与另一个自旋系统的低温展开在数学上等价并且两者通过一个优美的公式β̃ -1/2 log(tanh β)联系起来。传统上发现这样的对偶关系依赖于物理学家深刻的物理直觉和精巧的数学变换比如高温展开、傅里叶变换或规范理论。然而随着模型变得复杂——比如引入次近邻耦合、多体相互作用或在非平凡格点上——寻找解析对偶变换的难度呈指数级增长。这时一个自然而然的想法是能否让机器来帮我们“发现”这些隐藏的对称性这正是我们这项工作的核心将对偶性发现形式化为一个机器学习优化问题。我们不再依赖先验的解析猜测而是构建一个框架让算法通过分析系统的蒙特卡洛模拟数据自动寻找可能的对偶映射和耦合参数。这项工作的价值在于它为我们探索未知的强关联系统、寻找新的普适类甚至理解拓扑序与对称性保护相提供了一个全新的、数据驱动的视角。对于从事统计物理、凝聚态理论或格点场论研究的同行以及希望将机器学习应用于基础科学发现的探索者本文将详细拆解我们如何将这一想法落地从核心思路到代码实现并分享一路走来的经验与教训。2. 核心思路与算法框架设计2.1 问题形式化什么是对偶性发现在传统机器学习中一个常见任务是“逆伊辛问题”给定一组自旋构型样本{σ_i}推断出产生这些样本的哈密顿量H(σ) -Σ K_ij σ_i σ_j中的耦合常数K_ij。这本质上是一个参数拟合问题。对偶性发现则是一个更宏大的“元问题”。我们不仅不知道耦合参数甚至连描述系统的基本变量、它们所居住的“空间”格点以及可观测量之间的映射关系都是未知的。具体来说我们需要同时学习对偶格点上的耦合参数β̃在对偶描述中相互作用强度是多少对偶格点的几何结构L̃自旋是放在原来的格点上还是其对偶格点如正方格子的对偶仍是正方格子三角格子的对偶是蜂窝格子可观测量映射G(σ̃)如何将原格点上的物理量如自旋关联函数用对偶格点上的变量σ̃表达出来我们的目标是找到一个三元组(β̃, L̃, G)使得用这个三元组定义的对偶理论其所有可观测量的统计分布或至少其低阶矩与原理论完全一致。这就构成了一个损失函数的基础衡量两个理论预测的关联函数之间的差异。2.2 两大技术路径机器学习方法与拓扑线方法在实际操作中我们探索了两种互补的技术路径它们各有优劣适用于不同场景。2.2.1 机器学习方法端到端的黑箱优化这种方法最为直接和通用。我们构建一个可微分的模型其输入是原系统的自旋样本输出是对偶系统的预测关联函数。模型内部包含可学习的参数分别对应β̃、L̃的某种参数化表示例如一个选择哪些链接有效的注意力权重以及映射G例如一个小型神经网络。损失函数通常定义为原系统与对偶系统预测的若干低阶关联函数如一阶、二阶矩之间的均方误差。通过梯度下降如Adam优化器同时优化所有参数当损失收敛到零附近时我们就找到了一个候选对偶。实操心得特征工程的重要性在最初的尝试中我们直接将原始自旋构型输入一个全连接网络效果很差。后来我们发现必须将物理先验知识编码进去。例如我们不是输入所有自旋而是输入一个局部窗口内所有可能链接上的自旋乘积即σ_i σ_j作为特征。这相当于告诉网络“对偶变换很可能只与局部自旋关联有关”。这大大提升了训练的稳定性和可解释性。链接的选择如图2中的7种可能成为了L̃学习的一部分。2.2.2 拓扑线方法利用对称性的白盒优化第二种方法限制更多但更高效、更可解释。它专门用于寻找Kramers-Wannier型的对偶即由Z2自旋翻转对称性所生成的拓扑缺陷线拓扑线所实现的对偶。这类对偶的特点是对偶格点就是原格点的几何对偶如正方格子的对偶仍是正方格子并且可观测量映射G具有确定的形式。具体来说对于原系统的一个Z2对称性算子我们可以定义其对应的拓扑缺陷线。穿过这条线的自旋关联函数会以确定的方式映射到对偶系统的自旋算子上。对于伊辛模型这个映射就是著名的Kramers-Wannier对偶映射对偶自旋σ̃_p位于原格点的 plaquette方格中心其期望值由穿过该 plaquette 边界的原自旋弦算子的乘积给出。在这个强约束下问题大大简化L̃和G已知只剩下对偶耦合参数β̃需要学习。这就将问题规约到了一个增强版的逆伊辛问题。3. 拓扑线方法的数值实现详解鉴于拓扑线方法在效率和物理清晰度上的优势我们将其作为本文深入剖析的重点。下面我们以包含次近邻相互作用的广义伊辛模型为例一步步拆解算法流程。3.1 模型与目标我们考虑的原模型哈密顿量为H̃[σ̃] -β̃ Σ_ĩj̃ σ̃_ĩ σ̃_j̃ - κ̃ Σ_[ĩj̃] σ̃_ĩ σ̃_j̃其中ĩj̃表示最近邻相互作用[ĩj̃]表示对角方向的次近邻相互作用见图6。β̃和κ̃是耦合常数。我们的目标是假设对偶模型也是一个定义在正方格点上的伊辛模型可能也包含次近邻耦合其哈密顿量为H[σ] -β Σ_ij σ_i σ_j - κ Σ_[ij] σ_i σ_j。我们要从原模型的蒙特卡洛样本出发找出最优的(β, κ)使得两个模型在可观测量上匹配。3.2 核心算法基于牛顿法的耦合学习算法的核心思想源于逆伊辛问题中的一个关键等式。对于一个广义伊辛模型H -Σ_{i,δ} K_{i,iδ} σ_i σ_{iδ}利用自旋的离散性σ_i ±1可以推导出一个自洽方程类似于系统的“运动方程”σ_j σ_{jδ} tanh( Σ_{δ} K_{j,jδ} σ_{jδ} ) * σ_{jδ} 这个等式的含义是两点关联函数可以用一个包含所有与σ_j耦合的自旋的非线性函数双曲正切的期望值来表示。如果我们的耦合参数K_{δ}是正确的那么这个等式应该严格成立。算法步骤数据准备使用蒙特卡洛方法如Metropolis-Hastings算法对原哈密顿量H̃进行采样得到大量自旋构型{σ̃_i}。计算对偶关联函数利用已知的拓扑线映射关系即Kramers-Wannier对偶公式从{σ̃_i}样本中计算出对偶自旋{σ_i}的所有偶数阶关联函数。注意在有限系统中奇数阶关联函数恒为零。关键难点我们无法直接得到对偶自旋σ_i的样本只能得到它们的关联函数。而上述自洽方程III.18的右边需要计算像tanh(...) σ_{jδ}这样的非线性关联函数这无法直接从关联函数的集合中简单获得。重构边际概率分布为了解决上述难点我们采用了一个“迂回”但精确的方法。由于耦合K_δ是局域的假设其作用范围为R那么要计算涉及中心自旋σ_j的非线性关联函数我们只需要知道以j为中心、边长为N 2R的窗口内所有自旋的联合概率分布p({σ_a})其中a标记窗口内的N^2个自旋。这个分布p({σ_a})是一个2^{N^2}维的向量因为每个自旋有±1两种状态。窗口内所有可能的自旋乘积的期望值即所有偶数阶关联函数构成了一个线性方程组其系数矩阵由自旋乘积的取值±1决定右边是步骤2中计算出的关联函数值。通过求解这个2^{N^2} × 2^{N^2}的线性方程组III.21我们就可以反解出完整的边际概率分布p({σ_a})。这一步在N较小如3或4时是可行的因为矩阵尺寸512×512或65536×65536在现代计算机上尚可处理。迭代优化耦合参数 a.计算梯度有了p({σ_a})任何非线性函数f({σ_a})的期望值都可以通过直接求和计算f Σ_{σ_a} p({σ_a}) f({σ_a})。利用此公式我们可以计算自洽方程III.18两边的值并定义残差。 b.牛顿法更新我们计算残差关于耦合参数K_δ的雅可比矩阵D_{δ,δ}即梯度公式III.19。然后求解线性方程残差 Σ_{δ} D_{δ,δ} ΔK_{δ}得到耦合参数的更新量ΔK_δ。 c.迭代以学习率ε更新参数K_δ ← K_δ ε ΔK_δ。重复步骤4a-4c直至残差收敛到零附近。3.3 针对次近邻模型的调整与挑战当我们处理包含次近邻耦合κ̃的模型时情况变得更有趣。传统的Kramers-Wannier对偶解析推导会变得异常复杂见附录A.2因为次近邻相互作用在“规范固定”时存在多种等价的连接方式导致对偶模型中会出现非局域或多体相互作用。在我们的数值方法中我们采取了一种实用主义和近似的视角假设我们限制对偶哈密顿量H[σ]的形式只允许它包含最近邻 (β) 和次近邻 (κ) 耦合。这是一个模型选择。目标在这个受限的模型空间内寻找最优的(β, κ)使得其对关联函数的预测与原模型通过拓扑线映射计算出的对偶关联函数最匹配。结果我们发现即使这个受限的对偶模型在严格意义上并非原模型的精确对偶解析推导表明精确对偶可能包含更复杂的相互作用但在一个相当宽的参数范围内特别是当κ̃较小时它提供了一个惊人的优秀近似。图8显示两个模型的关联函数在多个距离上几乎重合。深度解析为什么算法在无序相失效在图7中一个非常显著的现象是当原系统处于无序相β̃ β̃_c时我们的算法无法可靠地恢复出有意义的耦合参数。八个独立链接的耦合值不再聚类为预期的两组β和κ而是变得杂乱无章。 这背后有深刻的物理和算法原因。在无序相系统的关联长度很短自旋涨落很大。我们用于重构边际概率分布p({σ_a})的窗口大小N是固定的。如果关联长度远小于N那么窗口边缘的自旋与中心自旋几乎独立导致我们试图从关联函数反推联合概率分布时问题变得病态ill-posed。矩阵(III.21)接近奇异微小噪声会导致解的巨大波动。此外在无序相系统更接近高温展开其对偶应处于有序相β很大。学习大耦合参数本身就是逆伊辛问题中的难点因为此时概率分布非常尖锐梯度信号很弱。这提醒我们算法的成功应用依赖于物理系统的状态需要谨慎选择参数区域。4. 实操过程从代码到结果4.1 环境准备与依赖安装我们的实现主要基于 Python 的科学计算栈。以下是核心依赖# 核心计算与自动微分 pip install numpy torch # 蒙特卡洛模拟与数据分析 pip install scipy pandas # 可视化 pip install matplotlib seaborn注意为了复现论文中的精确结果需要确保随机数种子固定并且使用双精度浮点数进行计算特别是在求解线性方程组步骤3时以避免数值误差积累。4.2 蒙特卡洛采样模块实现我们采用标准的 Metropolis 算法生成自旋构型。关键在于高效计算能量差。import numpy as np import torch def ising_energy(spins, beta, kappa, latticesquare): 计算包含次近邻耦合的伊辛模型能量。 spins: 二维 numpy 数组或 torch 张量值为 ±1。 beta: 最近邻耦合。 kappa: 次近邻对角耦合。 lattice: 格点类型默认为正方。 # 最近邻相互作用 energy_nn -beta * (np.roll(spins, 1, axis0) np.roll(spins, -1, axis0) np.roll(spins, 1, axis1) np.roll(spins, -1, axis1)) * spins # 次近邻相互作用两个对角线方向 energy_diag -kappa * (np.roll(np.roll(spins, 1, axis0), 1, axis1) np.roll(np.roll(spins, 1, axis0), -1, axis1) np.roll(np.roll(spins, -1, axis0), 1, axis1) np.roll(np.roll(spins, -1, axis0), -1, axis1)) * spins return np.sum(energy_nn) np.sum(energy_diag) def metropolis_sweep(spins, beta, kappa, n_sweeps1000): 执行多次 Metropolis 扫描。 L spins.shape[0] for _ in range(n_sweeps * L * L): i, j np.random.randint(0, L, size2) delta_E 2 * spins[i, j] * ( beta * (spins[(i1)%L, j] spins[(i-1)%L, j] spins[i, (j1)%L] spins[i, (j-1)%L]) kappa * (spins[(i1)%L, (j1)%L] spins[(i1)%L, (j-1)%L] spins[(i-1)%L, (j1)%L] spins[(i-1)%L, (j-1)%L]) ) if delta_E 0 or np.random.rand() np.exp(-delta_E): spins[i, j] * -1 return spins4.3 核心算法边际分布重构与牛顿法优化这是整个项目的核心我们将其分解为几个函数。def compute_dual_correlations_from_original(spins_original, window_size3): 步骤2从原自旋样本通过Kramers-Wannier映射计算对偶自旋在窗口内的所有关联函数。 输入: spins_original - 原格点自旋样本 (L, L) 输出: corrs_dict - 字典键为关联函数索引元组值为估计的期望值。 例如corrs_dict[((0,0), (1,0))] 表示 σ_(0,0) σ_(1,0) L spins_original.shape[0] # 对偶自旋位于原格点plaquette中心索引可视为原格点链接的右上角顶点 # 简化实现将对偶自旋定义为周围四个原自旋的乘积 dual_spins np.ones((L, L)) for i in range(L): for j in range(L): # 注意边界条件这里使用周期性边界 dual_spins[i, j] (spins_original[i, j] * spins_original[(i1)%L, j] * spins_original[i, (j1)%L] * spins_original[(i1)%L, (j1)%L]) # 现在计算 window_size x window_size 窗口内所有对偶自旋的关联函数 # 这里需要计算所有可能的自旋乘积的期望值即所有 2^(N^2) 种组合。 # 为了效率我们通常不会显式计算所有组合而是通过后续的矩阵方程隐式处理。 # 更高效的做法是直接构建所有可能的自旋乘积模式并计算其在样本上的平均值。 # 这部分代码较长核心是遍历所有可能的 α_a ∈ {0,1}计算对应的关联函数。 # ... return corrs_dict def solve_marginal_distribution(corrs_dict, window_size3): 步骤3从关联函数反解边际概率分布 p({σ_a})。 输入: corrs_dict - 上一步计算的所有关联函数。 window_size - 窗口边长 N。 输出: p_vec - 长度为 2^(N^2) 的向量表示所有可能自旋构型的概率。 N2 window_size * window_size total_states 1 N2 # 2^(N^2) # 构建矩阵 M 和向量 b M np.zeros((total_states, total_states)) b np.zeros(total_states) # 为每个自旋构型索引 s (0 到 2^(N^2)-1) 和每个关联函数模式 α (也用一个整数表示) # 填充矩阵 M[s, α] Π_a (σ_a)^(α_a)其中 σ_a 根据构型s确定。 # 向量 b[α] corrs_dict 中对应的关联函数值。 # 这是一个大规模的但稀疏的因为很多关联函数为0线性方程组。 # 实践中由于N3时矩阵大小为512x512可以直接用np.linalg.solve求解。 # ... p_vec np.linalg.solve(M, b) # 确保概率和为1且所有元素非负可进行小幅裁剪和归一化 p_vec np.clip(p_vec, 1e-10, None) p_vec / p_vec.sum() return p_vec def compute_gradient_and_update(K_current, p_marginal, learning_rate0.01): 步骤4计算梯度并用牛顿法更新耦合参数 K。 K_current: 当前耦合参数向量例如 [β, κ, ...] 对应不同方向的链接。 p_marginal: 边际概率分布向量。 # 计算所有需要的期望值左边 σ_j σ_{jδ} 和右边 tanh(...) σ_{jδ} # 以及雅可比矩阵 D_{δ,δ} ∂σ_j σ_{jδ} / ∂K_{δ} # 这需要遍历边际分布中的所有构型计算对应的自旋值和双曲正切值。 left_expectations [] right_expectations [] jacobian [] # ... (具体计算涉及对 p_marginal 的加权求和) # 构建线性系统残差 J * ΔK residual np.array(left_expectations) - np.array(right_expectations) J np.array(jacobian) # 求解 ΔK。注意J可能奇异使用伪逆或添加正则项。 delta_K np.linalg.lstsq(J, residual, rcond1e-6)[0] # 更新参数 K_new K_current learning_rate * delta_K return K_new, np.linalg.norm(residual)4.4 主训练循环与结果可视化将上述模块串联起来形成完整的工作流。def main_training_loop(beta_original, kappa_original, lattice_size12, window_size3, epochs100): 主训练循环。 # 1. 生成原模型数据 print(fGenerating MC samples for (β̃, κ̃) ({beta_original}, {kappa_original})) spins_init np.random.choice([-1, 1], size(lattice_size, lattice_size)) spins_mc metropolis_sweep(spins_init.copy(), beta_original, kappa_original, n_sweeps5000) # 通常需要多个独立样本这里简化为使用一个平衡后的构型进行多次测量。 # 2. 初始化对偶耦合参数猜测 K np.random.randn(8) * 0.1 # 假设8个独立链接参数 losses [] for epoch in range(epochs): # 3. 计算对偶关联函数 dual_corrs compute_dual_correlations_from_original(spins_mc, window_size) # 4. 求解边际分布 p_marginal solve_marginal_distribution(dual_corrs, window_size) # 5. 牛顿法更新 K, loss compute_gradient_and_update(K, p_marginal, learning_rate0.05) losses.append(loss) if epoch % 10 0: print(fEpoch {epoch}: Loss {loss:.6f}, K {K}) if loss 1e-6: print(Converged!) break # 结果分析将学习到的8个参数聚类为β和κ # 对于正方格子的最近邻和次近邻模型8个链接可分为两组 # 4个最近邻链接应趋于相同的值 β4个次近邻链接应趋于相同的值 κ。 beta_learned np.mean(K[[0, 2, 4, 6]]) # 假设索引0,2,4,6对应最近邻 kappa_learned np.mean(K[[1, 3, 5, 7]]) # 假设索引1,3,5,7对应次近邻 print(fLearned parameters: β {beta_learned:.4f}, κ {kappa_learned:.4f}) return beta_learned, kappa_learned, losses运行此代码对于(β̃, κ̃) (0.45, 0.05)我们可能得到类似β ≈ 0.261, κ ≈ 0.019的结果。然后可以计算并绘制两个模型的关联函数进行对比如图8所示验证对偶近似的质量。5. 常见问题、调试技巧与经验总结在实际复现和扩展这项工作的过程中我们遇到了诸多挑战。以下是总结出的关键问题和解决方案。5.1 数值不稳定与病态问题问题表现在求解边际分布p({σ_a})时np.linalg.solve程序报错或解中出现巨大的负值或大于1的值。牛顿法更新步长ΔK异常大导致参数发散。在无序相 (β̃ β̃_c) 结果完全不可信。根源与解决方案矩阵条件数过大方程III.21的矩阵M可能接近奇异。这是因为当窗口内自旋关联很弱时许多关联函数近似线性相关。解决方案使用正则化技术。将np.linalg.solve(M, b)替换为np.linalg.lstsq(M, b, rcond1e-8)或添加 Tikhonov 正则项(M^T M λ I) Δp M^T b其中λ是一个小的正数如1e-6。采样噪声蒙特卡洛样本有限导致关联函数b的估计有误差噪声在反问题中被放大。解决方案增加蒙特卡洛采样量并使用更先进的采样技术如 Wolff 集群算法来减少临界点附近的临界慢化。对多个独立样本计算关联函数并取平均。无序相的根本困难如前所述在无序相问题本身是病态的。此时不应期望算法能给出稳定结果。一个实用的检查点是计算出的边际分布p({σ_a})是否接近均匀分布所有构型概率接近1/2^{N^2}如果是说明系统缺乏有效信息算法失效是预期的。5.2 收敛性与超参数选择问题表现损失函数震荡不下降。学习到的耦合参数β,κ不收敛到稳定值或者收敛到明显错误的解如所有耦合都趋于零。调试技巧学习率ε这是最关键的超参数。起始学习率建议设在0.01到0.1之间。如果损失震荡降低学习率如果下降过慢适当增加。可以实现学习率衰减例如每50轮乘以0.9。初始化耦合参数K的初始值不宜设为零或过大。零初始化可能导致梯度也为零。建议从N(0, 0.1)高斯分布中随机初始化。损失监控除了总损失还应监控自洽方程III.18中每个链接δ对应的残差。这有助于识别是某个特定方向的耦合学习有问题还是整体问题。梯度裁剪在牛顿法求解ΔK后如果||ΔK||过大如 1.0可以将其按比例缩放防止单步更新过大破坏稳定性。5.3 性能优化与扩展性瓶颈瓶颈分析 整个算法最耗时的部分是边际分布求解涉及O(2^{3N^2})复杂度的矩阵构建与求解N3时为512^3量级。N4时矩阵尺寸为65536 x 65536内存和计算需求激增。期望值计算每次迭代都需要遍历所有2^{N^2}个构型来计算非线性期望值和雅可比矩阵计算量为O(2^{N^2} * num_links^2)。优化建议利用对称性对于具有平移、旋转对称性的系统许多链接的耦合K_δ是等价的。可以强制约束这些链接参数相等从而减少待优化参数数量同时雅可比矩阵的维度也大大降低。稀疏性与快速变换矩阵M和计算期望值的求和具有高度结构化的模式。探索是否能用快速沃尔什-哈达玛变换来加速计算因为自旋取值为±1所有计算本质上是在布尔超立方体上的求和。替代方案对于更大的N直接矩阵求逆不可行。可以考虑使用变分推断或训练一个小型神经网络来直接从可测量的关联函数估计非线性期望值tanh(...) σ从而绕过显式求解p({σ_a})这一步。这正是原文未来工作展望中提到的方向。5.4 物理结果解读与验证如何判断算法成功了损失收敛损失函数应平滑下降并稳定在一个很小的值如1e-6以下。参数聚类对于具有对称性的模型学习到的多个独立链接参数应自动聚类到少数几个值如最近邻链接的β次近邻链接的κ。图7左图在κ̃0时κ被驱动到零就是成功的标志。关联函数匹配终极验证是比较原模型通过拓扑线映射计算和对偶模型用学到的β, κ进行蒙特卡洛模拟的关联函数。如图8所示在中等距离上应有良好重合。差异应随距离增大而增大由于有限尺寸效应和近似对偶的本质。与已知结果交叉验证在κ̃0时恢复的β应与解析公式β -1/2 log(tanh β̃)给出的曲线基本一致。近似对偶的意义 当κ̃ ≠ 0时我们得到的(β, κ)定义的对偶模型是一个有效理论。它可能在某个能标或某个关联长度范围内非常精确地描述了原系统的物理。这本身就有很大价值因为它提供了一个更简单仍然是短程耦合的模型来理解原系统。探究这个有效理论的适用范围随κ̃增大或系统尺寸增大如何退化本身就是一个有趣的物理问题。通过这套方法我们不仅自动化了经典对偶的复现更打开了一扇门系统性地扫描参数空间寻找新的、未知的近似对偶关系或许能帮助我们理解复杂模型背后的有效描述甚至发现新的普适类。这正体现了机器学习与理论物理结合的魅力——它不是要取代物理直觉而是为我们提供前所未有的计算显微镜和探索工具。
机器学习发现统计物理对偶性：从伊辛模型到拓扑线方法
1. 项目概述当统计物理遇见机器学习在统计物理的研究工具箱里对偶性一直是一个既强大又神秘的存在。简单来说它就像是为一个复杂的物理系统找到了另一套完全不同的“语言”来描述它但这套新语言描述的物理本质却完全一样。最经典的例子莫过于二维伊辛模型的Kramers-Wannier对偶一个在高温下无序、低温下有序的自旋系统其高温展开竟然与另一个自旋系统的低温展开在数学上等价并且两者通过一个优美的公式β̃ -1/2 log(tanh β)联系起来。传统上发现这样的对偶关系依赖于物理学家深刻的物理直觉和精巧的数学变换比如高温展开、傅里叶变换或规范理论。然而随着模型变得复杂——比如引入次近邻耦合、多体相互作用或在非平凡格点上——寻找解析对偶变换的难度呈指数级增长。这时一个自然而然的想法是能否让机器来帮我们“发现”这些隐藏的对称性这正是我们这项工作的核心将对偶性发现形式化为一个机器学习优化问题。我们不再依赖先验的解析猜测而是构建一个框架让算法通过分析系统的蒙特卡洛模拟数据自动寻找可能的对偶映射和耦合参数。这项工作的价值在于它为我们探索未知的强关联系统、寻找新的普适类甚至理解拓扑序与对称性保护相提供了一个全新的、数据驱动的视角。对于从事统计物理、凝聚态理论或格点场论研究的同行以及希望将机器学习应用于基础科学发现的探索者本文将详细拆解我们如何将这一想法落地从核心思路到代码实现并分享一路走来的经验与教训。2. 核心思路与算法框架设计2.1 问题形式化什么是对偶性发现在传统机器学习中一个常见任务是“逆伊辛问题”给定一组自旋构型样本{σ_i}推断出产生这些样本的哈密顿量H(σ) -Σ K_ij σ_i σ_j中的耦合常数K_ij。这本质上是一个参数拟合问题。对偶性发现则是一个更宏大的“元问题”。我们不仅不知道耦合参数甚至连描述系统的基本变量、它们所居住的“空间”格点以及可观测量之间的映射关系都是未知的。具体来说我们需要同时学习对偶格点上的耦合参数β̃在对偶描述中相互作用强度是多少对偶格点的几何结构L̃自旋是放在原来的格点上还是其对偶格点如正方格子的对偶仍是正方格子三角格子的对偶是蜂窝格子可观测量映射G(σ̃)如何将原格点上的物理量如自旋关联函数用对偶格点上的变量σ̃表达出来我们的目标是找到一个三元组(β̃, L̃, G)使得用这个三元组定义的对偶理论其所有可观测量的统计分布或至少其低阶矩与原理论完全一致。这就构成了一个损失函数的基础衡量两个理论预测的关联函数之间的差异。2.2 两大技术路径机器学习方法与拓扑线方法在实际操作中我们探索了两种互补的技术路径它们各有优劣适用于不同场景。2.2.1 机器学习方法端到端的黑箱优化这种方法最为直接和通用。我们构建一个可微分的模型其输入是原系统的自旋样本输出是对偶系统的预测关联函数。模型内部包含可学习的参数分别对应β̃、L̃的某种参数化表示例如一个选择哪些链接有效的注意力权重以及映射G例如一个小型神经网络。损失函数通常定义为原系统与对偶系统预测的若干低阶关联函数如一阶、二阶矩之间的均方误差。通过梯度下降如Adam优化器同时优化所有参数当损失收敛到零附近时我们就找到了一个候选对偶。实操心得特征工程的重要性在最初的尝试中我们直接将原始自旋构型输入一个全连接网络效果很差。后来我们发现必须将物理先验知识编码进去。例如我们不是输入所有自旋而是输入一个局部窗口内所有可能链接上的自旋乘积即σ_i σ_j作为特征。这相当于告诉网络“对偶变换很可能只与局部自旋关联有关”。这大大提升了训练的稳定性和可解释性。链接的选择如图2中的7种可能成为了L̃学习的一部分。2.2.2 拓扑线方法利用对称性的白盒优化第二种方法限制更多但更高效、更可解释。它专门用于寻找Kramers-Wannier型的对偶即由Z2自旋翻转对称性所生成的拓扑缺陷线拓扑线所实现的对偶。这类对偶的特点是对偶格点就是原格点的几何对偶如正方格子的对偶仍是正方格子并且可观测量映射G具有确定的形式。具体来说对于原系统的一个Z2对称性算子我们可以定义其对应的拓扑缺陷线。穿过这条线的自旋关联函数会以确定的方式映射到对偶系统的自旋算子上。对于伊辛模型这个映射就是著名的Kramers-Wannier对偶映射对偶自旋σ̃_p位于原格点的 plaquette方格中心其期望值由穿过该 plaquette 边界的原自旋弦算子的乘积给出。在这个强约束下问题大大简化L̃和G已知只剩下对偶耦合参数β̃需要学习。这就将问题规约到了一个增强版的逆伊辛问题。3. 拓扑线方法的数值实现详解鉴于拓扑线方法在效率和物理清晰度上的优势我们将其作为本文深入剖析的重点。下面我们以包含次近邻相互作用的广义伊辛模型为例一步步拆解算法流程。3.1 模型与目标我们考虑的原模型哈密顿量为H̃[σ̃] -β̃ Σ_ĩj̃ σ̃_ĩ σ̃_j̃ - κ̃ Σ_[ĩj̃] σ̃_ĩ σ̃_j̃其中ĩj̃表示最近邻相互作用[ĩj̃]表示对角方向的次近邻相互作用见图6。β̃和κ̃是耦合常数。我们的目标是假设对偶模型也是一个定义在正方格点上的伊辛模型可能也包含次近邻耦合其哈密顿量为H[σ] -β Σ_ij σ_i σ_j - κ Σ_[ij] σ_i σ_j。我们要从原模型的蒙特卡洛样本出发找出最优的(β, κ)使得两个模型在可观测量上匹配。3.2 核心算法基于牛顿法的耦合学习算法的核心思想源于逆伊辛问题中的一个关键等式。对于一个广义伊辛模型H -Σ_{i,δ} K_{i,iδ} σ_i σ_{iδ}利用自旋的离散性σ_i ±1可以推导出一个自洽方程类似于系统的“运动方程”σ_j σ_{jδ} tanh( Σ_{δ} K_{j,jδ} σ_{jδ} ) * σ_{jδ} 这个等式的含义是两点关联函数可以用一个包含所有与σ_j耦合的自旋的非线性函数双曲正切的期望值来表示。如果我们的耦合参数K_{δ}是正确的那么这个等式应该严格成立。算法步骤数据准备使用蒙特卡洛方法如Metropolis-Hastings算法对原哈密顿量H̃进行采样得到大量自旋构型{σ̃_i}。计算对偶关联函数利用已知的拓扑线映射关系即Kramers-Wannier对偶公式从{σ̃_i}样本中计算出对偶自旋{σ_i}的所有偶数阶关联函数。注意在有限系统中奇数阶关联函数恒为零。关键难点我们无法直接得到对偶自旋σ_i的样本只能得到它们的关联函数。而上述自洽方程III.18的右边需要计算像tanh(...) σ_{jδ}这样的非线性关联函数这无法直接从关联函数的集合中简单获得。重构边际概率分布为了解决上述难点我们采用了一个“迂回”但精确的方法。由于耦合K_δ是局域的假设其作用范围为R那么要计算涉及中心自旋σ_j的非线性关联函数我们只需要知道以j为中心、边长为N 2R的窗口内所有自旋的联合概率分布p({σ_a})其中a标记窗口内的N^2个自旋。这个分布p({σ_a})是一个2^{N^2}维的向量因为每个自旋有±1两种状态。窗口内所有可能的自旋乘积的期望值即所有偶数阶关联函数构成了一个线性方程组其系数矩阵由自旋乘积的取值±1决定右边是步骤2中计算出的关联函数值。通过求解这个2^{N^2} × 2^{N^2}的线性方程组III.21我们就可以反解出完整的边际概率分布p({σ_a})。这一步在N较小如3或4时是可行的因为矩阵尺寸512×512或65536×65536在现代计算机上尚可处理。迭代优化耦合参数 a.计算梯度有了p({σ_a})任何非线性函数f({σ_a})的期望值都可以通过直接求和计算f Σ_{σ_a} p({σ_a}) f({σ_a})。利用此公式我们可以计算自洽方程III.18两边的值并定义残差。 b.牛顿法更新我们计算残差关于耦合参数K_δ的雅可比矩阵D_{δ,δ}即梯度公式III.19。然后求解线性方程残差 Σ_{δ} D_{δ,δ} ΔK_{δ}得到耦合参数的更新量ΔK_δ。 c.迭代以学习率ε更新参数K_δ ← K_δ ε ΔK_δ。重复步骤4a-4c直至残差收敛到零附近。3.3 针对次近邻模型的调整与挑战当我们处理包含次近邻耦合κ̃的模型时情况变得更有趣。传统的Kramers-Wannier对偶解析推导会变得异常复杂见附录A.2因为次近邻相互作用在“规范固定”时存在多种等价的连接方式导致对偶模型中会出现非局域或多体相互作用。在我们的数值方法中我们采取了一种实用主义和近似的视角假设我们限制对偶哈密顿量H[σ]的形式只允许它包含最近邻 (β) 和次近邻 (κ) 耦合。这是一个模型选择。目标在这个受限的模型空间内寻找最优的(β, κ)使得其对关联函数的预测与原模型通过拓扑线映射计算出的对偶关联函数最匹配。结果我们发现即使这个受限的对偶模型在严格意义上并非原模型的精确对偶解析推导表明精确对偶可能包含更复杂的相互作用但在一个相当宽的参数范围内特别是当κ̃较小时它提供了一个惊人的优秀近似。图8显示两个模型的关联函数在多个距离上几乎重合。深度解析为什么算法在无序相失效在图7中一个非常显著的现象是当原系统处于无序相β̃ β̃_c时我们的算法无法可靠地恢复出有意义的耦合参数。八个独立链接的耦合值不再聚类为预期的两组β和κ而是变得杂乱无章。 这背后有深刻的物理和算法原因。在无序相系统的关联长度很短自旋涨落很大。我们用于重构边际概率分布p({σ_a})的窗口大小N是固定的。如果关联长度远小于N那么窗口边缘的自旋与中心自旋几乎独立导致我们试图从关联函数反推联合概率分布时问题变得病态ill-posed。矩阵(III.21)接近奇异微小噪声会导致解的巨大波动。此外在无序相系统更接近高温展开其对偶应处于有序相β很大。学习大耦合参数本身就是逆伊辛问题中的难点因为此时概率分布非常尖锐梯度信号很弱。这提醒我们算法的成功应用依赖于物理系统的状态需要谨慎选择参数区域。4. 实操过程从代码到结果4.1 环境准备与依赖安装我们的实现主要基于 Python 的科学计算栈。以下是核心依赖# 核心计算与自动微分 pip install numpy torch # 蒙特卡洛模拟与数据分析 pip install scipy pandas # 可视化 pip install matplotlib seaborn注意为了复现论文中的精确结果需要确保随机数种子固定并且使用双精度浮点数进行计算特别是在求解线性方程组步骤3时以避免数值误差积累。4.2 蒙特卡洛采样模块实现我们采用标准的 Metropolis 算法生成自旋构型。关键在于高效计算能量差。import numpy as np import torch def ising_energy(spins, beta, kappa, latticesquare): 计算包含次近邻耦合的伊辛模型能量。 spins: 二维 numpy 数组或 torch 张量值为 ±1。 beta: 最近邻耦合。 kappa: 次近邻对角耦合。 lattice: 格点类型默认为正方。 # 最近邻相互作用 energy_nn -beta * (np.roll(spins, 1, axis0) np.roll(spins, -1, axis0) np.roll(spins, 1, axis1) np.roll(spins, -1, axis1)) * spins # 次近邻相互作用两个对角线方向 energy_diag -kappa * (np.roll(np.roll(spins, 1, axis0), 1, axis1) np.roll(np.roll(spins, 1, axis0), -1, axis1) np.roll(np.roll(spins, -1, axis0), 1, axis1) np.roll(np.roll(spins, -1, axis0), -1, axis1)) * spins return np.sum(energy_nn) np.sum(energy_diag) def metropolis_sweep(spins, beta, kappa, n_sweeps1000): 执行多次 Metropolis 扫描。 L spins.shape[0] for _ in range(n_sweeps * L * L): i, j np.random.randint(0, L, size2) delta_E 2 * spins[i, j] * ( beta * (spins[(i1)%L, j] spins[(i-1)%L, j] spins[i, (j1)%L] spins[i, (j-1)%L]) kappa * (spins[(i1)%L, (j1)%L] spins[(i1)%L, (j-1)%L] spins[(i-1)%L, (j1)%L] spins[(i-1)%L, (j-1)%L]) ) if delta_E 0 or np.random.rand() np.exp(-delta_E): spins[i, j] * -1 return spins4.3 核心算法边际分布重构与牛顿法优化这是整个项目的核心我们将其分解为几个函数。def compute_dual_correlations_from_original(spins_original, window_size3): 步骤2从原自旋样本通过Kramers-Wannier映射计算对偶自旋在窗口内的所有关联函数。 输入: spins_original - 原格点自旋样本 (L, L) 输出: corrs_dict - 字典键为关联函数索引元组值为估计的期望值。 例如corrs_dict[((0,0), (1,0))] 表示 σ_(0,0) σ_(1,0) L spins_original.shape[0] # 对偶自旋位于原格点plaquette中心索引可视为原格点链接的右上角顶点 # 简化实现将对偶自旋定义为周围四个原自旋的乘积 dual_spins np.ones((L, L)) for i in range(L): for j in range(L): # 注意边界条件这里使用周期性边界 dual_spins[i, j] (spins_original[i, j] * spins_original[(i1)%L, j] * spins_original[i, (j1)%L] * spins_original[(i1)%L, (j1)%L]) # 现在计算 window_size x window_size 窗口内所有对偶自旋的关联函数 # 这里需要计算所有可能的自旋乘积的期望值即所有 2^(N^2) 种组合。 # 为了效率我们通常不会显式计算所有组合而是通过后续的矩阵方程隐式处理。 # 更高效的做法是直接构建所有可能的自旋乘积模式并计算其在样本上的平均值。 # 这部分代码较长核心是遍历所有可能的 α_a ∈ {0,1}计算对应的关联函数。 # ... return corrs_dict def solve_marginal_distribution(corrs_dict, window_size3): 步骤3从关联函数反解边际概率分布 p({σ_a})。 输入: corrs_dict - 上一步计算的所有关联函数。 window_size - 窗口边长 N。 输出: p_vec - 长度为 2^(N^2) 的向量表示所有可能自旋构型的概率。 N2 window_size * window_size total_states 1 N2 # 2^(N^2) # 构建矩阵 M 和向量 b M np.zeros((total_states, total_states)) b np.zeros(total_states) # 为每个自旋构型索引 s (0 到 2^(N^2)-1) 和每个关联函数模式 α (也用一个整数表示) # 填充矩阵 M[s, α] Π_a (σ_a)^(α_a)其中 σ_a 根据构型s确定。 # 向量 b[α] corrs_dict 中对应的关联函数值。 # 这是一个大规模的但稀疏的因为很多关联函数为0线性方程组。 # 实践中由于N3时矩阵大小为512x512可以直接用np.linalg.solve求解。 # ... p_vec np.linalg.solve(M, b) # 确保概率和为1且所有元素非负可进行小幅裁剪和归一化 p_vec np.clip(p_vec, 1e-10, None) p_vec / p_vec.sum() return p_vec def compute_gradient_and_update(K_current, p_marginal, learning_rate0.01): 步骤4计算梯度并用牛顿法更新耦合参数 K。 K_current: 当前耦合参数向量例如 [β, κ, ...] 对应不同方向的链接。 p_marginal: 边际概率分布向量。 # 计算所有需要的期望值左边 σ_j σ_{jδ} 和右边 tanh(...) σ_{jδ} # 以及雅可比矩阵 D_{δ,δ} ∂σ_j σ_{jδ} / ∂K_{δ} # 这需要遍历边际分布中的所有构型计算对应的自旋值和双曲正切值。 left_expectations [] right_expectations [] jacobian [] # ... (具体计算涉及对 p_marginal 的加权求和) # 构建线性系统残差 J * ΔK residual np.array(left_expectations) - np.array(right_expectations) J np.array(jacobian) # 求解 ΔK。注意J可能奇异使用伪逆或添加正则项。 delta_K np.linalg.lstsq(J, residual, rcond1e-6)[0] # 更新参数 K_new K_current learning_rate * delta_K return K_new, np.linalg.norm(residual)4.4 主训练循环与结果可视化将上述模块串联起来形成完整的工作流。def main_training_loop(beta_original, kappa_original, lattice_size12, window_size3, epochs100): 主训练循环。 # 1. 生成原模型数据 print(fGenerating MC samples for (β̃, κ̃) ({beta_original}, {kappa_original})) spins_init np.random.choice([-1, 1], size(lattice_size, lattice_size)) spins_mc metropolis_sweep(spins_init.copy(), beta_original, kappa_original, n_sweeps5000) # 通常需要多个独立样本这里简化为使用一个平衡后的构型进行多次测量。 # 2. 初始化对偶耦合参数猜测 K np.random.randn(8) * 0.1 # 假设8个独立链接参数 losses [] for epoch in range(epochs): # 3. 计算对偶关联函数 dual_corrs compute_dual_correlations_from_original(spins_mc, window_size) # 4. 求解边际分布 p_marginal solve_marginal_distribution(dual_corrs, window_size) # 5. 牛顿法更新 K, loss compute_gradient_and_update(K, p_marginal, learning_rate0.05) losses.append(loss) if epoch % 10 0: print(fEpoch {epoch}: Loss {loss:.6f}, K {K}) if loss 1e-6: print(Converged!) break # 结果分析将学习到的8个参数聚类为β和κ # 对于正方格子的最近邻和次近邻模型8个链接可分为两组 # 4个最近邻链接应趋于相同的值 β4个次近邻链接应趋于相同的值 κ。 beta_learned np.mean(K[[0, 2, 4, 6]]) # 假设索引0,2,4,6对应最近邻 kappa_learned np.mean(K[[1, 3, 5, 7]]) # 假设索引1,3,5,7对应次近邻 print(fLearned parameters: β {beta_learned:.4f}, κ {kappa_learned:.4f}) return beta_learned, kappa_learned, losses运行此代码对于(β̃, κ̃) (0.45, 0.05)我们可能得到类似β ≈ 0.261, κ ≈ 0.019的结果。然后可以计算并绘制两个模型的关联函数进行对比如图8所示验证对偶近似的质量。5. 常见问题、调试技巧与经验总结在实际复现和扩展这项工作的过程中我们遇到了诸多挑战。以下是总结出的关键问题和解决方案。5.1 数值不稳定与病态问题问题表现在求解边际分布p({σ_a})时np.linalg.solve程序报错或解中出现巨大的负值或大于1的值。牛顿法更新步长ΔK异常大导致参数发散。在无序相 (β̃ β̃_c) 结果完全不可信。根源与解决方案矩阵条件数过大方程III.21的矩阵M可能接近奇异。这是因为当窗口内自旋关联很弱时许多关联函数近似线性相关。解决方案使用正则化技术。将np.linalg.solve(M, b)替换为np.linalg.lstsq(M, b, rcond1e-8)或添加 Tikhonov 正则项(M^T M λ I) Δp M^T b其中λ是一个小的正数如1e-6。采样噪声蒙特卡洛样本有限导致关联函数b的估计有误差噪声在反问题中被放大。解决方案增加蒙特卡洛采样量并使用更先进的采样技术如 Wolff 集群算法来减少临界点附近的临界慢化。对多个独立样本计算关联函数并取平均。无序相的根本困难如前所述在无序相问题本身是病态的。此时不应期望算法能给出稳定结果。一个实用的检查点是计算出的边际分布p({σ_a})是否接近均匀分布所有构型概率接近1/2^{N^2}如果是说明系统缺乏有效信息算法失效是预期的。5.2 收敛性与超参数选择问题表现损失函数震荡不下降。学习到的耦合参数β,κ不收敛到稳定值或者收敛到明显错误的解如所有耦合都趋于零。调试技巧学习率ε这是最关键的超参数。起始学习率建议设在0.01到0.1之间。如果损失震荡降低学习率如果下降过慢适当增加。可以实现学习率衰减例如每50轮乘以0.9。初始化耦合参数K的初始值不宜设为零或过大。零初始化可能导致梯度也为零。建议从N(0, 0.1)高斯分布中随机初始化。损失监控除了总损失还应监控自洽方程III.18中每个链接δ对应的残差。这有助于识别是某个特定方向的耦合学习有问题还是整体问题。梯度裁剪在牛顿法求解ΔK后如果||ΔK||过大如 1.0可以将其按比例缩放防止单步更新过大破坏稳定性。5.3 性能优化与扩展性瓶颈瓶颈分析 整个算法最耗时的部分是边际分布求解涉及O(2^{3N^2})复杂度的矩阵构建与求解N3时为512^3量级。N4时矩阵尺寸为65536 x 65536内存和计算需求激增。期望值计算每次迭代都需要遍历所有2^{N^2}个构型来计算非线性期望值和雅可比矩阵计算量为O(2^{N^2} * num_links^2)。优化建议利用对称性对于具有平移、旋转对称性的系统许多链接的耦合K_δ是等价的。可以强制约束这些链接参数相等从而减少待优化参数数量同时雅可比矩阵的维度也大大降低。稀疏性与快速变换矩阵M和计算期望值的求和具有高度结构化的模式。探索是否能用快速沃尔什-哈达玛变换来加速计算因为自旋取值为±1所有计算本质上是在布尔超立方体上的求和。替代方案对于更大的N直接矩阵求逆不可行。可以考虑使用变分推断或训练一个小型神经网络来直接从可测量的关联函数估计非线性期望值tanh(...) σ从而绕过显式求解p({σ_a})这一步。这正是原文未来工作展望中提到的方向。5.4 物理结果解读与验证如何判断算法成功了损失收敛损失函数应平滑下降并稳定在一个很小的值如1e-6以下。参数聚类对于具有对称性的模型学习到的多个独立链接参数应自动聚类到少数几个值如最近邻链接的β次近邻链接的κ。图7左图在κ̃0时κ被驱动到零就是成功的标志。关联函数匹配终极验证是比较原模型通过拓扑线映射计算和对偶模型用学到的β, κ进行蒙特卡洛模拟的关联函数。如图8所示在中等距离上应有良好重合。差异应随距离增大而增大由于有限尺寸效应和近似对偶的本质。与已知结果交叉验证在κ̃0时恢复的β应与解析公式β -1/2 log(tanh β̃)给出的曲线基本一致。近似对偶的意义 当κ̃ ≠ 0时我们得到的(β, κ)定义的对偶模型是一个有效理论。它可能在某个能标或某个关联长度范围内非常精确地描述了原系统的物理。这本身就有很大价值因为它提供了一个更简单仍然是短程耦合的模型来理解原系统。探究这个有效理论的适用范围随κ̃增大或系统尺寸增大如何退化本身就是一个有趣的物理问题。通过这套方法我们不仅自动化了经典对偶的复现更打开了一扇门系统性地扫描参数空间寻找新的、未知的近似对偶关系或许能帮助我们理解复杂模型背后的有效描述甚至发现新的普适类。这正体现了机器学习与理论物理结合的魅力——它不是要取代物理直觉而是为我们提供前所未有的计算显微镜和探索工具。
