馬爾科夫鏈, PageRank演算法, 與代數圖論

馬爾科夫鏈是一種時間離散的隨機過程。因為時間是離散的，所以我們可以用整數來表示時間指標。在每一個時刻，馬爾科夫鏈都可以用一個狀態向量 $oldsymbol{xi}_{n}$ 來表示。隨著時間的推移，狀態向量也在演化。我們可以用一個矩陣將時刻到時刻狀態向量聯繫起來:

$oldsymbol{xi}_{n + 1} = M oldsymbol{xi}_{n}$

這個矩陣叫做馬爾科夫矩陣。因為它表示不同時刻狀態間的轉移，所以又叫做轉移矩陣。該矩陣應該滿足一系列的性質，這些性質將會在下面舉例說明。

馬爾科夫鏈的一個特徵是無記憶性，也就是馬爾科夫鏈的演化只依賴於當前狀態，而不依賴於過去的歷史。現實世界中很多過程都可以簡化為一個馬爾科夫鏈。一個例子就是一維的隨機遊走。假設有一個被限制在一維空間的粒子，這個粒子每次以概率向左移動一步，以概率向右移動一步，並且每次粒子移動都不依賴於過去的歷史，那麼這個粒子的運動就是一個馬爾科夫過程。進一步假設粒子的運動步長是固定的，也就是粒子的位置取值只能是離散的，那麼這又是一個馬爾科夫鏈。

第一節：一個簡單的馬爾科夫鏈

我們很有可能在高中時就已經接觸到馬爾科夫鏈了。這裡有一道高中數學問題：

某城市有兩個超市，分別記作A和B. 每次去A超市的人有20%下次仍然去A超市，其餘的80%去B超市；每次去B超市的人有40%下次仍然去B超市，其餘的60%去A超市。請問，經過足夠長的時間後，去這A超市的人口和去B超市的人數比例會不會趨近於一個常數？如果會，請問這個常數為多少？

我們可以記第次去A超市的人數為，去B超市的人數為 . 於是，我們可以得到一個遞歸關係：

$egin{align} & a_{n+1} = 0.2 a_{n} + 0.6 b_{n} \ & b_{n+1} = 0.8 a_{n} + 0.4 b_{n} end{align}$

用矩陣表示為

$egin{pmatrix} a_{n+1} \ b_{n+1} end{pmatrix} = egin{pmatrix} 0.2 & 0.6 \ 0.8 & 0.4 end{pmatrix} egin{pmatrix} a_{n} \ b_{n} end{pmatrix} , n ge 1$

這裡得到了一個轉移矩陣 $M = egin{pmatrix} 0.2 & 0.6 \ 0.8 & 0.4 end{pmatrix}$ .

使用遞歸，很容易得到

$egin{pmatrix} a_{n+1} \ b_{n+1} end{pmatrix} = egin{pmatrix} 0.2 & 0.6 \ 0.8 & 0.4 end{pmatrix}^{n} egin{pmatrix} a_{1} \ b_{1} end{pmatrix}$

這是一個簡單的馬爾科夫鏈。對於每一個，馬爾科夫鏈的狀態都可以用兩個數字 $a_{n}, b_{n}$ 完全表示出來。時刻的狀態可以用時刻的狀態唯一地表示出來，而不依賴於此前的歷史，這正是馬爾科夫鏈的性質。很容易證明，去A超市的人數與去B超市的人數的和是一個常數，因為：

$egin{pmatrix} 1 & 1 end{pmatrix} egin{pmatrix} a_{n+1} \ b_{n+1} end{pmatrix} = egin{pmatrix} 1 & 1 end{pmatrix} egin{pmatrix} 0.2 & 0.6 \ 0.8 & 0.4 end{pmatrix} egin{pmatrix} a_{n} \ b_{n} end{pmatrix} = egin{pmatrix} 1 & 1 end{pmatrix} egin{pmatrix} a_{n} \ b_{n} end{pmatrix}$

這裡已經涉及到了馬爾科夫矩陣的一個重要性質，就是該矩陣的每一列的元素求和都是1. 因此，馬爾科夫矩陣滿足這個條件：

$oldsymbol{1}^{T} M = oldsymbol{1}^{T}$

