使用王朗道蒙特卡羅算法和機器學習研究伊辛模型相變

曹紅英，田華*，侯杰，張萬舟

（1.太原理工大學 信息與計算機學院，山西 晉中 030600；2.太原理工大學 物理學院，山西 晉中 030600）

0 引言

研究物質的相和相變在統(tǒng)計物理領域中是一個傳統(tǒng)且重要的任務，目前仍然是一個十分活躍的研究方向［1］。近年來，隨著機器學習的發(fā)展，人們將機器學習引入到統(tǒng)計物理領域用于識別相變［2］。開始人們使用機器學習中的監(jiān)督學習結合經(jīng)典的模型來對相變進行研究，包括伊辛模型［3］、具有拓撲相變的 XY 模型［4］、Potts模型［5］等。后來，研究發(fā)現(xiàn)雖然監(jiān)督學習可以得到比較好的結果，但是需要提前定義數(shù)據(jù)標簽，不適于探索未知物理系統(tǒng)的相變。相較于監(jiān)督學習，非監(jiān)督學習不需要提前設置數(shù)據(jù)標簽［6］也能獲得系統(tǒng)的相變點。例如：通過擴散映射可以對XY 模型的相變進行研究［7-8］。

對于機器學習來說，數(shù)據(jù)集是否完備直接決定了結果的準確性。在研究物質相變時，數(shù)據(jù)集是均勻覆蓋整個相空間的構型集合，這些構型由蒙特卡羅方法獲得。最常用的蒙特卡羅方法是Metropolis 算法［9］，但該算法存在不足之處：低溫下自旋演化概率低［10］，存在臨界慢化問題，在一級相變中存在亞穩(wěn)態(tài)問題［11］，這些不足導致數(shù)據(jù)集制備不完整。

低溫下自旋難以演化是指，使用 Metropolis算法模擬物理模型時，自旋的演化概率正比于exp(?ΔE T)，其中T是溫度，ΔE是自旋翻轉前后系統(tǒng)的能量差。低溫極限下，演化概率很小，系統(tǒng)很難遍歷到一些可能的典型狀態(tài)。如模擬伊辛模型時［12］，狀態(tài)可能會局限在兩種可能的鐵磁態(tài)中的一種，只生成自旋全向上或自旋全向下狀態(tài)，而真實情況下的狀態(tài)為兩者共存，經(jīng)過數(shù)據(jù)降維后，反映在二維空間散點圖中是一個只有橫縱坐標的數(shù)據(jù)點，無法運用聚類算法研究和分析相變。

臨界慢化是指在相變溫度附近，Metropolis方法需要更新足夠長的步數(shù)才能獲得獨立的或關聯(lián)較弱的構型。亞穩(wěn)態(tài)中的“亞”是狀態(tài)對應的并不是真正的全局能量最低，但是接近最低能量，“穩(wěn)”是指狀態(tài)對應的能量是局部能量最小，達到局部穩(wěn)定的狀態(tài)。亞穩(wěn)態(tài)是局部最小能量對應的狀態(tài)，并不是全局最小能量對應的狀態(tài)，系統(tǒng)一旦進入亞穩(wěn)態(tài)，就很難跳過能量壁壘演化到其他狀態(tài)，即一旦處于亞穩(wěn)態(tài)，Metropolis方法很難幫助系統(tǒng)遍歷到所有態(tài)。

為了克服上述不足導致的構型數(shù)據(jù)集不完備的情況，本文采用王朗道算法（Wang-Landau algorithm）［13］產(chǎn)生模型的構型。與 Metropolis 算法不同的是，該算法直接用自旋翻轉前后的能量對應的態(tài)密度之比作為演化概率，其中Eold表示翻轉之前的系統(tǒng)能量，Enew表示翻轉之后的系統(tǒng)能量，該演化概率與溫度無關，低溫下可以遍歷到所有構型。王朗道算法也不存在臨界慢化，同時可以避免系統(tǒng)被束縛在亞穩(wěn)態(tài)。

本文使用王朗道算法生成模型的構型數(shù)據(jù)集后，使用非監(jiān)督學習的方法對數(shù)據(jù)集進行降維和聚類，最后使用ch指標（Calinski-Harabasz Index）［14］準確得到了含磁場伊辛模型的一級相變點和不含磁場伊辛模型的二級相變點。

1 相關方法

本文采用王朗道蒙特卡羅方法和機器學習中的非監(jiān)督學習方法研究經(jīng)典模型的相變。首先采用王朗道算法生成模型的構型，并將這些構型狀態(tài)作為數(shù)據(jù)集；之后使用非監(jiān)督學習方法的主成分分析法［15］、K-means算法［16］對數(shù)據(jù)集進行處理；最后通過聚類評估指標ch找出模型的相變點。算法流程圖如圖1所示。

