李航,第十一章,條件隨機場html
攜代碼:用 Python 經過馬爾可夫隨機場(MRF)與 Ising Model 進行二值圖降噪【推薦!】網絡
CRF:http://www.jianshu.com/p/55755fc649b1app
模型
------dom
首先什麼是隨機場呢,一組隨機變量,他們樣本空間同樣,那麼就是隨機場。當這些隨機變量之間有依賴關係的時候,對咱們來講纔是有意義的。分佈式
咱們利用這些隨機變量之間的關係建模實際問題中的相關關係,實際問題中咱們可能只知道這兩個變量之間有相關關係,但並不知道具體是多少,咱們想知道這些依賴關係具體是什麼樣的,因而就把相關關係圖畫出來,而後經過實際數據訓練去把具體的相關關係訓練出來嵌入到圖裏,而後用獲得的這個圖去進行預測、去進行reference等不少事情。函數
那麼爲了簡化某些爲問題來講,也爲了這個圖畫出來能用,咱們會在畫圖的時候要遵循一些假設和規則,好比馬爾科夫獨立性假設。按照這個假設和規則來畫圖,畫出來的圖會知足一系列方便的性質便於使用。post
馬爾可夫獨立性假設是說:對一個節點,在給定他所鏈接的全部節點的前提下,他與外接是獨立的。就是說若是你觀測到了這個節點直接鏈接的那些節點的值的話,那他跟那些不直接鏈接他的點就是獨立的。形式上,咱們是想把他設計成這個樣子的,邊能夠傳遞信息,點與點之間經過邊相互影響,若是觀測到一個節點的取值或者這個節點的取值是常量,那麼別的節點就沒法經過這個節點來影響其餘節點。因此對一個節點來講,若是用來鏈接外界的全部節點都被鎖住了,那他跟外界就沒法傳遞信息,就獨立了。這比貝葉斯網絡就直觀多了,貝葉斯網絡要判斷兩點之間獨立還要看有沒有v-structure,還要看邊的指向。學習
吶,知足馬爾可夫獨立性的隨機場,就叫馬爾可夫隨機場。它不只具備我剛纔說的那些性質,除此以外,還等價於吉布斯分佈。優化
這些邊具體是如何建模的呢,以什麼形式記錄這些機率信息的?貝葉斯網絡每一條邊是一個條件機率分佈,P(X|Y),條件是父節點、結果是子節點。他有一個問題,就是當我知道A、B、C三個變量之間有相關關係,可是不知道具體是誰依賴誰,或者我不想先假設誰依賴誰,這個時候貝葉斯就畫不出來圖了。由於貝葉斯網絡是經過變量之間的條件分佈來建模整個網絡的,相關關係是經過依賴關係(條件分佈)來表達的。而馬爾可夫隨機場是這樣,我不想知道這三個變量間究竟是誰依賴誰、誰是條件誰是結果,我只想用聯合分佈直接表達這三個變量之間的關係。好比說兩個變量A、B,這兩個變量的聯合分佈是:
| A, B | P(A, B) |
|--------------+---------|
| A = 0, B = 0 | 100 |
| A = 0, B = 1 | 10 |
| A = 1, B = 0 | 20 |
| A = 1, B = 1 | 200 |
【這裏與有向圖不同的地方是,P(a,b)再也不是用機率表示,但這並不會帶來什麼影響】
這個分佈表示,這條邊的功能是使它鏈接的兩點(A和B)趨同,當A = 0的時候B更可能等於0不太可能等於1,當A = 1的時候B更可能等於1不太可能等於0。這樣一來你知道了三個變量之間的聯合分佈,那他們兩兩之間的條件分佈天然而然就在裏面。
這樣出來的圖是等價于吉布斯分佈的,就是說,你能夠只在每一個最大子團上定義一個聯合分佈(而不須要對每一個邊定義一個聯合分佈),整個圖的聯合機率分佈就是這些最大子團的聯合機率分佈的乘積。固然這裏最大子團的聯合機率並非標準的聯合機率形式,是沒歸一化的聯合機率,叫factor(因子),整個圖的聯合機率乘完以後下面再除一個歸一化因子和就歸一化了,最終是一個聯合機率,每一個子團記載的都是因子,是沒歸一化的機率,嚴格大於零,能夠大於一。但關鍵是依賴關係、這些相關關係已經encode在裏面了。
條件隨機場是指這個圖裏面一些點我已經觀測到了,求,在我觀測到這些點的前提下,整張圖的分佈是怎樣的。就是given觀測點,你去map inference也好你去作之類的事情,你可能不求具體的分佈式什麼。這裏還要注意的是,馬爾科夫隨機場跟貝葉斯網絡同樣都是產生式模型,條件隨機場纔是判別式模型。
這是條件隨機場,NER(命名實體識別)這個任務用到的是線性鏈條件隨機場。
線性鏈條件隨機場的形式是這樣的,觀測點是你要標註的這些詞自己和他們對應的特徵,例如說詞性是否是專有名詞、語義角色是否是主語之類的。隱節點,是這些詞的標籤,好比說是否是人名結尾,是否是地名的開頭這樣。這些隱節點(就是這些標籤),依次排開,相鄰的節點中間有條邊,跨節點沒有邊(線性鏈、二階)。而後全部觀測節點(特徵)同時做用於全部這些隱節點(標籤)。至於觀測節點之間有沒有依賴關係,這些已經不重要了,由於他們已經被觀測到了,是固定的。
這是線性鏈條件隨機場的形式。
吶,這些特徵是怎麼表達的呢?是這樣,他有兩種特徵:
特徵表達形式比較簡單,就是你是否知足我特徵所說的這個配置,是就是1,不是就是0。好比說,上一個狀態是地名的中間,且當前詞是'國'(假設他把中國分詞 拆成兩個了),且當前詞的詞性是專有名詞、且上一個詞的詞性也是專有名詞,若是知足這個配置、輸出就是一、不知足就輸出0。而後這些特徵每一個都有一個權重,咱們最後要學的就是這些權重。特徵跟權重乘起來再求和,外面在套個exp,出來就是這個factor的形式。這是一個典型的對數線性模型的表達方式。這種表達方式很是常見,有不少好處,好比爲何要套一個exp呢?
推斷
------
線性鏈的條件隨機場跟線性鏈的隱馬爾科夫模型同樣,通常推斷用的都是維特比算法。這個算法是一個最簡單的動態規劃。
首先咱們推斷的目標是給定一個X,找到使P(Y|X)最大的那個Y嘛。而後這個Z(X),一個X就對應一個Z,因此X固定的話這個項是常量,優化跟他不要緊(Y的取值不影響Z)。而後 exp也是單調遞增的,也不帶他,直接優化exp裏面。因此最後優化目標就變成了裏面那個線性和的形式,就是對每一個位置的每一個特徵加權求和。好比說兩個狀態的話,它對應的機率就是從開始轉移到第一個狀態的機率加上從第一個轉移到第二個狀態的機率,這裏機率是隻exp裏面的加權和。那麼這種關係下就能夠用維特比了,首先你算出第一個狀態取每一個標籤的機率,而後你再計算到第二個狀態取每一個標籤得機率的最大值,這個最大值是指從狀態一哪一個標籤轉移到這個標籤的機率最大,值是多 少,而且記住這個轉移(也就是上一個標籤是啥)。而後你再計算第三個取哪一個標籤機率最大,取最大的話上一個標籤應該是哪一個。以此類推。整條鏈計算完以後, 你就知道最後一個詞去哪一個標籤最可能,以及去這個標籤的話上一個狀態的標籤是什麼、取上一個標籤的話上上個狀態的標籤是什麼,醬。這裏我說的機率都是 exp裏面的加權和,由於兩個機率相乘其實就對應着兩個加權和相加,其餘部分都沒有變。
學習
------
這是一個典型的無條件優化問題,基本上全部我知道的優化方法都是優化似然函數。典型的就是梯度降低及其升級版(牛頓、擬牛頓、BFGS、L-BFGS),這裏版本最高的就是L-BFGS了吧,因此通常都用L-BFGS。除此以外EM算法也能夠優化這個問題。
【有點難理解的東西,仍是直接實戰的好】
數學表示:
既然噪聲圖是從原圖添加噪聲而來,咱們擁有了先驗知識 1:
yi和xi有很強的聯繫。
通常圖片裏,每一個像素和與它相鄰的像素值應當較爲接近,好比上圖中的黑色筆畫和白色負空間,除了邊緣之外,黑色的像素周圍都是黑色像素,白色像素的周圍都是白色像素(連成一片))。這樣咱們就獲得了先驗知識 2:
xi和與它相鄰的其餘像素也存在較強的聯繫
若是咱們狠一點,假設原圖像素只與它的直接相鄰像素有聯繫(即具有條件獨立性質),咱們就能夠獲得一個具有局部馬爾可夫性質(Local Markov property)的圖模型
如此,找」最大團「也相對簡單,都是成對的結點。在這樣一個圖模型裏,咱們有兩種團(clique):
這兩種團合併起來,獲得的 {xi, yi, xj} 顯然是一個最大團(Maximal Clique),此時咱們能夠利用它來對這個馬爾可夫隨機場進行 factorization,即求得其聯合機率分佈關於最大團 xC = {xi, yi, xj} 的函數。
【這裏又最大團又變成了三個結點構成】
其中 Z 爲 partition function,是 p(x) 的歸一化常數(normalization constant),求法參見 PRML 8.3.2。由於與咱們的實現不相關,這裏不贅述。
ψC(xC) 即所謂的 potential function,爲了方便咱們一般只求它的對數形式 E(xC)(按照其物理意義稱爲 energy function)
關於 factorization 的過程和推導能夠參見 PRML 8.3.2,這裏咱們須要作的是定義一個 energy function:
Step 1. 由於咱們須要處理的是二值圖,首先咱們定義 xi ∈ {-1, +1},假設這裏白色爲1,黑色爲-1。
對於原圖像素與噪聲圖像素構成的團 {xi, yi},咱們定義一項 −ηxiyi,其中 η 爲一個非負的權重。當二者同號時【1*1 or (-1)*(-1),表示噪聲圖像素與原圖像素較爲接近時】,這項的值爲−η,異號時爲 η。這一項可以下降 energy(由於爲負)。
Step 2. 對於噪聲圖中相鄰像素構成的團 {xi, xj},咱們定義一項 −βxixj,其中 β 爲一個非負的權重。這樣,當處理過的圖像裏相鄰的像素較爲接近時,這一項可以下降 energy(由於爲負)。
Step 3. 最後,咱們再定義一項 hxi,使處理過的圖像總體偏向某一個像素值。
對圖像中的每個像素,將這三項加起來,就獲得咱們的 energy function:
對應聯合機率:
顯然 energy 越低,降噪過的圖像與原圖一致的機率越高。
(注意由於咱們這裏求的 E 已經對整個矩陣求和,即對應 potential function 的積,因此計算聯合機率分佈的時候不須要再求積)
代碼分析:
使用 Python 實現這個 energy function 時,咱們可使用一個 closure 來實現一個 function factory,經過傳遞beta
(β),eta
(η)和 h
參數,生成對應的 energy function。
此外爲了方便,咱們假設傳入的x
和y
不是一維向量,而是對應圖像的二維矩陣(注意是np.ndarray
而不是nd.matrix
,前者的*
纔是array multiplication即逐個元素相乘,後者的*
是矩陣乘法)。
import numpy as np def E_generator(beta, eta, h): """Generate energy function E. Usage: E = E_generator(beta, eta, h) Formula: E = h * \sum{x_i} - beta * \sum{x_i x_j} - eta * \sum{x_i y_i} """
def E(x, y): """Calculate energy for matrices x, y. Note: the computation is not localized, so this is quite expensive. """
# sum of products of neighboring paris {xi, xj},參考「技巧解釋」
# 由於邊界元素不必定有四個鄰居,−βxixj這項存在邊界問題,咱們須要特別處理,
# 利用 numpy 的 fancy index,寫起來並不困難,以下。
xxm = np.zeros_like(x) xxm[ :-1,: ] = x[1: , : ] # down
xxm[1:, : ] += x[ :-1, : ] # up
xxm[ :, :-1] += x[ : ,1: ] # right
xxm[ :, 1: ] += x[ : , :-1] # left
xx = np.sum(xxm * x) xy = np.sum(x * y) xsum = np.sum(x)
return h * xsum - beta * xx - eta * xy return E
技巧:注意到若是用 xi0 ~ xi3 表示 xi 的四個鄰居,則 xi * xi0 + xi * xi1 + xi * xi2 + xi * xi3 = + xi * (xi1 + ... + xi3),即乘法結合律,所以咱們能夠先將鄰居相加,再與 x
相乘。
注意這裏生成的E
每次都要對矩陣中的全部元素進行運算,因此即便有 numpy 加持,開銷依然較大。後面咱們會按照需求進行優化。
基本思想:
注意若是咱們固定 y 做爲先驗知識(假設噪聲圖不變),咱們所求的機率就變成了 p(x|y),這種模型叫作 Ising Model,【這不就成了CRF了麼?】
在統計物理中有普遍的應用。這樣咱們的問題就成了以 y 爲基準,想辦法不斷變更 x,而後找出最接近原圖的 x。
一種最簡單的辦法是:(ICM)
這種方法稱爲 Iterated Conditional Modes(ICM),由 Julian Besag 在 1986 年的論文 On the Statistical Analysis of Dirty Pictures 中提出(這篇論文在 80 年代英國數學家所著論文裏引用數排名第一……)。
由於 ICM 的每一步實際上固定住了其餘元素,只變更當前遍歷到的那個元素,因此咱們能夠將 E的計算 localize,只對受影響的那一小片區域從新計算。
咱們可讓 function factory E_generator
返回兩個版本的 E:
這裏添加了兩個函數:
def E_generator(beta, eta, h): def E(x, y): ... # as shown before
def is_valid(i, j, shape): """Check if coordinate i, j is valid in shape."""
return i >= 0 and j >= 0 and i < shape[0] and j < shape[1] def localized_E(E1, i, j, x, y): """Localized version of Energy function E. Usage: old_x_ij, new_x_ij, E1, E2 = localized_E(Ecur, i, j, x, y) """ oldval = x[i, j] newval = oldval * -1 # flip
# local computations
E2 = E1 - (h * oldval) + (h * newval) E2 = E2 + (eta * y[i, j] * oldval) - (eta * y[i, j] * newval) adjacent = [(0, 1), (0, -1), (1, 0), (-1, 0)] neighbors = [x[i + di, j + dj] for di, dj in adjacent if is_valid(i + di, j + dj, x.shape)] E2 = E2 + beta * sum(a * oldval for a in neighbors) E2 = E2 - beta * sum(a * newval for a in neighbors) return oldval, newval, E1, E2 return E, localized_E
ICM實現:
def ICM(y, E, localized_E): """Greedy version of simulated_annealing().""" x = np.array(y) Ebest = Ecur = E(x, y) # initial energy
initial_time = time.time() energy_record = [[0.0, ], [Ebest, ]] accept, reject = 0, 0 for idx in np.ndindex(y.shape): # for each pixel in the matrix
old, new, E1, E2 = localized_E(Ecur, idx[0], idx[1], x, y) if (E2 < Ebest): Ecur, x[idx] = E2, new Ebest = E2 # update Ebest
else: Ecur, x[idx] = E1, old # record time and Ebest of this iteration
end_time = time.time() energy_record[0].append(end_time - initial_time) energy_record[1].append(Ebest) return x, energy_record
能夠看到大約 96% 的像素與原圖一致。不過光從視覺上看,降噪事後的圖依然有很多明顯的噪點,這是由於 ICM 採起的相似貪心的策略使得它容易陷入局部最優(local optimum)。
若是想要獲得一個更好的降噪結果,咱們顯然要採起一種可以獲得全局最優(global optimum)的策略。
故,可使用模擬退火的辦法。
其實這個問題能夠當作一個搜索問題:
「在全部 x 的兩種狀態組成的狀態空間裏(假設有 n 個像素,那麼狀態空間大小爲 2n),找到能使 energy 最低的狀態。「
因爲狀態空間大小呈指數級增加,僅僅是對於一個有 240 × 180 = 43,200 個像素的圖片來講,這個狀態空間就已經不可能使用暴力搜索解決了(其實是 NP-Hard 的),因此咱們須要考慮其餘的搜索策略。
這裏咱們能夠嘗試使用模擬退火(Simulated annealing),在有限的時間內找到儘量好的解。本方法由 Stuart Geman 與 Donald Geman (兄弟喲) 在 1984 年的論文 Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images 中提出,而且文中對模擬退火可收斂至全局最優進行了詳細的證實(這篇論文和以前 Julian Besag 的那篇都是 ISI Highly Cited Paper,引用數至關魔性……)。
模擬退火須要一個控制溫度的 schedule,具體怎樣能夠本身調,只須要在迭代次數 k 接近最大迭代次數 kmax 時逼近 0 便可,這裏設定爲:
def temperature(k, kmax): """Schedule the temperature for simulated annealing."""
return 1.0 / 500 * (1.0 / k - 1.0 / kmax)
模擬退火的核心在於:當局部變化致使情況惡化(這裏爲能量變大)時,依據當前溫度 t,設該變化對「全局最優」有利的機率爲 e△E/t,按照這個機率來肯定是否保留這個變化,即所謂的 transition probabilities,以下
def prob(E1, E2, t): """Probability transition function for simulated annealing."""
return 1 if E1 > E2 else np.exp((E1 - E2) / t)
模擬退火實現:
def simulated_annealing(y, kmax, E, localized_E, temp_dir): """Simulated annealing process for image denoising. Parameters ---------- y: array_like The noisy binary image matrix ranging in {-1, 1}. kmax: int The maximun number of iterations. E: function Energy function. localized_E: function Localized version of E. temp_dir: path Directory to save temporary results. Returns ---------- x: array_like The denoised binary image matrix ranging in {-1, 1}. energy_record: [time, Ebest] records for plotting. """ x = np.array(y) Ebest = Ecur = E(x, y) # initial energy
initial_time = time.time() energy_record = [[0.0, ], [Ebest, ]] for k in range(1, kmax + 1): # iterate kmax times
start_time = time.time() t = temperature(k, kmax + 1) print "k = %d, Temperature = %.4e" % (k, t) accept, reject = 0, 0 # record the acceptance of alteration
for idx in np.ndindex(y.shape): # for each pixel in the matrix
old, new, E1, E2 = localized_E(Ecur, idx[0], idx[1], x, y) p, q = prob(E1, E2, t), random() if p > q: accept += 1 Ecur, x[idx] = E2, new if (E2 < Ebest): Ebest = E2 # update Ebest
else: reject += 1 Ecur, x[idx] = E1, old # record time and Ebest of this iteration
end_time = time.time() energy_record[0].append(end_time - initial_time) energy_record[1].append(Ebest) print "k = %d, accept = %d, reject = %d" % (k, accept, reject) print "k = %d, %.1f seconds" % (k, end_time - start_time) # save temporary results
temp = sign(x, {-1: 0, 1: 255}) temp_path = os.path.join(temp_dir, 'temp-%d.png' % (k)) Image.fromarray(temp).convert('1', dither=Image.NONE).save(temp_path) print "[Saved]", temp_path return x, energy_record
Result:
不管是 ICM 仍是模擬退火,都是經過最小化能量來找到原圖的近似。
後來 Greig, Porteous 和 Seheult 提出了使用 graph cuts 的模型,將能量值的最小化轉化爲最大化解的後驗估計(a posteriori estimate),進而轉爲計算機科學裏常見的 max-flow/min-cut 的問題,求解後可以獲得更好的效果
(參考 D.M. Greig, B.T. Porteous and A.H. Seheult (1989), Exact maximum a posteriori estimation for binary images, Journal of the Royal Statistical Society Series B, 51, 271–279.)。
From: http://blog.echen.me/2012/01/03/introduction-to-conditional-random-fields/
舉個例子,假若有一張小明閉着嘴的照片,怎麼分類?顯然難以直接判斷,須要參考閉嘴以前的照片,若是以前的照片顯示小明在吃飯,那這個閉嘴的照片極可能是小明在咀嚼食物準備下咽,能夠給它打上吃飯的標籤;若是以前的照片顯示小明在唱歌,那這個閉嘴的照片極可能是小明唱歌瞬間的抓拍,能夠給它打上唱歌的標籤。
因此,爲了讓咱們的分類器可以有更好的表現,在爲一張照片分類時,咱們必須將與它相鄰的照片的標籤信息考慮進來。這——就是條件隨機場(CRF)大顯身手的地方!
給一個句子中的每一個單詞註明詞性。好比這句話:「Bob drank coffee at Starbucks」,註明每一個單詞的詞性後是這樣的:「Bob (名詞) drank(動詞) coffee(名詞) at(介詞) Starbucks(名詞)」。
下面,就用條件隨機場來解決這個問題。
以上面的話爲例,有5個單詞,咱們將:(名詞,動詞,名詞,介詞,名詞)做爲一個標註序列,稱爲l,可選的標註序列有不少種,咱們要在這麼多的可選標註序列中,挑選出一個最靠譜的做爲咱們對這句話的標註。
假如咱們給每個標註序列打分,打分越高表明這個標註序列越靠譜,咱們至少能夠說,凡是標註中出現了動詞後面仍是動詞的標註序列,要給它負分!!
上面所說的動詞後面仍是動詞就是一個特徵函數,咱們能夠定義一個特徵函數集合,用這個特徵函數集合來爲一個標註序列打分,並據此選出最靠譜的標註序列。
也就是說,每個特徵函數均可以用來爲一個標註序列評分,把集合中全部特徵函數對同一個標註序列的評分綜合起來,就是這個標註序列最終的評分值。
特徵函數,就是這樣的函數,它接受四個參數:【輸入】
輸出值是0或者1:【輸出】
定義好一組特徵函數後,咱們要給每一個特徵函數fj賦予一個權重λj。如今,只要有一個句子s,有一個標註序列l,咱們就能夠利用前面定義的特徵函數集來對l評分。
m: 特徵函數集中的個數
n: 句子單詞個數
對這個分數進行指數化和標準化,咱們就能夠獲得標註序列l的機率值p(l|s),以下所示:
當l_i是「副詞」而且第i個單詞以「ly」結尾時,咱們就讓f1 = 1,其餘狀況f1爲0。不難想到,f1特徵函數的權重λ1應當是正的。並且λ1越大,表示咱們越傾向於採用那些把以「ly」結尾的單詞標註爲「副詞」的標註序列。
若是i=1,l_i=動詞,而且句子s是以「?」結尾時,f2=1,其餘狀況f2=0。一樣,λ2應當是正的,而且λ2越大,表示咱們越傾向於採用那些把問句的第一個單詞標註爲「動詞」的標註序列。
當l_i-1是介詞,l_i是名詞時,f3 = 1,其餘狀況f3=0。λ3也應當是正的,而且λ3越大,說明咱們越認爲介詞後面應當跟一個名詞。
一個條件隨機場就這樣創建起來了,讓咱們總結一下:
爲了建一個條件隨機場,咱們首先要定義一個特徵函數集,每一個特徵函數都以整個句子s,當前位置i,位置i和i-1的標籤爲輸入。
而後爲每個特徵函數賦予一個權重,而後針對每個標註序列l,對全部的特徵函數加權求和,必要的話,能夠把求和的值轉化爲一個機率值。
HMM和CRF的關係:
每個HMM模型都等價於某個CRF
可是,CRF要比HMM更增強大,緣由主要有兩點:
CRF能夠定義數量更多,種類更豐富的特徵函數。HMM模型具備自然具備局部性,就是說,在HMM模型中,當前的單詞只依賴於當前的標籤,當前的標籤只依賴於前一個標籤。這樣的局部性限制了HMM只能定義相應類型的特徵函數,咱們在上面也看到了。可是CRF卻能夠着眼於整個句子s定義更具備全局性的特徵函數,如這個特徵函數:
![]()
若是i=1,l_i=動詞,而且句子s是以「?」結尾時,f2=1,其餘狀況f2=0。
CRF可使用任意的權重 將對數HMM模型看作CRF時,特徵函數的權重因爲是log形式的機率,因此都是小於等於0的,並且機率還要知足相應的限制,如
但在CRF中,每一個特徵函數的權重能夠是任意值,沒有這些限制。