物理雜記3 - Ising模型與模擬退火

Ising模型主要是用來描述物質的鐵磁性 (ferromagnetism)，是由德國物理學家Ernst Ising在1920左右提出了處理一維的解析解。該模型包含了每個晶格上原子的自旋 (spin) 參數，兩兩相鄰的原子彼此間可以是偏好處在共同的態上 (都是自旋向上)  或是處在相反的態上 (一上一下)，而外部的磁場會對該晶格上的自旋狀態施加能量，導致其偏向某一個方向。在這篇文章裡，我們不重複Ising當年複雜的數學解析解，而且該解也只能適用在1D的狀態，對於2D以上的系統目前並沒有解析解，所以我們將要另闢蹊徑，採用MCMC模擬的方式來決定一個系統終態每個晶格上自旋的排列情形。

一、Ising模型簡介

考慮一個由圖1所構成的晶格系統，晶格上的點都代表一顆原子，而黑白則表原子的自旋狀態，我們用黑表示向下而白表示向上，所以這個binary系統可以用更簡潔的\(\pm 1\)來表示每個原子可能的自旋 (即這個系統的樣本空間\(\Omega=\{+1,-1\}\))，另外為了方便，我們採用物理學的習慣將自旋表示為希臘字母\(\sigma_{n}\)，下標\(n\)可以想像成在n處原子的自旋態。原則上討論的都是一個由費米子 (fermions) 構成的系統，所以自旋應該用\(\pm \frac{1}{2}\)來表述，但在計算機上我們用\(\pm 1\)比較乾脆且簡潔。另外補充一下冷知識自然界已知的物質幾乎由費米子構成，僅有媒介作用力的粒子是玻色子 (bosons)，它們的差別在於費米子有半整數 (即不能約化成整數的分數) 自旋而玻色子則是整數自旋。

圖1. 晶格系統，黑白點分別代表該處原子的自旋狀態，黑表向下，白表向上。紅色的箭頭指出中心原子的4位鄰居\(m\)

　Ising模型假定了每個晶格上的原子的自旋都會受其鄰居 (neighbors \(\mathcal{N}\)) 的狀態影響，例如圖1中間白色的原子雖然他在這一刻是處在\(+1\)的狀態，但根據其週邊四個鄰居 (分別是3白1黑) 的況，它仍有一定的機會會在下一刻翻轉自旋的狀態。
　對一個給定位置n的原子，它和它的周圍鄰居構成的局部區域 (local area) 形成的一個系統的能量由Ising模型給出 \[ E_n = -\sum_{m\in \mathcal{N}} J \sigma_n \sigma_m - h\sigma_n. \tag{1} \] 當然對整個區域來說，我們需要也把n加總起來才會得到全部系統的總能量，但我們目前只看一個特定的n附近的狀態即可 (\(\sum_m\)表示加總所有n的鄰居m)。在Eq. (1)裡\(J=\pm 1\)，\(+1\)表示相鄰兩原子喜歡處在同一種狀態而\(-1\)比表示不喜歡處在同一態，\(h\)表示外加場 (external field) 的強度 (通常指的是磁場)。
　對n而言，在下一刻它可以處在\(\sigma_n =\pm 1\)，至於在\(1\)或在\(-1\)的機率與它此刻的狀態無關，僅與週邊的鄰居有關 (conditional on \(\mathcal{N}\))，如果還記得統計力學中系統處在某處的機率與他在該處的能量存在關係式 \[ P(\sigma_n|\mathcal{N}) = \frac{e^{-E_n/k_B T}}{\mathcal{Z}},\tag{2} \] 其中\(k_B\)是Boltzmann constant，不失一般性我們會取1，而\(T\)是系統溫度，配分函數\(\mathcal{Z}=P(+1|\mathcal{N})+P(-1|\mathcal{N})\)，或者想成歸一化係數。而上式表示當該自旋態會使局域系統的能量較低時它出現的機率會隨之上升。所以我們可以得到在n處的原子在下一刻處在\(\sigma_n=+1\)時的機率為 (MacKay 2003) \[ p(+1|\mathcal{N})=\frac{1}{1+e^{-2b_n/T}}\tag{3} \] 其中 \[ b_n = \sum_{m\in \mathcal{N}} J\sigma_m +h,\tag{4} \] 當然歸一化保證了\(p(-1|\mathcal{N})=1-p(+1|\mathcal{N})\)，所以我們只要知道處在\(\sigma_n=+1\)的機率就足夠了。
　舉例若我們計算出\(p(+1)=0.7\)，則在下一刻該處是否會翻成自旋為\(+1\)便是由一個隨機試驗決定，有0.7的機會會翻成\(+1\)而0.3的機率會變成\(-1\)。

