使用王朗道蒙特卡罗算法和机器学习研究伊辛模型相变

曹红英，田华*，侯杰，张万舟

（1.太原理工大学 信息与计算机学院，山西 晋中 030600；2.太原理工大学 物理学院，山西 晋中 030600）

0 引言

研究物质的相和相变在统计物理领域中是一个传统且重要的任务，目前仍然是一个十分活跃的研究方向［1］。近年来，随着机器学习的发展，人们将机器学习引入到统计物理领域用于识别相变［2］。开始人们使用机器学习中的监督学习结合经典的模型来对相变进行研究，包括伊辛模型［3］、具有拓扑相变的 XY 模型［4］、Potts模型［5］等。后来，研究发现虽然监督学习可以得到比较好的结果，但是需要提前定义数据标签，不适于探索未知物理系统的相变。相较于监督学习，非监督学习不需要提前设置数据标签［6］也能获得系统的相变点。例如：通过扩散映射可以对XY 模型的相变进行研究［7-8］。

对于机器学习来说，数据集是否完备直接决定了结果的准确性。在研究物质相变时，数据集是均匀覆盖整个相空间的构型集合，这些构型由蒙特卡罗方法获得。最常用的蒙特卡罗方法是Metropolis 算法［9］，但该算法存在不足之处：低温下自旋演化概率低［10］，存在临界慢化问题，在一级相变中存在亚稳态问题［11］，这些不足导致数据集制备不完整。

低温下自旋难以演化是指，使用 Metropolis算法模拟物理模型时，自旋的演化概率正比于exp(−ΔE T)，其中T是温度，ΔE是自旋翻转前后系统的能量差。低温极限下，演化概率很小，系统很难遍历到一些可能的典型状态。如模拟伊辛模型时［12］，状态可能会局限在两种可能的铁磁态中的一种，只生成自旋全向上或自旋全向下状态，而真实情况下的状态为两者共存，经过数据降维后，反映在二维空间散点图中是一个只有横纵坐标的数据点，无法运用聚类算法研究和分析相变。

临界慢化是指在相变温度附近，Metropolis方法需要更新足够长的步数才能获得独立的或关联较弱的构型。亚稳态中的“亚”是状态对应的并不是真正的全局能量最低，但是接近最低能量，“稳”是指状态对应的能量是局部能量最小，达到局部稳定的状态。亚稳态是局部最小能量对应的状态，并不是全局最小能量对应的状态，系统一旦进入亚稳态，就很难跳过能量壁垒演化到其他状态，即一旦处于亚稳态，Metropolis方法很难帮助系统遍历到所有态。

为了克服上述不足导致的构型数据集不完备的情况，本文采用王朗道算法（Wang-Landau algorithm）［13］产生模型的构型。与 Metropolis 算法不同的是，该算法直接用自旋翻转前后的能量对应的态密度之比作为演化概率，其中Eold表示翻转之前的系统能量，Enew表示翻转之后的系统能量，该演化概率与温度无关，低温下可以遍历到所有构型。王朗道算法也不存在临界慢化，同时可以避免系统被束缚在亚稳态。

本文使用王朗道算法生成模型的构型数据集后，使用非监督学习的方法对数据集进行降维和聚类，最后使用ch指标（Calinski-Harabasz Index）［14］准确得到了含磁场伊辛模型的一级相变点和不含磁场伊辛模型的二级相变点。

1 相关方法

本文采用王朗道蒙特卡罗方法和机器学习中的非监督学习方法研究经典模型的相变。首先采用王朗道算法生成模型的构型，并将这些构型状态作为数据集；之后使用非监督学习方法的主成分分析法［15］、K-means算法［16］对数据集进行处理；最后通过聚类评估指标ch找出模型的相变点。算法流程图如图1所示。

1.1 王朗道蒙特卡罗算法

