从样本数据估计费舍尔信息矩阵：MCMC与Lanczos方法在相变探测中的应用

1. 项目概述：从样本数据中“看见”统计流形的曲率

费舍尔信息矩阵（Fisher Information Matrix, FIM）在统计推断和信息几何中，扮演着“度量衡”的角色。简单来说，它衡量了当你微调模型参数时，其对应的概率分布会发生多大变化。这个变化率，或者说敏感度，直接决定了你从数据中估计这些参数的精度上限——这就是著名的克拉美-罗下界。在物理领域，特别是研究多体系统的相变时，系统的基态或热态概率分布会随着某个物理参数（如温度、磁场强度）的连续变化而发生剧烈改变，这种剧变往往就体现在FIM的奇异性（如发散或尖锐峰值）上。因此，FIM成为了一个探测相变边界的强大理论工具。

然而，理论上的FIM定义需要我们知道精确的概率分布形式。在现实世界的物理模拟或实验中，我们往往无法获得完整的解析分布，只能通过数值模拟（如蒙特卡洛采样）或实验测量，得到一系列来自该分布的样本。这就引出了本文的核心挑战与主题：如何仅从有限的、离散的样本数据中，稳健且高效地估计出那个理论上连续、光滑的FIM？

传统的做法可能需要先进行密度估计，再求导计算，这在复杂高维系统中几乎不可行。我们这里探讨的，是两条更“聪明”的路径：一条基于经典的马尔可夫链蒙特卡洛（MCMC），适用于经典统计模型（如伊辛模型）；另一条基于Lanczos算法，用于求解量子多体系统的低能本征态，进而得到其概率分布。这两种方法都绕开了显式的密度函数，直接从样本或态矢量出发，通过巧妙的统计处理或线性代数计算，逼近FIM。无论你是从事统计物理、量子计算还是机器学习的研究者，只要面临从样本中推断模型内在结构的问题，理解这些估计方法都将大有裨益。接下来，我将拆解这两种方法的核心思路、实操细节，并分享在实现过程中积累的经验与避坑指南。

2. 核心原理：为何FIM能揭示相变，以及如何从样本逼近它

2.1 费舍尔信息矩阵的几何直观与物理意义

首先，我们得抛开复杂的公式，从几何角度理解FIM。想象一个“统计流形”，这个流形上的每一个点，都对应着一组特定参数下的概率分布。当我们在这个流形上移动时（即改变参数），分布之间的距离如何度量？FIM就是这个流形上的黎曼度量张量。它定义了参数空间局部邻域内的“距离”概念：两个无限接近的参数点所对应的分布之间的KL散度，可以近似表示为参数差值的二次型，而其中的系数矩阵就是FIM。

在物理相变点附近，系统处于临界状态，其概率分布对参数的微小变化极端敏感。对应到统计流形上，就是流形的曲率在此处变得极大。因此，FIM的矩阵元（或其标量不变量，如迹或行列式）会在相变点附近出现尖锐的峰值或发散。这就是我们用FIM探测相变的根本原因。它不依赖于任何序参量的先验定义，是一种直接从系统概率分布中提取的、内禀的探针。

2.2 从KL散度到样本估计：核心桥梁公式

理论上的FIM定义为对数似然函数得分向量（Score）的协方差矩阵。但对于我们从样本估计的场景，一个更实用的出发点是它与KL散度的关系。对于两个无限接近的参数点 $\lambda$ 和 $\lambda + d\lambda$，它们的KL散度可以展开为： $D_{KL}(P_{\lambda} || P_{\lambda + d\lambda}) \approx \frac{1}{2} d\lambda^T G(\lambda) d\lambda$ 其中 $G(\lambda)$ 就是FIM。

我们的目标就是从样本中估计这个二次型。假设我们有两个接近的参数点 $\lambda^+$ 和 $\lambda^-$，其中心点为 $\lambda^0$，参数差为 $\delta\lambda = \lambda^+ - \lambda^-$。一个关键的关系式是（推导详见附录，但思想很重要）： $g(\lambda^0)(\delta\lambda) \approx 4 \times \left[ 1 - \mathbb{E}{x \sim q} \left( \frac{2}{\sqrt{r_x} + 1/\sqrt{r_x}} \right) \right]$ 这里 $g(\lambda^0)(\delta\lambda)$ 是FIM在方向 $\delta\lambda$ 上的二次型值，$q$ 是 $P{\lambda^+}$ 和 $P_{\lambda^-}$ 的混合分布，而 $r_x = P_{\lambda^+}(x) / P_{\lambda^-}(x)$ 是似然比。