二、MCMC計算Ising模型

請參閱本文的python腳本，請參閱以得到更多的編程細節，另外在這個腳本內我們需要引入比之前多的模組，例如在scipy.stats底下的bernoulli函數，這我們稍候再解釋。首先我們先建立一個方形的陣列系統，陣列上元素的值都是在\(\pm 1\)間取值，代表該系統上每個原子的自旋態。可以透過如下的方式來生成： 
np.random.seed(8787)
ising = np.random.choice([1,-1],p=[0.5,0.5],size=(100,100))
生成的結果如圖2所示，是一個\(100\times 100\)的陣列，但我們加上了seed random，以控制每次的初始狀況都是相同的，這方便後面比較。

圖2. 一個\(100\times 100\)的隨機自旋陣列

　然或我們必須建立一個函數用來計算在給定陣列的位置它在下一刻跳至\(+1\)的機率為何，這可以從Eq. (4)直接套用，因為是計算到\(+1\)，所以取名叫prob1函數： 
def prob1(s,J=1,T=1,H=0):
    """
    Generate the probability of a spin state being +1 at
    any location when the contributions from its neighbors
    are given.
    """
    return 1/(1+np.e**(-2*(J*s+H)/T))
這裡我拿掉腳本內大部份docstr的內容，有興趣請自己開來看，原則上我們將給定位置鄰近的狀態s餵進prob1裡，接著便返回在下一刻處在\(+1\)的機率。預設我們使得該處原子喜歡和鄰居處在相同的狀態即\(J=1\)，\(T\)和\(H\)則可以之後視計算需要再調整。
　那麼鄰居的狀態要怎麼算呢？從Eq. (5)我們可以馬上看出便是將鄰居的自旋態整個加起來即可，當然還需要乘上\(J\)和加上\(h\)。當給定位置時，計算周圍鄰居自旋的總和可以用scipy.signal裡的correlate2d來完成 (見IA6有更多的資訊)，我們考慮了週邊8位鄰居，所以每個位置周圍鄰居的自旋總和為： 
ising_fin = deepcopy(ising)
# Calculate the contributions from the neighborhoods at every locations
nei = correlate2d(ising_fin,[[1,1,1],\
                             [1,0,1],\
                             [1,1,1]],mode="same")
最後返回的nei在每個位置的值即是該處鄰居狀態的總和，和初始陣列有相同的大小。另外ising_fin和系統的初始陣列ising是一樣的，我們用deepcopy複製一個全新的初始態，因為ising_fin會隨著之後迭代計算不斷更新，漸漸變得和初始態愈來愈不一樣，但我們需要保留初始態以做後面的參考。
　我們接著來計算下一刻該系統的每個元素會變成什麼樣子，所以對初始陣列ising_fin做逐元素計算 (element-wise)： 
# Given the params
J,T,H = 1,1,0
row,col = ising_fin.shape

for r in range(row):
    for c in range(col):
        # Calculate the probability of jumping to +1 at loc (r,c)
        jump = prob1(nei[r,c],H=H,J=J,T=T)
        # Peforming random choice between two possible states +1 and -1,
        # although the probiblity of being at +1 is given
        ising_fin[r,c] = np.random.choice([1,-1],p=[jump,1-jump])