1.1 王朗道蒙特卡羅算法

與Metropolis算法對樣本的抽樣方式不同，王朗道算法通過隨機翻轉自旋狀態(tài)在能量空間隨機行走，與自旋組態(tài)相關的能量E以與其態(tài)密度g(E)的倒數(shù)成比例的概率被接受，模擬過程中更新記錄狀態(tài)密度g(E)。隨機行走中，為了遍歷到所有能量空間，引入累計直方圖H(E)，能量E每被訪問一次，H(E)對應值增加1，最終得到每個能量E的訪問次數(shù)及訪問分布情況。得到平坦的能量直方圖，意味著每個能量都被均勻訪問。按照能量由小到大的順序給出能級編號，記錄下能級對應的多種構型就可以制備出包含完整構型的數(shù)據(jù)集。

王朗道算法的具體步驟為：（1）設置整個晶格的初始自旋狀態(tài)，為所有可能的能量E設置態(tài)密度g(E)=1，設置修正因子，如：f=2。（2）隨機翻轉一個自旋態(tài)形成試驗態(tài)，如果自旋翻轉前后的Eold和Enew能量發(fā)生變化，則從能量Eold到Enew的躍遷概率為，

如果g(Enew)≤g(Eold)， 則翻轉后的狀態(tài)就被接受，否則以的概率被隨機接受。如果翻轉后的狀態(tài)被接受，則更新g(Enew)→g(Enew)×f，更新能量直方圖中的H(Enew)，即 H(Enew)→H(Enew)+1；如果隨機行走的狀態(tài)被拒絕接受并保持相同的能量Eold，就 乘以 存在 的態(tài) 密度g(Eold)，即g(Eold

)→g(Eold)×f，更 新 能 量 直 方 圖 ，即H(Eold)→H(Eold)+1。因為g(E)會變得很大，所以在算法中采用的是態(tài)密度的對數(shù)，即ln[g(E)]→ln[g(E)]×f。（3）繼 續(xù) 執(zhí) 行 步 驟2，直到直方圖平坦，使用函數(shù)減小修正系數(shù)。（4）重復步驟 2和步驟 3，直到 f=fmin，得到最后一次的態(tài)密度g(E)。

隨機行走過程中，記錄每個能量對應的不同構型。根據(jù)不同能量Ei對應的態(tài)密度求出能量權重，按照能級對應的能量由小至大的順序對權重進行歸一化處理，即，其中 n為能級的個數(shù)。之后根據(jù)權重生成概率區(qū)間

通過不同區(qū)間的能量差對應的能級，選出能級對應的構型制備數(shù)據(jù)集。例如：隨機取一個［0，1］之 間 的 隨 機 數(shù) γ，如 果 γ 落 在[P(E1)，P(E1)+P(E2)]這一區(qū)間內，則任意選取能量E2對應的能級構型中的一種記錄下來，重復上述取樣的過程，制備出數(shù)據(jù)集。

根據(jù)態(tài)密度，內能U（T）可以通過下式計算，

其中，E是能量，g（E）是態(tài)密度，T是溫度，kB是玻爾茲曼常數(shù)。比熱C（T）可以由內能的漲落來確定：

絕對磁化強度M（T）：

1.2 非監(jiān)督學習相關方法

使用王朗道方法得到模型尺寸L的構型數(shù)據(jù)集X，可表示為

其中，數(shù)據(jù)集的每一行有n個數(shù)據(jù)，代表L×L=n個格點自旋狀態(tài)組成的構型數(shù)據(jù)，m表示采樣的構型數(shù)有m行。

使用主成分分析法對數(shù)據(jù)進行降維。首先將數(shù)據(jù)集中的每一行去中心化后計算其協(xié)方差矩陣，通過奇異值分解求出協(xié)方差矩陣的最大兩個特征值對應的特征向量pca1和pca2，組成二維空間，最后將原始數(shù)據(jù)映射到該二維空間中。

使用K-means方法對降維后的數(shù)據(jù)進行聚類，在不給定任何標簽的情況下，讓數(shù)據(jù)自行根據(jù)特征聚合在一起。由于聚類后，簇內數(shù)據(jù)差異性小，簇間數(shù)據(jù)差異性大［17］。ch指標是通過計算簇間平方誤差和與簇內平方誤差和的比值獲得的。簇內數(shù)據(jù)距離越小，簇間數(shù)據(jù)距離越大，則ch值越大。ch由以下公式給出