注意：这个公式的美妙之处在于，它将FIM的估计问题，转化为了对似然比 $r_x$ 的估计问题。我们不需要知道完整的概率 $P_{\lambda}(x)$，只需要知道它们在样本点上的比值。这对于MCMC这类只能给出未归一化概率 $\tilde{P}{\lambda}(x) \propto P{\lambda}(x)$ 的方法至关重要。

2.3 两种技术路线的分野：MCMC vs. Lanczos

基于上述原理，我们发展出两条技术路线，它们的分野根植于所研究系统的本质：

1. 基于MCMC的经典统计物理路径：适用于经典格点模型，如伊辛模型（Ising400, IsNNN400）。其概率分布是经典的吉布斯分布 $P_{\lambda}(x) \propto \exp(-\beta H_{\lambda}(x))$。MCMC（如Metropolis-Hastings算法）可以高效地从该分布中采样，并且我们能轻松计算每个样本 $x$ 对应的未归一化概率$\tilde{P}{\lambda}(x) = \exp(-H{\lambda}(x))$（这里设 $\beta=1$，温度影响可吸收进哈密顿量）。利用公式(G3)和(G4)，我们可以用两个邻近参数点的样本和未归一化概率，来估计似然比 $r_x$ 以及关键的配分函数比值 $\sqrt{Z_{\lambda^-}/Z_{\lambda^+}}$，最终拼凑出FIM。

2. 基于Lanczos的量子多体系统路径：适用于量子晶格模型，如费米-哈伯德模型（Hubbard12）或海森堡模型（XXZ300）。我们关心的是其基态（或有限温吉布斯态）在计算基矢下的概率分布 $p_z(\lambda) = |\langle z|\psi_0(\lambda)\rangle|^2$。这里，Lanczos算法（或其变种如隐式重启的Lanczos-Arnoldi方法）被用来数值求解大型稀疏哈密顿量的前几个本征态。一旦得到基态波函数 $|\psi_0(\lambda)\rangle$，我们就能直接计算出任意计算基矢 $|z\rangle$ 下的概率幅，进而得到精确的 $p_z(\lambda)$。有了精确的概率分布，FIM就可以通过有限差分法直接计算（公式G10）。

实操心得：选择哪条路径，首先取决于你的系统是经典的还是量子的。对于经典系统，MCMC是唯一现实的选择，因为状态空间太大无法枚举。对于中小规模的量子系统（~24个格点），Lanczos直接对角化是可行的；对于更大的一维量子系统（如300个格点），则需采用密度矩阵重整化群（DMRG）等张量网络方法先获得近似基态（MPS），再从中采样或计算概率。

3. 实操详解（一）：基于MCMC的FIM估计（以Ising400为例）

让我们以经典的二维伊辛模型（Ising400，即20x20格点，400个自旋）为例，深入MCMC估计FIM的每一步。我们的参数是温度 $T$（或归一化后的 $\lambda = T/4$），目标是计算FIM随温度的变化，并在临界温度 $T_c$ 附近观察到峰值。

3.1 数据生成：平衡采样与温度扫描

生成可靠样本是估计的基石。对于伊辛模型，我们采用模拟退火结合并行回火的思路来确保每个温度下的样本都充分平衡。

1. 初始化与预热：从高温随机状态开始，运行足够多的MCMC步数使系统在初始温度下平衡。每个MCMC“步”通常指尝试更新所有格点一次（400个子步）。
2. 温度扫描与自适应步长：从高到低扫描一系列温度点（例如1000个）。在每一个新温度 $T_i$ 下，不能直接从上一个温度的末态开始演化，因为温度突变可能导致系统陷入非平衡。我们的策略是：
   • 根据当前温度与临界温度 $T_c$ 的接近程度，动态调整MCMC步数。定义一个函数 $d(T) = 1 + (0.1 + (T - T_c)^2)^{-1}$。接近 $T_c$ 时，关联长度发散，弛豫时间变长，因此需要更多步数 ($\lfloor 30d \rfloor$) 来重新平衡。
   • 平衡后，采样一个构型。然后继续运行 $\lfloor 9d \rfloor$ 步以降低样本间的自相关性，再次采样。这样每个温度点我们获得两个（近似）独立的样本。