与Metropolis算法对样本的抽样方式不同，王朗道算法通过随机翻转自旋状态在能量空间随机行走，与自旋组态相关的能量E以与其态密度g(E)的倒数成比例的概率被接受，模拟过程中更新记录状态密度g(E)。随机行走中，为了遍历到所有能量空间，引入累计直方图H(E)，能量E每被访问一次，H(E)对应值增加1，最终得到每个能量E的访问次数及访问分布情况。得到平坦的能量直方图，意味着每个能量都被均匀访问。按照能量由小到大的顺序给出能级编号，记录下能级对应的多种构型就可以制备出包含完整构型的数据集。

王朗道算法的具体步骤为：（1）设置整个晶格的初始自旋状态，为所有可能的能量E设置态密度g(E)=1，设置修正因子，如：f=2。（2）随机翻转一个自旋态形成试验态，如果自旋翻转前后的Eold和Enew能量发生变化，则从能量Eold到Enew的跃迁概率为，

如果g(Enew)≤g(Eold)， 则翻转后的状态就被接受，否则以的概率被随机接受。如果翻转后的状态被接受，则更新g(Enew)→g(Enew)×f，更新能量直方图中的H(Enew)，即 H(Enew)→H(Enew)+1；如果随机行走的状态被拒绝接受并保持相同的能量Eold，就 乘以 存在 的态 密度g(Eold)，即g(Eold

)→g(Eold)×f，更 新 能 量 直 方 图 ，即H(Eold)→H(Eold)+1。因为g(E)会变得很大，所以在算法中采用的是态密度的对数，即ln[g(E)]→ln[g(E)]×f。（3）继 续 执 行 步 骤2，直到直方图平坦，使用函数减小修正系数。（4）重复步骤 2和步骤 3，直到 f=fmin，得到最后一次的态密度g(E)。

随机行走过程中，记录每个能量对应的不同构型。根据不同能量Ei对应的态密度求出能量权重，按照能级对应的能量由小至大的顺序对权重进行归一化处理，即，其中 n为能级的个数。之后根据权重生成概率区间

通过不同区间的能量差对应的能级，选出能级对应的构型制备数据集。例如：随机取一个［0，1］之 间 的 随 机 数 γ，如 果 γ 落 在[P(E1)，P(E1)+P(E2)]这一区间内，则任意选取能量E2对应的能级构型中的一种记录下来，重复上述取样的过程，制备出数据集。

根据态密度，内能U（T）可以通过下式计算，

其中，E是能量，g（E）是态密度，T是温度，kB是玻尔兹曼常数。比热C（T）可以由内能的涨落来确定：

绝对磁化强度M（T）：

1.2 非监督学习相关方法

使用王朗道方法得到模型尺寸L的构型数据集X，可表示为

其中，数据集的每一行有n个数据，代表L×L=n个格点自旋状态组成的构型数据，m表示采样的构型数有m行。

使用主成分分析法对数据进行降维。首先将数据集中的每一行去中心化后计算其协方差矩阵，通过奇异值分解求出协方差矩阵的最大两个特征值对应的特征向量pca1和pca2，组成二维空间，最后将原始数据映射到该二维空间中。

使用K-means方法对降维后的数据进行聚类，在不给定任何标签的情况下，让数据自行根据特征聚合在一起。由于聚类后，簇内数据差异性小，簇间数据差异性大［17］。ch指标是通过计算簇间平方误差和与簇内平方误差和的比值获得的。簇内数据距离越小，簇间数据距离越大，则ch值越大。ch由以下公式给出

其中，N为数据集样本数，k为样本簇数，BGSS（Between-Groups Sum of Squared Error）是簇间平方误差和，用来衡量簇间的分离度，BGSS越大，cht越大，簇越分散。WGSS（Within-Groups Sum of Squared Error）是簇内平方误差和，用来度量簇内的紧密度，WGSS越大，chb越大，簇内样本越分散。研究模型的相和相变时，同一相内数据差异小，不同相之间数据差异大，则可以通过计算出每个温度下的ch及其分量cht和chb的值，根据簇内和簇间的分散程度寻找模型的相变点。

2 分析和讨论

现在通过分析二维伊辛模型的相变，测试方法的可行性。伊辛模型［18］的哈密顿量为：

2.1 Metropolis算法和王朗道算法的对比分析

