Ising 模型及馬爾科夫鏈蒙特卡羅方法
這個月剛出爐的諾貝爾物理學獎讓大家重新認識了一個主題,復雜系統(tǒng)。喬治·帕里西(Giorgio Parisi)于上世紀 80 年代對自旋玻璃作了深入研究,發(fā)現(xiàn)了隨機現(xiàn)象背后的隱藏規(guī)則,被認為是對復雜系統(tǒng)理論最重要的貢獻之一。

自旋玻璃是一種特殊類型的金屬合金,其中鐵原子隨機混合進銅原子網(wǎng)格。每一個鐵原子的行為都像一個小磁鐵(或者說自旋),受到鄰近鐵原子的影響。在普通磁鐵中,所有自旋都指向同一方向,但在自旋玻璃中,一些自旋對想要指向同一方向,而另一些自旋對想要指向相反方向。那么,它們如何趨于平衡呢?
自旋玻璃的奇特性質為復雜系統(tǒng)提供了一個模型,在這些系統(tǒng)中,各個部分必須在各種反作用力間達到平衡。帕里西對自旋玻璃結構的基本發(fā)現(xiàn)非常深刻,使人們能夠理解和描述許多不同的、顯然完全隨機的材料和現(xiàn)象,不僅可用于物理學領域,而且在數(shù)學、生物學、神經科學和機器學習等領域也大顯身手。
帕里西還研究了許多其他現(xiàn)象,如:為什么冰河期會周期性地重復出現(xiàn)?對于混沌和湍流系統(tǒng),是否存在更具一般性的數(shù)學描述?數(shù)千只椋鳥的咕噥聲是如何形成特定模式的?這些問題看上去似乎與自旋玻璃相去甚遠。但帕里西說,他的大部分研究都關注簡單行為如何導致復雜的整體行為,無論對自旋玻璃還是對椋鳥群都同樣適用。
為了對所謂的復雜系統(tǒng)有所了解,我們來看一個這次拿下諾貝爾獎的自旋玻璃模型的簡化版本,即經典的 Ising 模型。

1Ising 模型簡史
Ising 模型的提出是為了解釋鐵磁物質的相變,即磁鐵在加熱到一定臨界溫度以上會出現(xiàn)磁性消失的現(xiàn)象,而降溫到臨界溫度以下又會表現(xiàn)出磁性。這種有磁性、無磁性兩相之間的轉變,是一種連續(xù)相變。
Ising 模型假設鐵磁物質是由一堆規(guī)則排列的小磁針構成,每個磁針只有上下兩個方向(自旋)。相鄰的小磁針之間通過能量約束發(fā)生相互作用,同時又會由于環(huán)境熱噪聲的干擾而發(fā)生磁性的隨機轉變。
漲落的大小由關鍵的溫度參數(shù)決定,溫度越高,隨機漲落干擾越強,小磁針越容易發(fā)生無序而劇烈地狀態(tài)轉變,從而讓上下兩個方向的磁性相互抵消,整個系統(tǒng)消失磁性。
如果溫度很低,則小磁針相對寧靜,系統(tǒng)處于能量約束高的狀態(tài),大量的小磁針方向一致,鐵磁系統(tǒng)展現(xiàn)出磁性。而當系統(tǒng)處于臨界溫度的時候,Ising 模型表現(xiàn)出一系列冪律行為和自相似現(xiàn)象。
由于 Ising 模型的高度抽象,人們可以很容易地將它應用到其他領域之中。
例如,人們將每個小磁針比喻為某個村落中的村民,而將小磁針上、下的兩種狀態(tài)比喻成個體所具備的兩種政治觀點(例如對 A,B 兩個不同候選人的選舉),相鄰小磁針之間的相互作用比喻成村民之間觀點的影響。環(huán)境的溫度比喻成每個村民對自己意見不堅持的程度。這樣,整個 Ising 模型就可以建模該村落中不同政治見解的動態(tài)演化(即觀點動力學 opinion dynamics)。在社會科學中,人們已經將 Ising 模型應用于股票市場、種族隔離、政治選擇等不同的問題。
另一方面,如果將小磁針比喻成神經元細胞,向上向下的狀態(tài)比喻成神經元的激活與抑制,小磁針的相互作用比喻成神經元之間的信號傳導,那么,Ising 模型的變種還可以用來建模神經網(wǎng)絡系統(tǒng),從而搭建可適應環(huán)境、不斷學習的機器,如 Hopfield 網(wǎng)絡和 Boltzmann 機之類。
Ising 模型不僅僅是一個統(tǒng)計物理模型,它更是一個建模各種復雜系統(tǒng)模型的典范。
2模型參數(shù)