3. 降低自相关：在采样时，额外引入随机全局自旋翻转和晶格平移操作，进一步破坏样本间的时空关联。

注意事项：临界点附近的采样是最大的挑战。d(T)函数中的常数（如0.1）需要根据具体模型和尺寸微调。一个实用的检查方法是计算不同样本间的能量自关联函数，确保在采样间隔内衰减到接近零。对于Ising400，我们使用了70条独立的MCMC链进行平均，以降低随机误差。

3.2 FIM估计的核心步骤：从样本到矩阵元

假设我们已经获得了参数点 $\lambda^+$（对应 $T^+$）和 $\lambda^-$（对应 $T^-$）处的两组样本 ${x^+i}$ 和 ${x^-j}$，以及每个样本的未归一化概率 $\tilde{P}{\lambda}(x) = \exp(-(H{\lambda}(x) - E_{\lambda}))$（减去平均能量 $E_{\lambda}$ 是为了数值稳定性）。

1. 估计配分函数比值：这是关键且容易出问题的一步。根据公式(G4)，我们有两种估计方式： $\frac{Z_{\lambda^-}}{Z_{\lambda^+}} \approx \frac{1}{N^+} \sum_{i=1}^{N^+} \frac{\tilde{P}{\lambda^-}(x^+i)}{\tilde{P}{\lambda^+}(x^+i)} \quad \text{或} \quad \left( \frac{1}{N^-} \sum{j=1}^{N^-} \frac{\tilde{P}{\lambda^+}(x^-j)}{\tilde{P}{\lambda^-}(x^-_j)} \right)^{-1}$ 理论上二者相等，但由于样本有限，估计值会有差异。一个稳健的做法是取两者的几何平均：$\hat{R} = \sqrt{\text{Est}_1 \times \text{Est}_2^{-1}}$。这有助于降低方差。

2. 计算样本似然比：对于来自混合分布 $q$ 的每个样本 $x$（我们可以将 $\lambda^+$ 和 $\lambda^-$ 的样本合并视为来自 $q$），计算： $\hat{r}x = \sqrt{ \frac{\tilde{P}{\lambda^+}(x)}{\tilde{P}_{\lambda^-}(x)} \times \hat{R} }$

3. 计算FIM二次型：代入公式(G2)的离散版本： $\hat{g}(\lambda^0)(\delta\lambda) \approx \frac{4}{N} \sum_{k=1}^{N} \left( 1 - \frac{2}{\sqrt{\hat{r}{x_k}} + 1/\sqrt{\hat{r}{x_k}}} \right)$ 其中 $N$ 是用于估计的合并样本总数。

4. 组装完整FIM：为了得到FIM矩阵 $G$，我们需要估计其在各个参数方向上的分量。对于单参数情况（如Ising400），$\delta\lambda$ 就是温度差，上述估计直接给出 $G_{TT}$。对于多参数情况（如IsNNN400有两个参数），我们需要选择一组线性无关的微小位移向量 ${\delta\lambda^{(µ)}}$（例如沿坐标轴方向），对每个方向执行上述估计，然后通过求解线性方程组来反推出矩阵 $G$ 的各个元素。

3.3 有限尺寸标度与误差分析

图7展示了Ising模型在不同尺寸（10x10, 20x20, 40x40）下的FIM（除以自旋数 $N$ 以进行标度对比）。我们可以看到，在低温区，标度后的FIM曲线发生了漂亮的“标度坍塌”，这是临界现象的特征。然而，有限尺寸效应使得FIM峰的位置 $T_{peak}$ 与热力学极限下的理论临界温度 $T_c$ 存在系统偏差。对于20x20的格点，我们的计算显示 $T_{peak}$ 略高于 $T_c$。

避坑技巧：为了可靠地估计峰值位置及其误差，我们进行了多次独立模拟（每次使用不同的随机数种子）。对于每个系统尺寸，我们计算了11次独立FIM估计的平均值和标准差（图7中的彩色条带）。通过双样本t检验，我们确认了 $T_{peak}$ 与 $T_c$ 的差异具有统计显著性。这提醒我们，用有限尺寸系统的FIM峰值外推热力学极限的临界点时，必须考虑并修正这种系统偏差。