为了证明王朗道算法的正确性，运用该算法对尺寸L=16的二维伊辛模型的比热和内能进行计算，将结果与通过Metropolis算法得出的结果进行对比。如图 2（a）和图 2（b）所示，二者结果一致。图2说明了本文所运用的王朗道算法以及代码是正确的。

但是使用机器学习研究相变时，需要得到同一能量下对应的多种构型来构建更完备的数据集。为了对比Metropolis算法和王朗道算法产生的构型数据集，本文采用两种算法对低温极限T=0.1下，尺寸大小L=16的二维伊辛模型进行构型采样。设置构型的初始状态为自旋全部向上，用两种算法各产生10 000个构型。将这些构型组成的数据集表示为10 000×256的矩阵，每一行是16×16=256个格点自旋状态组成的数据。其中，自旋向上和向下分别用1，−1表示。使用1.2小节的主成分分析法和K-means算法处理数据画出散点图。

图3是用两种算法得出的散点图。通过Metropolis算法，在低温极限下，自旋很难翻转。所有的自旋均处于设置的初始构型，即自旋全部向上。则在构型数据集中，矩阵元素全都是1，得出的协方差矩阵的最大的两个特征向量都是零向量，将原始矩阵元素映射到这两个零向量组成的二维空间中，只得到一个处于（0，0）点的数据，反映在图3（a）中只有一个位于（0，0）点的团簇。王朗道算法在获取构型时不受温度影响，在低温极限下可获得自旋向上和自旋向下的的构型，这是伊辛模型低温下完备的构型类型，反映在图3（b）中是两个相互分离的团簇。

王朗道算法和Metropolis算法均可产生构型数据集，但相较于Metropolis算法，王朗道算法不受低温的影响，可以采样到所有可能的重要构型来制备数据集，更加适用于非监督学习研究物理模型的相变。

2.2 无磁场伊辛模型的二级相变

通过王朗道算法得到伊辛模型的数据集时，需要求出态密度g（E）。设置整个晶格的初始自旋状态为全部向上，为所有可能的能量E设置态密度g（E）=1，初始修正因子设置为f =f0=2，通过开平方减小修正因子f，当修正因子小于 min_ f=exp（10−8） ≈ 1.000 000 01 时直方图平坦，停止模拟。实验环境为Windows 10操作系统，Intel E5-2680 CPU，64 GB内存，11 GB显存，经过5 452.02 s计算得到尺寸为L=16的二维伊辛模型的态密度，如图4所示。

使用王朗道算法得到伊辛模型的数据集之后，再使用主成分分析法和K-means算法进行数据处理，接着使用ch指标进行相变分析，计算出每个温度下的ch分量cht和chb，最后使用“0-1归一化”（将某一温度下的cht或chb的值与全部温度下该值的最小值的差值和该值最大值、最小值的差值作比，即x_norm=(x−x_min)/(x_max−x_min)，其 中 ，x_max，x_min分别代表x中的最大值和最小值，得到［0，1］之间无量纲的相对值）。如图5（a）—（b）所示，ch指标的分量chb在T =2.3左右处于尖峰位置，ch指标的分量cht在T =2.3左右有显著变化，与图6（a）—（b）所示的同尺寸下测得的比热尖峰C和平均磁化强度M得到的相变点的结果一致。chb和cht指示的相变点结果与热力学极限（系统尺寸无穷大时）理论相变点2.269存在一定的误差，误差与尺寸效应和定位精度相关。

为了进一步理解ch指标的结果，图7（a）−（d）给出温度 T=0.1、2.0、2.3、5.0时降维聚类后的散点图。可以观察到，低温下，由两种颜色标识的簇相互分离且簇内聚合；随着温度的升高，簇相互靠近，簇内分散；最后簇在高温下混合在一起，簇内聚合。