其中, 是元素全都是1的長度等於馬爾科夫矩陣列數的列向量。很容易證明， $M^{2}$ 也滿足這個條件。進一步可以證明，的任意次正數冪都滿足這個條件。將這條件兩邊取轉置，得到

$M^{T} oldsymbol{1} = oldsymbol{1}$

於是我們得到了矩陣的一個特徵值為1，它所對應的特徵向量為 . 因為和有相同的特徵值（特徵向量未必一樣），所以我們知道1也是馬爾科夫矩陣的特徵值。根據 Gershgorin circle theorem，我們知道馬爾科夫矩陣的特徵值不可能大於1. 因為該矩陣的每一個元素都大於0, 根據Perron-Frobenius theorem，我們知道該矩陣除了一個等於1的特徵值之外，其餘所有的特徵值的絕對值都小於1. 對於矩陣 $M = egin{pmatrix} 0.2 & 0.6 \ 0.8 & 0.4 end{pmatrix}$ ，我們可以將它對角化得到

$M = U egin{pmatrix} 1 & 0 \ 0 & -0.4end{pmatrix} U^{-1}$

其中 $U = egin{pmatrix}3 & 1 \ 4 & -1 end{pmatrix}$ .

馬爾科夫矩陣的次冪為

$M^n = U egin{pmatrix} 1 & 0 \ 0 & (-0.4)^nend{pmatrix} U^{-1}$

當冪為無窮大的時候，我們得到

$lim_{n ightarrowinfty} M^n = U egin{pmatrix}1 & 0 \ 0 & 0 end{pmatrix} U^{-1} = egin{pmatrix} 3/7 & 3/7 \ 4/7 & 4/7end{pmatrix}$

所以經過足夠長的時間後，去A超市的人數趨近於 $a_{infty} = frac{3}{7}(a_1 + b_1)$ ，去B超市的人數趨近於 $b_{infty} = frac{4}{7} (a_1 + b_1)$ . 我們稱狀態向量 $egin{pmatrix} a_{infty} \ b_{infty}end{pmatrix}$ 為馬爾科夫鏈的穩態，經過歸一化，我們知道穩態向量不依賴於初始條件。最後去A超市的人數與B超市的人數的比例為3:4.

第二節：馬爾科夫矩陣

上一節已經引入了馬爾科夫矩陣，並用這個矩陣求解了一個簡單的馬爾科夫鏈問題。這一節裏，我要將上一節的內容系統化。因為馬爾科夫矩陣表示了狀態間的轉移概率，所以馬爾科夫矩陣的每一個元素都應該是非負的。記馬爾科夫矩陣元素為 $M_{ij}$ . 該元素表示系統從狀態轉移到狀態的概率。例如，上一節中馬爾科夫矩陣為 $M = egin{pmatrix} 0.2 & 0.6 \ 0.8 & 0.4 end{pmatrix}$ . $M_{11}$ 為系統從狀態1(A超市)轉移到狀態1(A超市)的概率， $M_{12}$ 表示系統從狀態2(B超市)轉移到狀態1(A超市)的概率，等等。因為我們的狀態空間是完備的，也就是系統只能從一個狀態轉移到另外一個狀態，所以馬爾科夫矩陣的每一列求和都為1, 也就是 $sum_{i} M_{ij} = 1$ . 用矩陣表示為

$oldsymbol{1}^{T} M = oldsymbol{1}^{T}$

這個性質決定了馬爾科夫矩陣一定有一個值為1的特徵值。而且可以根據 Gershgorin circle theorem 得到馬爾科夫矩陣的任意一個特徵值的絕對值都不可能大於一。為了得到馬爾科夫矩陣的穩態解，我們需要計算馬爾科夫矩陣的冪。當冪次為無窮大時，如果得到的矩陣是收斂的話，我們就把得到的矩陣稱作是馬爾科夫鏈的穩態矩陣。如果馬爾科夫矩陣只有一個絕對值等於1的特徵值，那麼得到的冪矩陣一定是收斂的。相反，如果馬爾科夫矩陣還有其它絕對值等於1的特徵值，那麼這個馬爾科夫鏈就沒有穩態解。通過冪矩陣的方法計算馬爾科夫鏈穩態解的方法叫做power iteration method. Power iteration - Wikipedia 這個方法其實就是在計算馬爾科夫矩陣特徵值為1所對應的特徵向量，該特徵向量又叫做 Perron-Frobenius 向量。記時刻馬爾科夫鏈的狀態向量為 $oldsymbol{xi}_{n}$ ，我們有狀態向量的遞歸關係：