4. 实操详解（二）：基于Lanczos的FIM估计（以Hubbard12/FIL24为例）

对于量子多体系统，我们无法直接采样，但可以通过数值对角化获得精确的基态波函数。Lanczos算法是处理大型稀疏厄米矩阵部分本征值问题的利器。

4.1 基态求解与概率分布构建

我们的目标是计算在64x64的参数网格 $\lambda = (\lambda_0, \lambda_1)$ 上，系统的基态 $|\psi_0(\lambda)\rangle$。

1. 调用ARPACK-NG：我们通过SciPy的scipy.sparse.linalg.eigsh函数调用ARPACK-NG库，使用隐式重启的Lanczos-Arnoldi方法，求解每个 $\lambda$ 点哈密顿量 $H(\lambda)$ 的前4个最低本征态 ${|\psi_j(\lambda)\rangle}_{j=0}^3$ 和本征值 ${E_j(\lambda)}$。
2. 构建低温概率分布：严格来说，我们想要的是 $T=0$ 的基态概率 $|\langle z|\psi_0\rangle|^2$。但数值上，基态与第一激发态可能存在偶然简并或能隙极小。为了数值稳定性，我们构建一个极低温度（$\beta = 10^7$）的截断玻尔兹曼分布： $p_z(\lambda) = \frac{1}{C(\lambda)} \sum_{j=0}^{3} |\langle z|\psi_j\rangle|^2 e^{-\beta (E_j(\lambda)-E_0(\lambda))}$ 当能隙远大于 $1/\beta = 10^{-7}$ 时，该分布几乎完全由基态主导。$C(\lambda)$ 是归一化常数。
3. 处理简并：在Hubbard12模型的参数空间边界（$\lambda_0=0$ 线）上，我们发现了18个点存在数值简并（基态与第一激发态能量差在双精度下不可区分）。这导致Lanczos返回的“基态”是真实基态和激发态的任意线性组合，使得概率分布 $p_z(\lambda)$ 在这些点附近高度不确定。幸运的是，这些点位于参数空间边界，对内部相图影响有限，但需要在后续分析中留意。

4.2 维度约减技巧：应对指数爆炸

24个格点的希尔伯特空间维度是 $2^{24} \approx 1.68\times 10^7$，存储所有基矢系数不现实。我们利用对称性大幅压缩问题规模。

1. 粒子数守恒（Hubbard12）：在固定填充数（半满）下，系统总自旋向上和向下的粒子数分别守恒。这允许我们将问题限制在对应的子空间 $H_0$ 内。对于24位比特串，要求前12位（代表自旋向上）中恰好有6个1，后12位（自旋向下）也恰好有6个1。这样的比特串数量为 $\binom{12}{6}^2 = 924^2 = 853,776$，远小于全空间。我们预先建立索引与比特串的查找表，并在该约减基底下构建哈密顿量的稀疏矩阵。

2. 晶格对称性与项链计数（FIL24）：对于某些模型（如FIL24），其哈密顿量具有非正的非对角元（stoquastic性质），且晶格本身具有对称性（如12个格点的二面体群 $D_{12}$ 对称性）。我们可以将对称性相关的比特串归为同一个“轨道”。基态系数在同一个轨道内所有比特串上是相同的。因此，我们只需对轨道进行索引。12个格点、4种“颜色”（由自旋和轨道占据情况决定）的“手镯”（可翻转的项链）数量是A032275(12)=704,370。这样，我们需要用Lanczos对角化的矩阵维度就从1600万降到了70万，计算量和内存需求大大降低。

实操心得：在构建约减基底的哈密顿量矩阵时，预计算是关键。对于每个参数 $\lambda$，哈密顿量矩阵的非零元素模式（位置）是固定的，只有数值随 $\lambda$ 变化。因此，我们可以预先计算好所有非零元素的位置索引，在循环每个 $\lambda$ 点时，只需快速填充数值即可，这能极大提升计算效率。

4.3 从精确分布计算FIM与生成样本数据集

一旦获得了每个 $\lambda$ 点上的精确概率分布 ${p_v(\lambda)}$（其中 $v$ 是约减基底的索引），后续工作就变得直接。

