2D伊辛模型多重网格蒙特卡洛算法实现及粗化步骤咨询

算法可行性结论 #

多重网格蒙特卡洛完全可以应用于简单2D Ising模型，并非仅适配XY模型、高斯模型这类连续自由度模型。Sokal在提出该方法框架后，后续就有针对离散自旋系统的落地验证，不存在原理层面的障碍。
你遇到的粗化步骤卡点，本质是直接照搬了确定性多重网格求解偏微分方程的线性算子逻辑——确定性场景下网格值是连续的，加权平均、线性插值都可以直接用，但伊辛自旋是取值±1的离散变量，硬套连续场景的算子必然出现粗网格值非二值、反粗化赋值冲突的问题。
确定性多重网格算法的参考流程如下：

注意：实现时不需要刻意追求跨网格映射的信息无损，多重网格方法的加速效果恰恰来自于不同尺度下对长程、短程模式的分别处理，只要算子设计满足马尔可夫链细致平衡条件，采样分布就不会出现偏差。

可落地的粗化/反粗化实现方案 #

目前成熟的、能保证粗网格自旋为±1、反粗化无歧义的实现方案主要有两类：

• 固定块匹配重整化群粗化方案
  这是最容易用Python实现的入门方案：
  1. 粗化（限制）步骤：将细网格划分为互不重叠的2×2自旋块，每个块对应1个粗网格自旋，粗自旋取值为块内多数自旋的符号；如果块内正负自旋数量相等，等概率随机取+1或-1即可。
     2×2块粗化的快速numpy实现参考：

  import numpy as np
  def restrict(fine_grid):
      # 输入为尺寸2N*2N的细网格自旋数组，取值±1
      n = fine_grid.shape[0] // 2
      # 拆分为N*N个2*2子块
      blocks = fine_grid.reshape(n, 2, n, 2).transpose(0, 2, 1, 3)
      block_sum = blocks.sum(axis=(-1, -2))
      coarse_spin = np.sign(block_sum)
      # 处理块内自旋和为0的情况，随机赋值±1
      zero_mask = (coarse_spin == 0)
      coarse_spin[zero_mask] = np.random.choice([-1, 1], size=zero_mask.sum())
      return coarse_spin
  

  2. 粗网格有效哈密顿量：粗网格自旋间的耦合常数可以通过块自旋重整化的递推关系计算，初期调试可以先取粗网格耦合为细网格耦合的2倍，验证流程跑通后再替换为理论递推值。
  3. 反粗化（延拓）步骤：不要做确定性的自旋赋值，把粗网格自旋作为对应细网格块的背景外场，对细网格自旋做2-3次热浴更新（也就是原方法中替换Gauss-Seidel的松弛步骤）时，在自旋的条件概率计算里额外加上对应粗网格自旋的耦合项，自然不存在赋值歧义，也能自动满足细致平衡。
• 团簇匹配粗化方案
  如果追求临界点附近的更高效率，可以结合Swendsen-Wang团簇算法设计粗化规则：
  1. 粗化步骤：先在细网格上按Swendsen-Wang规则生成连通自旋团簇，每个连通团簇对应1个粗网格自旋，粗自旋取值和团簇自旋方向一致。
  2. 反粗化步骤：团簇内所有细网格自旋直接继承对应粗网格的自旋值，完全没有赋值冲突，额外做1-2次单点热浴松弛消除块边界的偏差即可。
     该方案的缺点是粗网格不再是规则正方形网格，需要额外维护粗自旋的邻接关系，代码实现复杂度更高。

实现避坑提示 #

• 不要直接照搬确定性多重网格的线性插值、加权限制算子，离散自旋模型的跨网格算子必须满足细致平衡，否则采样结果会出现系统性偏差。
• 每个网格层的松弛（热浴更新）不需要做太多次扫描，通常2-3次扫描就足够消除该尺度对应的高频模式偏差，扫描次数过多反而会降低整体效率。
• 初期实现可以先从高温、远离临界点的参数开始调试，确认磁化率、内能等物理量的计算结果和标准单点热浴/Metropolis算法结果一致后，再测试临界点附近的慢模式加速效果。

内容的提问来源于stack exchange，提问作者pter26