這兩個for回圈表示我們會掃過該陣列的每個元素，在位置[r,c]處時我們有鄰居狀態的總和nei[r,c]，透過它我們變可用prob1計算下一刻躍遷成\(+1\)的機率並存給變數jump。當然在下一刻的自旋到底是處在哪個態必須隨機選出來，這也就是我們使用np.random.choice的原因，選到\(+1\)的機率是jump而\(-1\)是1-jump。完成第一次迭代後得到圖3，其結果較初始狀態圖2稍微有序了一點。

圖3. 迭代一次後系統的狀態

　最終我們當然不可能只迭代一次，我們想看該系統在平衡時候，亦即即使再迭代更多的次數，系統中每顆原子的自旋也不再有很大的變化，我們將迭代30次的結果看初始解果放在一起比較，如下圖4，實現的代碼可以參見代碼腳本。

圖4. 迭代30次後系統每顆原子的自旋狀態 (右) 與初始狀態 (左) 比較

　不難發現系統狀態逐漸變得較有序，與其說他有序倒不如說是在該溫度 (預設\(T=1\))下這是這個系統最有可能的最大熵的分佈。初始狀態雖然熵比穩定系統的還要大，但由於要成為這麼亂的狀態還需要有夠高的溫度，這和系統當前的溫度不合，所以會自然降到在該溫度下熱平衡後所允許的最大熵的系統自旋分佈狀態。
　另外每次迭代中雖然原子在哪個態的機率可以算，但它下一刻會處在什麼態仍然是透過隨機來選擇的，並不是說到\(+1\)的機率大於到\(-1\)的機率我就強迫該原子翻轉成\(+1\)，他還是可以變成\(-1\)，端看機會和命運，這種給定機率然後隨機選擇的方式便稱做蒙地卡羅 (Monte Carlo, MC) 方法，而完成所有的迭代後，每一次迭代所得到的系統狀態按順序排列便構成一個隨機過程 (stochastic process)，由於當前的狀態只基於上一步鄰居的狀態，而與更早的狀態無關，亦即無記憶性 (memory-less)，這種隨機過程所構成的一條鏈便稱做馬可夫鏈 (Markov Chain)，兩個字合稱馬可夫鏈蒙地卡羅法 (Markov Chain Monte Carlo, MCMC)。基本上這是一個通用 (general purpose) 計算的概念，數學上有嚴謹的定義，但操作按問題的狀況有一定的彈性。

三、更有效率的代碼

第二節我們討論了如何計算系統的自旋態的方法，只不過在python上用for非常的緩慢與沒有效率，除非是非得計算這一步後才能計算下一步 (例如兩相鄰步驟之間的相依性) 這類的迭代我們才會用for，而前面對同一次迭代中逐元素計算陣列內每個位置[r,c]翻轉的機率prob1用for實在太不智，再來不知怎麼地np.random.choice也非常的慢，執行一次還不會發現，當像前一節那個\(100\times 100\)的矩陣完成一次整個陣列的迭代就必須執行10000次，便會發現其效率不彰。像如果要看30次的迭代就必須執行np.random.choice共30萬次，利用time模組計數需要費時6秒左右的時間，可想當陣列變得更大，或迭帶次數更多時間的增長會更可觀。
　這裡我們要稍微改寫一下，我們不再用for來做整個矩陣的逐元素計算，所以也不需要用np.random.choice來隨機選擇在某個位置自旋的狀態，我們要改引入scipy.stats模組內的bernoulli函數來計算。注意我們前面寫的prob1函數的輸入並不限於是純量，假設我們直接輸入整個nei陣列而非個別的元素nei[r,c]，那麼它回傳一個和nei大小相同的機率矩陣，每一個位置的值都代表著該處在下一刻自旋要翻轉成\(+1\)的機率值。這太好了，在python內只要可以向量化的計算都應該向量化以增進效能。
　所以既然prob1回傳的是一個機率陣列，那每個位置都是變成\(+1\)的機率\(p\)，那麼決定它是否跳到\(+1\)變可透過一個白努利試驗 (Bernoulli trial) 來決定。比如說當\(p=0.7\)，那麼bernoulli.rvs(p=0.7)回傳1表示它得自旋會變成\(+1\)，若回傳0表示自旋要變成\(-1\)。而bernoulli.rvs(p)的p接受陣列的形式，回傳陣列的大小與p相同，而回傳陣列的元素都是0或1，代表在該處所對應p的值的白努利試驗，元素為1表示該處自旋應取\(+1\)而元素0表自旋應取\(-1\)。我們可以透過自定義的spins函數來實現： 
def spins(p):
    """
    Randomly samples the spin state by giving the probability
    matrix.
    """
    spin = bernoulli.rvs(p=p)
    spin[spin==0] = -1
    return spin