1. 计算真实FIM（Ground Truth）：我们采用有限差分法。例如，为了估计 $G_{00}$ 在参数点 $((l_0+1/2)/r, l_1/r)$ 的值，使用公式(G10)： $[g]{00} \approx 8r^2 \left{ 1 - \sum_v \sqrt{p_v(l_0/r, l_1/r) \cdot p_v((l_0+1)/r, l_1/r)} \right}$ 这个公式本质上是计算相邻两点概率分布的Bhattacharyya系数（一种相似度度量），其与1的差距反映了分布的变化率，即FIM。类似地可以估计 $G{11}$ 和 $G_{01}$。最后，通过插值将这些在不同中心点的估计值统一到网格中心点 $M''$ 上，得到完整的FIM张量场。

2. 生成样本数据集 $D$：为了模拟“实验测量”，我们从每个 $\lambda$ 点的概率分布 ${p_v(\lambda)}$ 中，随机抽取140个样本。具体步骤是：首先根据概率分布抽样得到轨道索引 $v$，然后在该轨道 $o_v$ 内均匀随机选取一个具体的比特串 $x$。这样，我们就得到了一个包含 $64\times64\times140 \approx 57.3$ 万个 $(\lambda, x)$ 对的数据集。将其90%随机划分为训练集 $D_{train}$，10%作为测试集用于最终评估。

注意事项：有限差分法估计的FIM并非精确值，但对于在离散参数网格上评估估计方法的性能而言，它比真实的连续FIM更合适。因为真实的FIM可能在网格点之间有一个非常尖锐甚至发散的峰，如果峰宽小于网格间距，仅凭网格点上的值可能会完全错过它。有限差分估计反映了在给定网格分辨率下我们能捕捉到的变化。

5. 实操详解（三）：DMRG在更大规模量子系统中的应用（XXZ300）

对于更大规模的一维量子系统（如300个格点的XXZ模型），精确对角化（Lanczos）不再可行。这时，密度矩阵重整化群（DMRG）成为获取近似基态（以矩阵乘积态MPS表示）的标准工具。

5.1 DMRG基态制备流程

我们的目标是获得64x64参数网格上XXZ模型的基态MPS。直接对每个点从头运行DMRG成本高昂，且可能陷入局部极小。我们采用了一个鲁棒的初始化与优化流程：

1. 多初态并行尝试：选取三个初始态：两个最大键维为16的随机MPS，以及一个与系统边缘钉扎项对齐的Néel态（如|0101...01>）。分别用较小的键维（如16）进行1次DMRG扫描优化。
2. 淘汰与选择：计算优化后态之间的重叠（点积）和能量。如果某个态与另一个更低能量的态几乎平行（点积接近±1），则舍弃它；否则，舍弃能量最高的态。保留两个最有希望的候选态。
3. 逐步增加键维优化：对剩余的两个态，进行多轮（如5轮）DMRG扫描，并逐轮增加最大允许键维（如32, 64, 128, 256, 512），以逐步逼近更精确的基态。
4. 对称性投影与最终遴选：将得到的MPS投影到 $\otimes_i Z_i$ 算子的本征态上（因为XXZ模型在某些参数区具有Z2对称性）。舍弃范数过小或与其他态过于接近的投影。选择投影后能量最低的态作为当前最佳基态。
5. 精细优化与激发态试探：对最佳基态进行更多轮（如10轮）精细DMRG扫描以收敛。然后，尝试使用激发态DMRG（通常通过将已找到的态作为“禁止”子空间），看是否能找到能量更低的态。如果找到，则替换当前基态并再次精细优化。

避坑技巧：DMRG的结果严重依赖于初始化和超参数（如扫描次数、截断误差、键维）。图8展示了我们验证DMRG结果的过程：对比不同随机种子运行的结果（左上图显示能量差异），监控最大键维的使用情况（中上图），并检查最后一轮扫描的能量提升是否已可忽略（右上图）。强烈建议对关键参数点用不同随机种子和设置重复计算，以确保结果的可靠性和可重复性。

5.2 从MPS获取概率与样本