假設第
表示磁針朝上或者朝下。網(wǎng)格上相鄰的兩個小磁針可以發(fā)生耦合作用。
+幾股勢力
用于統(tǒng)攝所有小磁針的一個磁場,表示整體的影響或約束,本文姑且稱它為正能量(這里的正不是 +,是以它的方向為正統(tǒng))。
相鄰小磁針的相互影響,表示局部的影響或約束,本文姑且稱為耦合能量。
外界的干擾勢力,比如外界溫度越高,個體小磁針越亢奮,正能量和耦合能量對它的統(tǒng)攝和約束能力就越小。
一個系統(tǒng)在這三股勢力的綜合影響下趨于平衡,那么怎么將這三股勢力統(tǒng)一用公式表達出來呢?
+總能量
我們先看前兩股勢力:
如果小磁針方向與正能量方向一致,則總能量降低;
如果兩個相鄰小磁針的狀態(tài)一致,則系統(tǒng)的總能量減
個單位,否則就增 個單位。
綜合這兩點,將總能量定義為,
其中,
為一個能量耦合常數(shù), 表示系統(tǒng)處于狀態(tài)組合 下的總能量。對于邊 ,如果 ,則總能量減少 。 參數(shù)
表示正能量,某個小磁針的方向如果與正能量一致,則總能量減少 。
我們看個例子,假設系統(tǒng)中只有