所以輸入prob1的結果，回傳一個下一刻系統每顆原子自旋的狀態，最後我們把一次迭代的整個步驟濃縮進一個函數裡 (即上面代碼for r in... for c in的內容) 得到： 
def Ising(ini,J=1,T=1,H=0):
    """
    Generate the possible next spin states according to the
    current states.
    """
    # Calculate the contributions from neighbors
    nei = correlate2d(ini,[[1,1,1],\
                           [1,0,1],\
                           [1,1,1]],mode="same")    
    # Calculate the probability of being +1 according to the
    # neighborhood contributions
    p = prob1(nei,J=J,T=T,H=H)    
    # Samples the spin states according to p
    s = spins(p)    
    return s
便大功告成。
　接著我們重複前面的30次迭代，然後用新的方式計算並計時如下： 
ising_fin2 = deepcopy(ising)
start = time.time()
for i in range(30):
    ising_fin2 = Ising(ising_fin2)
end = time.time()
print("30 iterations costs %.3f seconds."%(end-start))
結果應會和圖4的右圖相仿，而耗費時間約是0.05秒左右，相比之前的方法，利用向量化和bernoulli.rvs的方式來隨機決定自旋方向的方式快了\(\mathcal{O}(100)\)倍左右！

圖5. \(1000\times 1000\)自旋陣列迭代100次的結果

　我們另外生成一個\(1000\times 1000\)的初始陣列，並計算100次迭代，結果如圖5所示，用向量化的方式計算所耗費時間約是14秒，難以想像如果使用for回圈估計要23分鐘才能計算完成！

四、模擬退火

回顧一下我們在Eq. (3)裡決定自旋狀態的機率其實有必要考慮到溫度的變量，所以溫度不同的時候，即使週邊的狀態一樣，落到\(\pm 1\)的機率也會發生變化，假設我們讓溫度一開始處在很高的狀態下 (eg. \(T> 10\))，並且隨著迭代次數的增加逐漸降低\(T\)的溫度，模擬冶金方法內逐步降溫的快慢來得到結晶晶粒不同的大小與對稱性，這就稱做模擬退火 (simulated annealing)，在最佳化理論內，利用模擬退火的方式預期可以找到全域極值 (global maximum) 而不會被truncated在某個局域極值 (local maximum) 內，缺點是computationally expensive。
　在下面的代碼裡，我們用迭代的次數來決定降溫的速度，起始溫度都是\(T=11\)。對於快速退火 (rapid cooling) 而言，每隔兩次迭代，溫度變成指數下降，即當迭代到第7次時溫度就降到1然後變不再降；另外緩慢退火 (slow cooling)，我們讓每經過10次迭代後才降1度，最終降到\(T=1\)時停止。在這個計算裡我還讓外場\(h=0.1\)，表示當最終平衡時，系統每顆原子的自旋會有很高的機率趨向外場的方向，若\(h\)取正，則自旋傾向變成\(-1\)。 
# T changes rapidly
ann_rapi = deepcopy(ising2)
T=11
for i in range(100):
    ann_rapi = Ising(ann_rapi,T=T,H=0.1)
    if i==2: T=8
    elif i==4: T=4
    elif i==6: T=2
    elif i==7: T=1
    else: pass