其中，N為數(shù)據(jù)集樣本數(shù)，k為樣本簇數(shù)，BGSS（Between-Groups Sum of Squared Error）是簇間平方誤差和，用來衡量簇間的分離度，BGSS越大，cht越大，簇越分散。WGSS（Within-Groups Sum of Squared Error）是簇內平方誤差和，用來度量簇內的緊密度，WGSS越大，chb越大，簇內樣本越分散。研究模型的相和相變時，同一相內數(shù)據(jù)差異小，不同相之間數(shù)據(jù)差異大，則可以通過計算出每個溫度下的ch及其分量cht和chb的值，根據(jù)簇內和簇間的分散程度尋找模型的相變點。

2 分析和討論

現(xiàn)在通過分析二維伊辛模型的相變，測試方法的可行性。伊辛模型［18］的哈密頓量為：

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

為了證明王朗道算法的正確性，運用該算法對尺寸L=16的二維伊辛模型的比熱和內能進行計算，將結果與通過Metropolis算法得出的結果進行對比。如圖 2（a）和圖 2（b）所示，二者結果一致。圖2說明了本文所運用的王朗道算法以及代碼是正確的。

但是使用機器學習研究相變時，需要得到同一能量下對應的多種構型來構建更完備的數(shù)據(jù)集。為了對比Metropolis算法和王朗道算法產(chǎn)生的構型數(shù)據(jù)集，本文采用兩種算法對低溫極限T=0.1下，尺寸大小L=16的二維伊辛模型進行構型采樣。設置構型的初始狀態(tài)為自旋全部向上，用兩種算法各產(chǎn)生10 000個構型。將這些構型組成的數(shù)據(jù)集表示為10 000×256的矩陣，每一行是16×16=256個格點自旋狀態(tài)組成的數(shù)據(jù)。其中，自旋向上和向下分別用1，?1表示。使用1.2小節(jié)的主成分分析法和K-means算法處理數(shù)據(jù)畫出散點圖。

圖3是用兩種算法得出的散點圖。通過Metropolis算法，在低溫極限下，自旋很難翻轉。所有的自旋均處于設置的初始構型，即自旋全部向上。則在構型數(shù)據(jù)集中，矩陣元素全都是1，得出的協(xié)方差矩陣的最大的兩個特征向量都是零向量，將原始矩陣元素映射到這兩個零向量組成的二維空間中，只得到一個處于（0，0）點的數(shù)據(jù)，反映在圖3（a）中只有一個位于（0，0）點的團簇。王朗道算法在獲取構型時不受溫度影響，在低溫極限下可獲得自旋向上和自旋向下的的構型，這是伊辛模型低溫下完備的構型類型，反映在圖3（b）中是兩個相互分離的團簇。

王朗道算法和Metropolis算法均可產(chǎn)生構型數(shù)據(jù)集，但相較于Metropolis算法，王朗道算法不受低溫的影響，可以采樣到所有可能的重要構型來制備數(shù)據(jù)集，更加適用于非監(jiān)督學習研究物理模型的相變。

2.2 無磁場伊辛模型的二級相變

通過王朗道算法得到伊辛模型的數(shù)據(jù)集時，需要求出態(tài)密度g（E）。設置整個晶格的初始自旋狀態(tài)為全部向上，為所有可能的能量E設置態(tài)密度g（E）=1，初始修正因子設置為f =f0=2，通過開平方減小修正因子f，當修正因子小于 min_ f=exp（10?8） ≈ 1.000 000 01 時直方圖平坦，停止模擬。實驗環(huán)境為Windows 10操作系統(tǒng)，Intel E5-2680 CPU，64 GB內存，11 GB顯存，經(jīng)過5 452.02 s計算得到尺寸為L=16的二維伊辛模型的態(tài)密度，如圖4所示。

使用王朗道算法得到伊辛模型的數(shù)據(jù)集之后，再使用主成分分析法和K-means算法進行數(shù)據(jù)處理，接著使用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指示的相變點結果與熱力學極限（系統(tǒng)尺寸無窮大時）理論相變點2.269存在一定的誤差，誤差與尺寸效應和定位精度相關。

為了進一步理解ch指標的結果，圖7（a）?（d）給出溫度 T=0.1、2.0、2.3、5.0時降維聚類后的散點圖?？梢杂^察到，低溫下，由兩種顏色標識的簇相互分離且簇內聚合；隨著溫度的升高，簇相互靠近，簇內分散；最后簇在高溫下混合在一起，簇內聚合。