$oldsymbol{xi}_{n+1} = M oldsymbol{xi}_{n}$

假設 $oldsymbol{xi}_{infty} = lim_{n ightarrow infty} oldsymbol{xi}_{n}$ 存在，對兩邊同時求極限，得到

$oldsymbol{xi}_{infty} = M oldsymbol{xi}_{infty}$

所以穩態向量 $oldsymbol{xi}_{infty}$ 為馬爾科夫矩陣對應著特徵值為1的特徵向量，該向量又可以看做是映射 $M: mathbb{R}^{n} ightarrow mathbb{R}^{n}$ 的不動點。因此，只要我們知道了特徵值為1的特徵向量，我們就可以得到馬爾科夫鏈的穩態解(如果這個穩態解存在的話)。下面我將要簡單介紹一下power iteration method.

Power iteration method 為，對於一個可對角化的矩陣，為了計算它絕對值為最大的那個特徵值所對應的特徵向量，我們可以使用這個遞歸公式：

$oldsymbol{v}_{k+1} = frac{A oldsymbol{v}_{k}}{|| A oldsymbol{v}_{k} ||}, k ge 0$

實際用這個演算法計算特徵向量的時候，我們需要先選取一個任意非零的初始向量 $oldsymbol{v}_{0}$ ，然後代入上面的式子做迭代，直到連續兩次的迭代結果差異在誤差範圍內。用Python實現如下:

import numpy as np

def power_iteration(A):
(row, col) = A.shape
assert(row == col)
dimension = row
v = np.random.rand(dimension)
counter = 0
iteration_max = 1000
eps = 1.0e-15
while(counter < iteration_max):
counter += 1
v_updated = A.dot(v)
v_updated = v_updated/np.linalg.norm(v_updated)
error = min(np.linalg.norm(v - v_updated), np.linalg.norm(v + v_updated))
print "counter = " + str(counter) + ", error = " + str(error)
if (error < eps):
break
v = v_updated
eigenvalue = v_updated.dot(A).dot(v_updated)/v_updated.dot(v_updated)
return (eigenvalue, v_updated)

經過測試，該演算法可以給出正確的最大特徵值和它所對應的特徵向量。測試代碼如下：

import numpy as np def generate_matrix(dimension): assert(dimension >= 2) A = np.zeros((dimension, dimension)) for i in range(dimension): for j in range(dimension): A[i, j] = np.random.uniform() return A

def find_leading_eig(A):
(eigenvalues, eigenvectors) = np.linalg.eig(A)
leading_eigenvalue = abs(eigenvalues[0])
leading_index = 0
for i in range(1, len(eigenvalues)):
if (leading_eigenvalue < abs(eigenvalues[i])):
leading_index = i
leading_eigenvalue = abs(eigenvalues[i])
return (eigenvalues[leading_index], eigenvectors[:, leading_index])

def is_proportional(a, b):
assert(len(a) == len(b))
eps = 1.0e-15
result = []
for i in range(len(a)):
if (abs(a[i]) < eps and abs(b[i]) < eps):
result.append(None)
continue
elif (abs(a[i]) < eps or abs(b[i]) < eps):
return False
else:
result.append(a[i]/b[i])
#print result
assert(len(result) == len(a))
result = filter(lambda ele: ele is not None, result)
if (len(result) == 0):
return True
for i in range(1, len(result)):
if (abs(result[i] - result[0]) > 1.0e-12):
return False
return True

def main():
import sys
if (len(sys.argv) != 2):
print "dimension = sys.argv[1]. "
return -1
dimension = int(sys.argv[1])
A = generate_matrix(dimension)
(eigenvalue, eigenvector) = power_iteration(A)
print "leading eigenvalue = " + str(eigenvalue)
print "leading eigenvector = "
print eigenvector
(leading_eigenvalue, leading_eigenvector) = find_leading_eig(A)
print leading_eigenvalue
print leading_eigenvector
print "eigenvalue error = " + str(abs(eigenvalue - leading_eigenvalue))
print "consistent leading eigenvector? " + str(is_proportional(eigenvector, leading_eigenvector))
return 0