获得基态MPS $|\psi\rangle$ 后，计算其在计算基矢下的概率 $p_z = |\langle z|\psi\rangle|^2$ 是直接的，但需要收缩整个张量网络，对于300个格点，这仍然是 $2^{300}$ 量级的求和，不可能全部计算。但我们只需要采样。

1. 概率采样：对MPS进行从左到右的逐点采样。在第一个格点，根据约化密度矩阵计算处于|0>或|1>的概率，依此概率随机采样该格点状态。将该状态固定，更新MPS，然后对下一个格点重复此过程，直至采样出整个比特串。这可以高效地生成服从 $p_z$ 分布的样本。
2. 计算FIM：有了MPS，我们可以用类似有限差分的思想，但更精确地计算FIM。不需要生成大量样本再估计，可以直接计算两个邻近参数点 $\lambda$ 和 $\lambda+\delta\lambda$ 的基态MPS之间的保真度 $F = |\langle \psi(\lambda)|\psi(\lambda+\delta\lambda)\rangle|$。对于纯态，FIM与保真度的关系为：$1 - F \approx \frac{1}{8} \delta\lambda^T G(\lambda) \delta\lambda$。通过计算不同方向 $\delta\lambda$ 下的保真度，即可反解出FIM。这比基于样本的估计更精确，但需要有能力计算不同参数点MPS之间的重叠积分。

6. 常见问题、挑战与解决方案实录

在实际操作中，从样本估计FIM会遇到诸多挑战。以下是我在复现和实验过程中遇到的一些典型问题及其解决思路。

6.1 MCMC估计中的方差与偏差问题

问题：在估计配分函数比值 $\hat{R}$ 和似然比 $\hat{r}x$ 时，特别是当 $P{\lambda^+}$ 和 $P_{\lambda^-}$ 重叠很小时，估计量的方差会变得极大，导致FIM估计不稳定。

排查与解决：

1. 诊断：观察 $\hat{r}_x$ 的分布。如果大量样本的 $\hat{r}_x$ 接近0或无穷大，说明两个分布差异太大，有限差分步长 $\delta\lambda$ 可能过大了。
2. 调整步长：缩小 $\delta\lambda$，使两个分布更接近。但这会使得FIM的有限差分估计本身误差增大。需要在方差和差分误差之间取得平衡。一个经验法则是让 $\mathbb{E}[\hat{r}_x]$ 在0.1到10之间。
3. 重要性采样：公式(G4)本质上是一种重要性采样估计。如果使用 $\lambda^+$ 的样本估计 $Z_{\lambda^-}/Z_{\lambda^+}$，那么权重 $w_i = \tilde{P}_{\lambda^-}(x^+i) / \tilde{P}{\lambda^+}(x^+_i)$ 的方差可能很大。可以尝试使用桥接采样或热力学积分方法，在 $\lambda^-$ 和 $\lambda^+$ 之间插入中间点，逐步估计。
4. 多次独立运行：对于关键参数点（如相变点附近），进行多次独立MCMC模拟，分别计算FIM，然后取平均值和标准差，用误差棒来反映估计的不确定性。

6.2 Lanczos对角化的数值稳定性与简并处理

问题：对于某些参数点，Lanczos算法可能无法收敛，或返回的本征值/态不准确，尤其是在存在简并或近简并的情况下。

排查与解决：

1. 检查收敛性：ARPACK的eigsh函数会返回收敛信息。务必检查返回的info参数。对于未收敛的情况，可以增加迭代次数 (maxiter) 或重启次数。
2. 处理近简并：当基态与第一激发态能隙很小时，Lanczos可能难以区分它们。我们的解决方案是计算前 $k>1$ 个本征态（如k=4），并用低温玻尔兹曼分布混合它们。这虽然“污染”了纯基态，但给出了一个在能隙小于数值精度时更稳健的概率分布。关键是要意识到，在这种情况下，系统本身可能处于量子相变点或具有拓扑简并，FIM的定义本身可能需要修正（例如考虑简并子空间）。
3. 验证基态：对于关键点，可以用不同的初始向量或不同的对角化方法（如全对角化用于小系统）进行交叉验证。对于DMRG，则对比不同随机种子和优化路径的结果。

6.3 高维参数空间与计算成本

问题：对于多参数模型（如2D参数空间），需要在密集网格上计算FIM，计算量随参数维度指数增长。

排查与解决：