由1.2小节对ch指标的介绍可知，簇内越分散，chb值越大。在低温阶段，伊辛模型整体处于有序状态，只有自旋向上和自旋向下两种构型状态，构成了散点图的两个团簇，簇内点聚合，chb值小；在相变点T=2.3附近，伊辛模型处于有序和无序共存的状态，样本涨落最大，簇内点分散程度最大，chb值最大；在高温下自旋状态杂乱无章，整体处于无序的状态，样本涨落小，簇内的点聚合，chb值随着远离相变点而逐步减小。类似地，簇间越分散，cht值越大。在低温下，簇间分散程度最大，cht值最大。温度升高，簇间距离逐渐缩小，在相变点T=2.3附近，cht值有显著变化；在高温下，两簇逐渐融合，cht值变小。结果表明，ch指标的分量cht和chb在检测伊辛模型的相变过程中表现良好。

经过上述分析，通过王朗道算法制备出伊辛模型的数据集后，经过非监督学习，可以确定伊辛模型二级相变时的临界点，结果与理论相变点吻合。

2.3 含磁场伊辛模型的一级相变

2.2 小节的研究表明本文的方法对于研究伊辛模型的二级相变是可行的。Metropolis算法在模拟含磁场伊辛模型时，系统很容易进入亚稳态，不利于机器学习。接下来讨论本文提出的方法能否克服这一困难，得到含磁场伊辛模型的一级相变点。

磁场作用下的二维伊辛模型的哈密顿量为

其中，h为磁场强度，N是总格点数。相应的配分函数计算公式如下

其中，E′是能量，M′是磁化强度，h是磁场强度。内能U（T）可以通过下式计算

绝对磁化强度

图8是经过公式（10）计算的尺寸L=16磁场作用下的伊辛模型在不同温度下的平均磁化强度M（T，h）。考虑磁场的情况下，温度T=1.0时，伊辛模型表现出强一级相变。本节研究在T=1.0时含磁场的伊辛模型。

如图 10 （a）—（b）是在温度 T=1.0时磁场作用下的伊辛模型采用王朗道算法得到构型数据集，使用非监督学习方法得到的结果图，可以发现ch指标的分量cht和chb在h = 0.0左右处于尖峰位置。与图10（c）所示的相同尺寸下L=16的测量的比热尖峰的结果和图9所示的磁化强度跳变的对应磁场大小一致。此处计算出的cht和chb值，同2.2小节一样使用了“0-1归一化”得到了在［0，1］之间无量纲的相对值。

为了进一步理解上述ch指标得到的结果，图 11（a）—（c） 给出在磁场强度 h = −1.0、0.0、1.0降维聚类后的散点图。可以观察到，在h = −1.0时，两簇聚合；在h = 0.0相变点处两簇分离，而且团簇比较大；在h = 1.0时，两簇聚合。

簇越分散，cht值越大。在h = 0.0相变点处，自旋全向上和自旋全向下两种构型，对应能量相等，整体处于二重简并的状态，两种构型对应的两个团簇分离程度最大，cht值最大；在其他负磁场和正磁场下，自旋被磁化，和磁场保持一致，数据集再次为全部相同的自旋构型，在散点图中表现为聚合在（0， 0）附近的团簇，cht值就小。cht值在检测伊辛模型相变过程中表现良好。类似地，簇内越分散，chb值越大。在h = 0.0相变点处团簇较大，chb值最大。结果表明，ch指标的分量cht和chb在检测磁场作用下伊辛模型的相变过程中表现良好。

经过上述分析，通过王朗道算法可以避免伊辛模型在发生一级相变时局限在亚稳态，从而制备出磁场作用下二维伊辛模型完备的构型数据集，之后经过非监督学习的处理，可以确定含磁场伊辛模型一级相变时的临界点，且结果与理论相变点吻合。

3 结论

本文运用伊辛模型作为研究对象，采用王朗道算法构建出完备的数据集，使用主成分分析法和K-means算法对数据处理后，结合非监督学习的ch指标分量cht和chb表征出模型在相变时的突变情况，确定了二维伊辛模型相变点的位置。研究结果表明采用王朗道算法可以得出较为完整的数据集，有助于机器学习方法研究相变，并可以准确得到相变点。本文为机器学习方法研究相变提供了一种新思路，促进机器学习方法在统计物理领域中的应用。