# T changes slowly
ann_slow = deepcopy(ising2)
T = 11
for i in range(100):
    ann_slow = Ising(ann_slow,T=T,H=0.1)
    if i%10 == 0: T-=1
    else: pass
結果如圖6。

圖6. 快速 (左) 與緩慢 (右) 退火 (從\(T=11\)降至\(T=1\)) 經過100次迭代後的結果。

　他們的初始態都是圖5(左)，不難發現快速退火的結果保有較多初始態那種凌亂的狀態，就像被快速凍結一樣，但慢速退火因為在降溫前整個系統都有機會與周圍的溫度達成熱平衡，或者說成該系統有足夠的時間落到該溫度下的最小的能量狀態 (若是一個逐溫度變化的凸函數，那麼在變化前系統能量都可以落到最小值的點去)，這變是為何緩慢退火每顆原子的自旋態都較為一致。這份模擬的結果也符合冶金學的結果，若要結晶漂亮且大顆那麼降溫的速度則要慢！

圖7.  棋盤格狀初始態經過退火後的情況

圖8. 加菲貓初始態經過退火後的情況

　計算機模擬最令人著迷的地方便是不受限於自然界的狀態，除了初始狀態是任意在自旋向上向下挑外，我還可以構造稀奇古怪的初始狀態，圖7和圖8我則是初始狀態排列為一棋盤格狀和加菲貓的造型，我們用同樣的快速和慢速退火得到的結果，另外在腳本內還有用ipywidgets模組的interact控制bar，可以讓使用者拖曳查看不同時刻系統的狀態，有興趣可以自己玩玩看！
　最後補充一個冷知識，這邊講的是退火，但如果我們讓溫度從低溫升回高溫的狀態，這叫作回火 (tempering)，在冶金學裡由於不同的回火速度會改變原子的排列，所以最終材料的係數 (例如韌性、剛性) 也會不一樣，不過這個概念似乎沒有應用到computer science裡。

五、小結

在這次的文章裡簡單的介紹了Ising模型，並用simulation的方式來看系統自旋態的演化過程，也比較了不同退火速度所造成的影響。Ising模型雖然是近百年前便被提出，但至今二維以上的系統狀態，如我們在這裡所模擬的結果仍然沒有解析的方程，另外模擬退火的概念也在最佳化理論、計算機視覺、機器學習等領域裡也扮演了舉足輕重的角色。



Source on github: PH3_Ising_model.ipynb

留言

張貼留言

這個網誌中的熱門文章

PS7 - 小專題: 為何高斯分佈擁有最大熵?