1. 自适应网格：不要在整个参数空间使用均匀网格。可以先在粗网格上运行，定位FIM可能较大的区域（相变边界），然后在这些区域进行网格加密。这需要算法能动态调整采样点。
2. 并行计算：不同参数点 $\lambda$ 的计算是完全独立的，这是令人愉快的并行。可以将参数网格任务分配到多个CPU核心甚至多个计算节点上。我们的Hubbard12和FIL24计算就消耗了约3000 CPU小时，没有并行化是不可想象的。
3. 利用对称性：如前所述，利用模型的对称性（粒子数守恒、点群对称性）可以极大降低每个点上的计算成本（希尔伯特空间维度）。这是实现中等规模系统精确计算的关键。

6.4 从FIM到相图：峰值识别与边界提取

问题：得到了FIM标量场（如 $\text{Tr}(G)$ 或 $\det(G)$ ）后，如何自动、鲁棒地识别峰值并勾勒相变边界？

排查与解决：

1. 平滑与去噪：由于估计误差，FIM场可能很嘈杂。在寻找峰值前，应用一个高斯滤波器进行适度的平滑是必要的。但要注意，平滑可能模糊尖锐的峰或移动峰的位置。
2. 峰值检测算法：使用成熟的峰值检测算法（如scipy.signal.find_peaks），并设置合适的参数：prominence（突出度）用来过滤掉小的波动，width（宽度）可以过滤掉过于尖锐的噪声尖峰，distance（距离）可以防止在同一个宽峰上检测到多个峰。
3. 二维边界提取：对于二维参数空间，FIM的峰值可能形成一条“脊线”。可以将FIM场视为地形高度，相变边界就是山脊线。可以使用图像处理中的骨架化或分水岭算法，或者更简单地，设定一个阈值，将FIM值高于阈值的区域标记出来，其边界即为推测的相变线。阈值可以选择为FIM最大值的某个比例（如50%）。
4. 与理论/其他方法对比：始终将自动识别的结果与已知的理论相图或其他探测方法（如序参量、纠缠熵）的结果进行对比校准。这是验证FIM估计方法有效性的最终标准。

7. 与其他方法的对比与定位

我们的工作（ClassiFIM）聚焦于从样本直接估计FIM。在量子多体相变的无监督机器学习领域，还有其他著名方法，如SPCA和W方法（Confusion Scheme）。理解它们的异同有助于定位我们的工作。

7.1 与SPCA的对比

SPCA是一种应用于经典阴影数据的核PCA方法，旨在将参数空间中的点嵌入到一个低维空间，使得不同相的点形成不同的簇。

• 数据需求不同：SPCA设计用于处理经典阴影数据，即每个参数点 $\lambda$ 对应 $T$ 个随机泡利基测量结果。而ClassiFIM处理的是计算基矢测量得到的比特串样本。这是两种不同类型的数据。
• 目标不同：SPCA的目标是获得一个揭示相结构的嵌入，其物理意义（如与FIM的关系）不直接明确。ClassiFIM的目标是直接估计FIM，这是一个有明确信息几何和统计意义的量。
• 内核适配：为了将SPCA应用于我们的比特串数据，我们提出了一个适配的内核公式(I5)。但需要警惕，这个内核与原始SPCA用于阴影数据的核(I3)在数学上并不等价，如附录I.3中的例子所示。它们捕捉的是概率分布之间不同的相似性度量。

7.2 与W方法（Confusion Scheme）的对比

W方法通过训练在错误标签数据上的分类器，其准确率会形成一个“W”形曲线，中间峰值对应相变点。

• 数据复杂度：在Kitaev链的对比中，W方法使用了截断纠缠谱（前10个约化密度矩阵本征值）作为输入数据。如图9所示，这种数据本身已经非常强大，甚至无需机器学习，直接计算其分布的FIM就能清晰显示相变点。相比之下，ClassiFIM使用的比特串测量数据要“弱”得多，它只包含局域粒子占据信息，不包含任何纠缠信息。因此，ClassiFIM在更简单、更易获取的数据上完成了任务。
• 鲁棒性与范围：W方法对参数扫描范围敏感。在我们的测试中，当参数范围从 $[-4, 0]$ 扩展到对称的 $[-4, 4]$ 时，由于数据对称性，W方法可能失效（产生V形而非W形曲线）。而ClassiFIM在整个范围内都能稳定工作，并正确反映出对称性。
• 输出形式：W方法输出的是一个标量曲线（分类准确率），其峰值位置指示相变。ClassiFIM输出的是FIM张量场，它不仅能指示相变位置（峰值），还能通过其特征向量指示相变驱动的主导参数方向，并提供更多的几何信息。