每個小磁針有
那么這
由能量項可知,小磁針相互之間以及與
沿用村民的比喻來說,系統(tǒng)的能量相當于村民觀點存在的沖突的數(shù)量。如果兩個相鄰的村民意見不一致,總沖突數(shù)就
+外部干擾
但是,系統(tǒng)的演化并不完全由總能量決定。由于小磁針暴露于外部環(huán)境中,熱漲落又會引起小磁針的狀態(tài)隨機反轉。可以用溫度
于是,每個小磁針就掙扎于這兩種不同的勢力之間:
假如外部干擾勢力溫度
趨于 ,則每個小磁針都會與外場相一致,那么,最終系統(tǒng)將處于全是 或者全是 的狀態(tài)。 假如
特別高,而 和 又特別小,則鄰居間的作用可以忽略,每個小磁針都將趨于隨機取值。
因此,整個 Ising 模型就有兩個參數(shù)
3統(tǒng)一建模
Ising 模型中并沒有針對每個小磁針設定其狀態(tài)轉變的規(guī)則,而是先讓每個小磁針的狀態(tài)發(fā)生隨機變化,再根據(jù)能量來依概率接受這種狀態(tài)變化。在 Ising 模型中,小磁針的自旋狀態(tài)分布取決于溫度,并遵循 Boltzmann 分布:
其中
簡寫為,
通過在不同溫度下多次運行 MCMC 分析來確定作為溫度函數(shù)的平均磁化強度。可以通過迭代每個可能的狀態(tài),計算給定磁化強度的能量,然后根據(jù)其頻率
+隨機選狀態(tài)
因此,與其遍歷所有可能的狀態(tài),不如隨機選擇一組狀態(tài)并按其頻率加權計算平均值。這樣做的問題是效率低下。玻爾茲曼分布并不均勻:它偏向于較低的能量。理想情況下,對于更頻繁出現(xiàn)的狀態(tài),我們希望用一種方法來更頻繁地找出它的自旋狀態(tài)。
4Metropolis-Hastings
而這正是核物理學家 Nicholas Metropolis 和他的合作者的貢獻。該方法最初發(fā)表于 1953 年。他們在曼哈頓項目的工作中率先引入了蒙特卡羅方法和馬爾可夫鏈蒙特卡羅方法。
他們使用他們新創(chuàng)建的技術在狀態(tài)空間中隨機游走,更頻繁地訪問更頻繁出現(xiàn)的自旋狀態(tài),并大大加快計算速度。用他們的話來說:
我們不是隨機選擇狀態(tài),然后用
對其進行加權; 而是以概率
隨機地選擇狀態(tài),然后均勻加權。
不過,他們并沒有解釋為什么這樣做可行。至于被命名為 Metropolis-Hastings 的原因,算法名字中加入 Hastings 意味著是承認統(tǒng)計學家 W. K. Hastings 在 1970 年的貢獻,他使得該技術背后的理論最終得到解釋和概括。從具有密度
計算依賴于
的地方僅僅是以 這樣的比率形式,其中 和 是樣本點。因此,不需要知道歸一化常數(shù),不需要對 進行因式分解,并且這些方法很容易在計算機上實現(xiàn)。 通過模擬馬爾可夫鏈獲得樣本序列。因此,產生的樣本是相關的,對標準差的估計和對誤差的估計可能需要比獨立樣本更加小心。
這個想法用到 Ising 模型上涉及到的物理學是自旋狀態(tài)遵循玻爾茲曼分布,即
具體取值如下,
:小磁針晶格的狀態(tài) :狀態(tài)的頻率分布 :無量綱溫度 : 所有相鄰自旋 的總和
用上面這些符號更新公式,
這里只能說它們之間是成正比的,而要像知道某個特定自旋狀態(tài)發(fā)生的概率,需要將它除以所有可能的狀態(tài),以便概率加起來為
但問題是我們并不知道分母,這也是最初困擾貝葉斯統(tǒng)計學家之處。幸運的是,我們可以使用 MCMC。正如 Hastings 所說,
計算僅僅以
這個比率形式依賴于 ,其中 和 是樣本點。因此,不需要知道歸一化常數(shù),不需要對 進行因式分解,并且這些方法很容易在計算機上實現(xiàn),
+具體實現(xiàn)
對于 Ising 模型,Metropolis-Hasting 算法遵循以下步驟:
1、選擇一個起始自旋狀態(tài)
2、選擇一個新的自旋狀態(tài)
3、根據(jù)新狀態(tài)與舊狀態(tài)的相對概率
4、重復步驟 2 和 3,直到收斂。實際上這里往往是設定一個
其中,關于
這意味著對于第 3 步,似然比是
這是一個比值,因此避免了對所有狀態(tài)求和
得到似然比為,
算法按照上述比率來選擇新舊狀態(tài)
5簡單實現(xiàn)
針對 Ising 模型的簡單模擬網(wǎng)上 Python 代碼很多,比如 https://rajeshrinet.github.io/blog/2014/ising-model/。
import?numpy?as?np
from?numpy.random?import?rand
import?matplotlib.pyplot?as?plt
%matplotlib?inline
#?Simulating?the?Ising?model
class?Ising():
????'''?Simulating?the?Ising?model?'''????
????
????##?monte?carlo?moves
????def?mcmove(self,?config,?N,?beta):
????????'''?This?is?to?execute?the?monte?carlo?moves?using?
????????Metropolis?algorithm?such?that?detailed
????????balance?condition?is?satisified'''
????????for?i?in?range(N):
????????????for?j?in?range(N):????????????
????????????????????a?=?np.random.randint(0,?N)
????????????????????b?=?np.random.randint(0,?N)
????????????????????s?=??config[a,?b]
????????????????????nb?=?config[(a+1)%N,b]?+?config[a,(b+1)%N]?+?config[(a-1)%N,b]?+?config[a,(b-1)%N]
????????????????????cost?=?2??s??nb
????????????????????if?cost?0:
????????????????????????s?*=?-1
????????????????????elif?rand()?????????????????????????s?*=?-1
????????????????????config[a,?b]?=?s
????????return?config
????
????def?simulate(self):???
????????'''?This?module?simulates?the?Ising?model'''
????????N,?temp?????=?64,?.4????????#?Initialse?the?lattice
????????config?=?2*np.random.randint(2,?size=(N,N))-1
????????f?=?plt.figure(figsize=(15,?15),?dpi=80);????
????????self.configPlot(f,?config,?0,?N,?1);
????????
????????msrmnt?=?1001
????????for?i?in?range(msrmnt):
????????????self.mcmove(config,?N,?1.0/temp)
????????????if?i?==?1:???????self.configPlot(f,?config,?i,?N,?2);
????????????if?i?==?4:???????self.configPlot(f,?config,?i,?N,?3);
????????????if?i?==?32:??????self.configPlot(f,?config,?i,?N,?4);
????????????if?i?==?128:?????self.configPlot(f,?config,?i,?N,?5);
????????????if?i?==?256:????self.configPlot(f,?config,?i,?N,?6);
??????????????????
????def?configPlot(self,?f,?config,?i,?N,?n_):
????????'''?This?modules?plts?the?configuration?once?passed?to?it?along?with?time?etc?'''
????????X,?Y?=?np.meshgrid(range(N),?range(N))
????????sp?=??f.add_subplot(3,?3,?n_?)??
????????plt.setp(sp.get_yticklabels(),?visible=False)
????????plt.setp(sp.get_xticklabels(),?visible=False)??????
????????plt.pcolormesh(X,?Y,?config,?cmap=plt.cm.RdBu);
????????plt.title('Time=%d'%i);?plt.axis('tight')????
????plt.show()
rm?=?Ising()
rm.simulate()
運行結果如下,