if __name__ == "__main__":
import sys
sys.exit(main())

使用這個演算法計算矩陣的最大特徵值有它的侷限性。為了知道它的侷限性，我們首先需要知道這個演算法為什麼成立。對於一個可對角化的矩陣，它有個線性無關的特徵向量，這個特徵向量可以張滿整個維歐氏空間。也就是任意一個維向量都可以唯一地表示為這個特徵向量的線性組合：

$oldsymbol{v} = sum_{i = 1}^{n}alpha_{i} oldsymbol{u}_{i}$

其中 $alpha_{i}$ 為表示係數， $oldsymbol{u}_{i}$ 為矩陣的特徵向量。這些特徵向量線性獨立，但是未必彼此正交。將作用到這個向量上，得到

$A^noldsymbol{v} = sum_{i = 1}^{n} alpha_{i}lambda_{i}^{n}oldsymbol{u}_{i}$

其中 $Aoldsymbol{u}_{i} = lambda_{i}oldsymbol{u}_{i}$ . 記絕對值最大的矩陣特徵值為 $lambda_{max}$ . 將上面的式子兩邊同時除以 $lambda_{max}^{n}$ ，得到

$lambda_{max}^{-n} A^noldsymbol{v} = sum_{i = 1}^{n}alpha_{i}Big( frac{lambda_{i}}{lambda_{max}} Big)^{n}oldsymbol{u}_{i}$

對於任意的特徵值，我們都有 $Biggvertfrac{lambda_{i}}{lambda_{max}}Biggvert le 1$ . 進一步假設除了 $lambda_{max}$ 之外，矩陣所有其他的特徵值的絕對值都小於 $|lambda_{max}|$ . 那麼當時，我們有

$lim_{n ightarrow infty} lambda_{max}^{-n} A^noldsymbol{v} = alpha oldsymbol{u}_{max}$

這裡， $oldsymbol{u}_{max}$ 為特徵值 $lambda_{max}$ 所對應的特徵向量。因為我們只關心特徵向量的方向，所以特徵向量乘以任意一個非零係數後得到的結果仍然是特徵向量。實際使用該演算法時，我們很有可能不知道 $lambda_{max}$ 的值，所以我們不能通過計算 $lambda_{max}^{-n} A^noldsymbol{v}$ 來得到矩陣的特徵向量，而是每計算一次都把它歸一化一次，這樣我們就可以在不知道最大特徵值的前提下就可以得到矩陣的最大特徵向量。知道了最大特徵向量後，我們就可以根據定義得到最大特徵值為:

$lambda_{max} = frac{oldsymbol{u}_{max}^T A oldsymbol{u}_{max} }{oldsymbol{u}_{max}^T oldsymbol{u}_{max}}$

在推導這個演算法的時候，我們做了兩個假設.

假設一：矩陣是可對角化的。
假設二：矩陣的絕對值為最大的特徵值是唯一的。

這兩個假設未必成立。有的矩陣是不可對角化的。不可對角化的矩陣稱作是defective matrix. 對於不可對角化的矩陣，它沒有個線性無關的特徵向量。這種矩陣的特徵多項式一定存在重根。記特徵多項式為 . 如果 $lambda_{1}$ 是該特徵多項式的一個重根，那麼特徵多項式必定滿足這個條件：

$f(lambda) propto (lambda - lambda_1)^{m_1}$

其中是特徵值的代數重數，方程組的非零解的個數稱作是該特徵值的幾何重數。幾何重數不可能大於代數重數。如果矩陣所有特徵值的幾何重數都等於它的代數重數，那麼這個矩陣就是可對角化的，否則這個矩陣是不可對角化的，這時它只能通過相似變換變到若爾當標準型，而不能變到對角型矩陣。如果一個矩陣是不可對角化的，那麼它的特徵向量就不能張滿整個歐氏空間，於是我們的基本假設 $oldsymbol{v} = sum_{i = 1}^{n}alpha_{i} oldsymbol{u}_{i}$ 就失效了。為了使用這個演算法，首先我們假設馬爾科夫矩陣是可對角化的。數值上這個假設並沒有給我們的問題加上一個特彆強的限制，因為我們知道可對角化的矩陣是稠密的，因此我們可以假設在機器精度範圍內，所有的矩陣都是可對角化的。在使用power iteration method的時候，我們需要進一步假設絕對值為最大的特徵值是唯一的，這個條件給我們的演算法加了一個比較強的限制。之所以需要這個條件，是因為我們使用這個演算法的時候，需要假設 $Biggvert frac{lambda_{i}}{lambda_{max}} Biggvert < 1$ 對任意的 $lambda_{i} e lambda_{max}$ 都成立。如果這個條件不成立，那麼這個演算法只能給出所有特徵值的絕對值等於 $vert lambda_{max} vert$ 的特徵向量的線性組合，而不能得到這些獨立的特徵向量。另外，該演算法的收斂速度取決於 $Biggvert frac{lambda_2}{lambda_{max}}Biggvert$ 的值，其中為絕對值為第二大的特徵值。因為第二大的特徵值決定了power iteration method的收斂速度，所以有很多專門的文章研究這個第二大的特徵值。