由1.2小節(jié)對ch指標的介紹可知，簇內越分散，chb值越大。在低溫階段，伊辛模型整體處于有序狀態(tài)，只有自旋向上和自旋向下兩種構型狀態(tài)，構成了散點圖的兩個團簇，簇內點聚合，chb值??；在相變點T=2.3附近，伊辛模型處于有序和無序共存的狀態(tài)，樣本漲落最大，簇內點分散程度最大，chb值最大；在高溫下自旋狀態(tài)雜亂無章，整體處于無序的狀態(tài)，樣本漲落小，簇內的點聚合，chb值隨著遠離相變點而逐步減小。類似地，簇間越分散，cht值越大。在低溫下，簇間分散程度最大，cht值最大。溫度升高，簇間距離逐漸縮小，在相變點T=2.3附近，cht值有顯著變化；在高溫下，兩簇逐漸融合，cht值變小。結果表明，ch指標的分量cht和chb在檢測伊辛模型的相變過程中表現(xiàn)良好。

經(jīng)過上述分析，通過王朗道算法制備出伊辛模型的數(shù)據(jù)集后，經(jīng)過非監(jiān)督學習，可以確定伊辛模型二級相變時的臨界點，結果與理論相變點吻合。

2.3 含磁場伊辛模型的一級相變

2.2 小節(jié)的研究表明本文的方法對于研究伊辛模型的二級相變是可行的。Metropolis算法在模擬含磁場伊辛模型時，系統(tǒng)很容易進入亞穩(wěn)態(tài)，不利于機器學習。接下來討論本文提出的方法能否克服這一困難，得到含磁場伊辛模型的一級相變點。

磁場作用下的二維伊辛模型的哈密頓量為

其中，h為磁場強度，N是總格點數(shù)。相應的配分函數(shù)計算公式如下

其中，E′是能量，M′是磁化強度，h是磁場強度。內能U（T）可以通過下式計算

絕對磁化強度

圖8是經(jīng)過公式（10）計算的尺寸L=16磁場作用下的伊辛模型在不同溫度下的平均磁化強度M（T，h）?？紤]磁場的情況下，溫度T=1.0時，伊辛模型表現(xiàn)出強一級相變。本節(jié)研究在T=1.0時含磁場的伊辛模型。

如圖 10 （a）—（b）是在溫度 T=1.0時磁場作用下的伊辛模型采用王朗道算法得到構型數(shù)據(jù)集，使用非監(jiān)督學習方法得到的結果圖，可以發(fā)現(xiàn)ch指標的分量cht和chb在h = 0.0左右處于尖峰位置。與圖10（c）所示的相同尺寸下L=16的測量的比熱尖峰的結果和圖9所示的磁化強度跳變的對應磁場大小一致。此處計算出的cht和chb值，同2.2小節(jié)一樣使用了“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相變點處，自旋全向上和自旋全向下兩種構型，對應能量相等，整體處于二重簡并的狀態(tài)，兩種構型對應的兩個團簇分離程度最大，cht值最大；在其他負磁場和正磁場下，自旋被磁化，和磁場保持一致，數(shù)據(jù)集再次為全部相同的自旋構型，在散點圖中表現(xiàn)為聚合在（0， 0）附近的團簇，cht值就小。cht值在檢測伊辛模型相變過程中表現(xiàn)良好。類似地，簇內越分散，chb值越大。在h = 0.0相變點處團簇較大，chb值最大。結果表明，ch指標的分量cht和chb在檢測磁場作用下伊辛模型的相變過程中表現(xiàn)良好。

經(jīng)過上述分析，通過王朗道算法可以避免伊辛模型在發(fā)生一級相變時局限在亞穩(wěn)態(tài)，從而制備出磁場作用下二維伊辛模型完備的構型數(shù)據(jù)集，之后經(jīng)過非監(jiān)督學習的處理，可以確定含磁場伊辛模型一級相變時的臨界點，且結果與理論相變點吻合。

3 結論

本文運用伊辛模型作為研究對象，采用王朗道算法構建出完備的數(shù)據(jù)集，使用主成分分析法和K-means算法對數(shù)據(jù)處理后，結合非監(jiān)督學習的ch指標分量cht和chb表征出模型在相變時的突變情況，確定了二維伊辛模型相變點的位置。研究結果表明采用王朗道算法可以得出較為完整的數(shù)據(jù)集，有助于機器學習方法研究相變，并可以準確得到相變點。本文為機器學習方法研究相變提供了一種新思路，促進機器學習方法在統(tǒng)計物理領域中的應用。