圖片
高爾頓板，或稱做quincunx，是由Galton爵士 (1822-1911) 所發明用來驗證中央極限定理 (CLT) 的裝置。它也是第一位討論均值回歸 (regression to the mean) 現象的科學家，但同時亦是優生學 (eugenics) 這樣錯誤觀念的推廣者 在 上篇雜記 的結尾我們提到了對於一個趨於穩態的系統擁有最大的熵，而這也是統計力學的直接應用，用物理的角度來講便是當系統平衡時 Gibbs 自由能不再變化 (\(\Delta G =0\))，而在此時會對應到系統的熵有全域極大值 (global maximum)。故大自然很多現象都已經處於穩態的狀況下，我們對這個現象的母體做抽樣，所得到的觀測樣本應也有最大熵，則描述這些觀測樣本最好的分佈便是符合該條件下擁有最大熵的機率分佈。而指數族的分佈符合最大熵的特點，我們在這個小專題內便嘗試證明高斯與二項式分佈在連續和離散的隨機變量下的條件擁有最大熵的性質。 一、機率熵 在 上篇雜記 裡裡我們探討了夏儂熵，或者稱作資訊量\(I\)的意含，所以一個系統內 (我們這裡就不再用系綜這種彆扭的歷史名詞了) 每個成員出現的機率是\(p_i\)，則對該系統的熵可以用夏儂熵公式來表示 \[ I(p) = -\sum_i p_i \ln p_i. \tag{7.1} \] 所以對某一組 連續 機率分佈函數\(p(x)\)，其中\(x\)是隨機變量，則將(7.1)的求和符號改成積分後有 \[ I(p) = -\int p(x)\ln p(x) dx, \tag{7.2} \] 稱為機率分佈\(p(x)\)的 機率熵 ，而積分的上下限則看隨機變量的範圍，例如高斯分佈就是\(\pm \infty\)。如果今天我們有兩組機率分佈\(p(x)\)和\(q(x)\)，我們便可以定義它們之間的 相對熵 ，或稱作 KL散度 (Kullback-Leibler divergence)，定義為 \[ D_{KL}(q,p)\equiv\int q(x)\ln \left(\frac{q(x)}{p(x)}\right) dx = I(q,p)-I(q)\geq 0, \tag{7.3} \] 是 正定的 (positively defined)，其中 \[ I(q,p)=-\int q(x)\ln p(x) d
閱讀完整內容

IA9a - K-means分群與EM算法: 理論

圖片
市場內分門別類的產品 原則上 K-means分群 是 EM算法 (expectation-maximization algorithm) 的一個特例，相較於EM從似然函數 (likelihood) 的方式出發並導出不同資料點對不同群 (或稱聚類、簇，cluster) 的貢獻度，K-means單純用歐氏距離來分類，雖然這加快了K-means算法的計算效率，但顯然的這麼naïve的方法會產生明顯的缺陷。但仍讓我們會先從簡單的K-means開始，並探討K-means算法的結果在某些情況下和人類肉眼辨識的結論會產生明顯的歧義，接著再推進到通用的EM算法，除了 python 的demo外，我也會證明在數學上如何從EM得到K-means的極限 (僅限高斯混合的情況)，而在這個part內，讓我們先從數學理論的角度來面。 一、K-means分群 今天如果我有一大坨資料共120個點如圖1所示，從資料的散佈趨勢來看，很明顯的分成兩群，那我們是不是可以有一個演算法來自動化將資料點如同我們臆測的一樣分成兩袋？一個最簡單的發法便是使用K-means分群，而這個大寫K表示群的數量，以圖1而言，我們猜測K=2，所以是個2-means分群問題。 圖1. 資料散佈圖 (普魯士藍) 與初始化的紅綠兩點 　所以為了達成分群，我們需要先初始化兩個點，這兩個點代表兩群的中心。事實上，初始化的位置不必要一定剛好在兩群的中心，可以是任意挑選的地方，但為了計算方便，我們會將起始化的兩點開始在靠近資料點的平均值附近，如同圖1中紅綠兩點。我們將這兩個初始化的點叫作\(p_k\)其中\(k=1,2\)表示第\(k\)群。另外將資料點的集合用\(\mathbf{D}=\{d_1,d_2\cdots d_n \cdots d_{120}\}\)表示，所以我們計算每個點\(d_n=(x_{d_n},y_{d_n})\)和\(p_k=(x_{p_k},y_{p_k})\)的歐氏距離 \[ \ell_{n,k}=\sqrt{(x_{d_n}-x_{p_k})^2+(y_{d_n}-y_{p_k})^2}. \tag{1} \] 若該點和\(p_1\)的距離小於和\(p_2\)的距離 (\(\ell_{n,1}<\ell_{n,2}\))，則我們說\(d_n\)比較靠近\(p_1\)且應該
閱讀完整內容