對於馬爾科夫矩陣，我們已經知道它最大的特徵值為1，其餘的特徵值的絕對值要麼小於1，要麼等於1. 如果我們假設馬爾科夫矩陣絕對值等於1的特徵值的只有一個，那麼我們就可以用power iteration method去計算它的特徵向量。因為我們知道 $lambda_{max} = 1$ ，所以我們可以直接計算 $oldsymbol{xi}_{n+1} = M oldsymbol{xi}_{n}$ 得到最大特徵向量，而不需要每次都對得到的向量做歸一化操作。這就等價於我們只需要計算 $M^{n}oldsymbol{xi}_{0}$ 就可以得到特徵值所對應的特徵向量，其中 $oldsymbol{xi}_{0}$ 是任意的初始非零向量。當馬爾科夫矩陣的只有一個絕對值為1的特徵值時，該馬爾科夫鏈的穩態解唯一，並且不依賴於初始向量，我們稱這樣的馬爾科夫鏈為各態遍歷的(ergodic)。如果馬爾科夫矩陣有超過一個絕對值等於1的特徵值，那麼這個馬爾科夫鏈就沒有穩態解。一個簡單的例子是 $M = egin{pmatrix} 0 & 1 \ 1 & 0 end{pmatrix}$ . 如果仍然用兩個超市的顧客人數為例，這個矩陣表示每次去A超市的顧客下次會全部去B超市，而每次去B超市的顧客每次都全部去A超市。這是一個振蕩的過程，這個馬爾科夫鏈沒有穩態解，也就是去A超市和B超市的人數比例不會趨近於一個常數。對應到矩陣 $M = egin{pmatrix} 0 & 1 \ 1 & 0 end{pmatrix}$ ，它有兩個特徵值，一個是 1, 另一個是 . 因為它有兩個絕對值為1的特徵值，所以power iteration method不能給出穩態解，這時也沒有穩態解。

一種特殊的馬爾科夫矩陣是對稱馬爾科夫矩陣。對稱馬爾科夫矩陣的所有列的和與所有行的和都是1，因此 $oldsymbol{1} = (1, 1, ..., 1)^{T}$ 是該矩陣對應特徵值為1的特徵向量，也是該矩陣的穩態解(如果這個矩陣存在穩態解的話)。因此，如果一個馬爾科夫鏈的轉移矩陣是對稱矩陣，並且存在穩態解，那麼經過足夠長的時間後，每個狀態都會被等概率地訪問到。

第三節：PageRank演算法

PageRank演算法最早被Google用來做網頁排序。這個演算法在很多領域都有著廣泛的應用。我們可以將網頁看作是一個有向圖，圖的節點是網頁，如果網頁A有一個鏈接指向網頁B，那麼在節點A和節點B之間我們就可以建立一條有向邊。Google最開始試圖發明一種演算法，該演算法可以給所有的網頁按照重要度排序。直覺上，我們可以認為如果一個網頁被很多網頁鏈接，那麼這個網頁應該比較重要。如果一個網頁被另一個重要的網頁鏈接了，那麼這個被鏈接的網頁也應該比較重要。PageRank演算法通過分析網頁之間的關聯來給網頁排序。這個演算法僅僅分析網頁之間的關聯，而不考慮網頁的內容，因此也有一定的侷限性。接下來，我要簡單的介紹一下這個演算法。我們很容易發現，這個演算法本質上就是馬爾科夫鏈的應用。

假設我們有如下的幾個網頁：