7.3 定量比较的挑战与新指标：PeakRMSE

由于不同方法处理的数据类型和目标输出不同，进行直接的定量比较非常困难。为此，我们提出了一个间接的、面向任务的评估指标：峰值均方根误差（PeakRMSE）。

核心思想：无论方法输出的是FIM曲线、分类准确率曲线还是嵌入后的聚类，我们最终关心的都是它们预测的相变边界。因此，我们可以在一维参数切片上，比较方法预测的峰值位置 $\hat{\lambda}{peak}$ 与已知的“真实”峰值位置 $\lambda^*{peak}$（通常来自有限尺寸精确对角化或高精度计算）。

计算步骤：

1. 对于参数空间的每个一维切片（固定 $\lambda_1$，扫描 $\lambda_0$），从方法输出中提取预测的峰值位置 $\hat{\lambda}_0(\lambda_1)$。
2. 计算预测值与真实值的平方误差 $(\lambda^*_0(\lambda_1) - \hat{\lambda}_0(\lambda_1))^2$。
3. 对所有切片求平均，再开方得到PeakRMSE。

注意事项：

• 需要处理多峰情况。如果真实情况有 $k$ 个峰，则让方法返回其置信度最高的 $k$ 个峰进行匹配。
• 需要设置合理的峰值检测阈值，以避免将噪声误判为峰。
• 这个指标使得比较ClassiFIM、W方法、SPCA等不同范式的技术成为可能，只要它们都能在一维切片上输出一个可识别峰值的标量函数。

8. 计算资源与实现建议

本项目的计算资源主要消耗在基态制备上，特别是对于量子模型。

• CPU密集型：Hubbard12和FIL24的Lanczos计算消耗了约3000 CPU小时。XXZ300的DMRG计算消耗了约12000 CPU小时。这里的“CPU小时”指单逻辑核心运行1小时。实际计算中，我们使用了多核并行（每个参数点独立）来大幅缩短墙钟时间。
• GPU训练：ClassiFIM中神经网络模型的训练相对便宜，折算到单块NVIDIA GTX 1060 6GB GPU上，大约需要200 GPU小时。
• 实现语言：核心高性能计算部分（MCMC、Lanczos调用、DMRG）使用C++或依赖高度优化的库（如ARPACK， ITensor）。上层协调、数据分析和机器学习部分使用Python（NumPy, SciPy, PyTorch/JAX）。通过ctypes或pybind11实现Python与C++的交互。
• 可复现性：所有关键算法步骤（如MCMC平衡方案、DMRG收敛标准、FIM估计公式）都应详细记录并设置随机种子。对于耗时的计算，建议保存中间结果（如每个 $\lambda$ 点的概率分布或样本集），以便后续分析无需重新计算。

从样本估计费舍尔信息矩阵是一座连接统计推断、信息几何和物理相变的桥梁。无论是通过MCMC在经典随机场中摸索，还是借助Lanczos/DMRG窥探量子世界的基态，其核心思想都是绕过难以触及的完整概率分布，从有限的、嘈杂的样本或数值近似中，提取出系统对参数变化的根本敏感性。这个过程充满了数值的陷阱：MCMC中失衡的链、Lanczos中不收敛的迭代、有限差分中过大的步长、峰值检测中嘈杂的背景……每一个环节都需要细致的调优和交叉验证。

我个人的体会是，成功的关键在于对每一层近似所带来的误差保持清醒的认识。MCMC估计的FIM带有抽样误差和自相关误差；Lanczos给出的“精确”分布受限于数值精度和截断的能级；有限差分法本身就在用离散逼近连续。将这些误差量化、控制，并通过多次独立运行来评估统计不确定性，比追求单一运行的“完美”结果更重要。最终，当你在相变点附近看到那条陡然耸立的FIM曲线时，你会觉得所有这些繁琐的步骤都是值得的——因为你正在用数据本身的语言，聆听系统在临界点发出的最尖锐的声响。