譜分群

拉普拉斯矩陣在實際應用上，常會被用來模擬一個變動的系統，如互相連接的儲水槽、彈黃系統、或是社群網路。因此每一條邊可能會有權重，來描述相對應的輸水量、彈力係數、熟識度等等。這些權重一般來說都是正的，而愈大的值表示兩個點的關係愈緊密。

一個 賦權圖（weighted graph） 指的是一個簡單圖 \(G\) 再搭配一個權重函數 \(w: E\rightarrow\mathbb{R}^+\)。每條邊 \(e = \{i,j\}\) 上的權重常簡寫為 \(w(e) = w_{i,j}\)。如此一來，我們可以定義一個賦權圖 \((G,w)\) 的拉普拉斯矩陣 \(L(G,w) = \begin{bmatrix} \ell_{i,j} \end{bmatrix}\) 為

\[ \ell_{i,j} = \begin{cases} -w_{i,j} & \text{ if } \{i,j\}\in E(G),\ i\neq j \\ 0 & \text{ if } \{i,j\}\notin E(G),\ i\neq j \\ \sum_{k: k\sim i}w_{i,k} & \text{ if } i = j. \end{cases} \]

這樣的定義下，無權重的簡單圖相當於視每條邊的權重均為 \(1\)；而另一方面在賦權圖裡，我們也把點的度數推廣成周圍各邊的權重相加。

舉例來說，如果 \(P_3\) 上的兩條邊權重分別為 \(0.1\) 和 \(2\)，則其拉普拉斯矩陣為

\[ L(P_3) = \begin{bmatrix} 0.1 & -0.1 & 0 \\ -0.1 & 2.1 & -2 \\ 0 & -2 & 2 \end{bmatrix}. \]

我們同樣可以觀察到

\[ \begin{bmatrix} 0.1 & -0.1 & 0 \\ -0.1 & 2.1 & -2 \\ 0 & -2 & 2 \end{bmatrix} = 0.1\begin{bmatrix} 1 & -1 & 0 \\ -1 & 1 & 0 \\ 0 & 0 & 0 \end{bmatrix} + 2\begin{bmatrix} 0 & 0 & 0 \\ 0 & 1 & -1 \\ 0 & -1 & 1 \end{bmatrix}. \]

因此其二次型可以寫為

\[ \bx\trans L(P_3)\bx = 0.1(x_1 - x_2)^2 + 2(x_2 - x_3)^2. \]

用類似的手法我們可以知道任一賦權圖 \((G,w)\) 的拉普拉斯矩陣二次型可寫為

\[ \bx\trans L(G)\bx = \sum_{\{i,j\}\in E(G)} w_{i,j}(x_i - x_j)^2. \]

而在所有邊權重都大於零的前提下，所有賦權拉普拉斯矩陣都保有相同的基本性質。

性質

令 \(L = L(G)\) 為賦權圖 \(G\) 的拉普拉斯矩陣，則

  1. \(L\) 是半正定矩陣，因此 \(\ker(L) = \{\bx : \bx\trans L\bx = 0\}\)。
  2. \(\bx\trans L\bx = 0\) 發生時，\(\bx\) 在 \(G\) 上同一個連通區塊的值都一樣。
  3. \(\nul(L)\) 等於 \(G\) 的連通區塊的個數。

值得注意的是，賦權圖的拉普拉斯矩陣可以用來描述所有滿足以下條件的矩陣 \(L\)：

  1. \(L = L\trans\)。
  2. \(L\) 的非對角線項非零即負。
  3. \(L\bone = \bzero\)。

是一個很大一類的矩陣。

一但加入權重之後，邊就不在是全有全無的選擇，我們可能在某處加入一段很細、權重很小的邊，而它只會對整個系統造成連續的細微改變。原先我們能用拉普拉斯矩陣來判斷圖的連通區塊，現在我們可以推廣成判斷圖的 群聚（cluster），差別在於群聚和群聚之間可能還是連通的，但相連的權重很微小。

我們考慮以下的例子。令 \(\Gamma_x\) 是一個 \(15\) 個點的賦權圖，把它的點分成三群 \(\{1,\ldots, 5\}\), \(\{6, \ldots, 10\}\), 以及 \(\{11, \ldots, 15\}\)，同一群內的任兩點相連，權重均為 \(1\)，之後我們加上兩條權重為 \(x\) 的邊 \(\{5,6\}\) 以及 \(\{10,11\}\)。當 \(x = 0\) 的時候我們就當作邊不存在。

根據拉普拉斯矩陣的基本性質，\(L(\Gamma_0)\) 的特徵值恰有三個 \(0\)，而 \(L(\Gamma_{0.1})\) 的特徵值只有一個 \(0\)。然而因為 \(L(\Gamma_{0.1})\) 是 \(L(\Gamma_0)\) 的輕微擾動，特徵值和特徵向量都不會改變太多，經過數值計算，以下列出兩個矩陣最小的三個特徵值以及對應的特徵向量，箭頭左側為 \(L(\Gamma_0)\) 的資訊、右側為 \(L(\Gamma_{0.1})\) 的資訊。（數值計算所用的程式碼附在本文末，有興趣的朋友可以動手試試看。）

\[ \{0,0,0\} \rightarrow \{0, 0.02, 0.06\} \]
\[ \begin{bmatrix} 0.45 & 0 & 0 \\ 0.45 & 0 & 0 \\ 0.45 & 0 & 0 \\ 0.45 & 0 & 0 \\ 0.45 & 0 & 0 \\ 0 & 0.45 & 0 \\ 0 & 0.45 & 0 \\ 0 & 0.45 & 0 \\ 0 & 0.45 & 0 \\ 0 & 0.45 & 0 \\ 0 & 0 & 0.45 \\ 0 & 0 & 0.45 \\ 0 & 0 & 0.45 \\ 0 & 0 & 0.45 \\ 0 & 0 & 0.45 \\ \end{bmatrix} \rightarrow \begin{bmatrix} 0.26 & 0.32 & -0.18 \\ 0.26 & 0.32 & -0.18 \\ 0.26 & 0.32 & -0.18 \\ 0.26 & 0.32 & -0.18 \\ 0.26 & 0.31 & -0.17 \\ 0.26 & 0.01 & 0.36 \\ 0.26 & -0.00 & 0.37 \\ 0.26 & -0.00 & 0.37 \\ 0.26 & -0.00 & 0.37 \\ 0.26 & -0.01 & 0.36 \\ 0.26 & -0.31 & -0.17 \\ 0.26 & -0.32 & -0.18 \\ 0.26 & -0.32 & -0.18 \\ 0.26 & -0.32 & -0.18 \\ 0.26 & -0.32 & -0.18 \\ \end{bmatrix} \]

我們發現 \(L(\Gamma_{0.1})\) 確實有三個很接近 \(0\) 的特徵值，所以我們可以用「接近 \(0\) 的特徵值的個數」來判斷「群聚的個數」。儘管這裡的「接近」、「群聚」都不是明確定義的觀念，但在實務應用上它還是給了不錯的資訊。

另一方面，我們發現特徵向量並沒有逐項地連續變動，這是因為對 \(L(\Gamma_{0})\) 來說，\(0\) 這個特徵值重數有 \(3\)，所以它的基底選取不唯一，如果基底選取適當的話還是可以把它看成連續變動。但正因為我們不知道會選到哪一個基底的連續變動，我們沒辦法像以前一樣用非零的位置來判斷群聚的範圍。然而，如果我們觀察 \(L(\Gamma_{0.1})\) 的特徵向量，把三個特徵向量擺一起，然後橫著看成 \(15\) 個三維空間的點，會發現前五個點的座標很接近、中間五個點很接近、最後五個點也很接近。自然而然也把圖的點分成三個群聚。

考慮任意一個賦權圖 \((G,w)\) 的拉普拉斯矩陣 \(L\)，以上的觀察整理如下：

  1. \(L\) 的特徵值中接近 \(0\) 的個數，可以用來估測 \((G,w)\) 的群聚數。
  2. 將 \(L\) 的小特徵值所對應到的特徵向量擺一起、再橫著看成樣本點，對這些樣本點運行任何的分群演算法（如 \(k\)-means clustering algorithm），則可以找到 \((G,w)\) 的群聚範圍。

這種把圖的每一點嵌入到歐氏空間中的方法，被稱為 圖嵌入（spectral embedding 或是 Laplacian eigenmap），而將圖嵌入及分群演算法合在一起的動作叫 譜分群（spectral clustering）。

值得一提的是，\(L\) 永遠有一個特徵值 \(0\)，其對應到的特徵向量是 \(\frac{1}{\sqrt{n}}\bone\)。由於這個特徵向量每一項都是一樣的，在圖嵌入中沒有提供任何有用的資訊，因此一般來說會從第二小的特徵值的特徵向量開始選取。

想想以下問題：

  1. 說明為什麼 \(\bx\trans L(G)\bx = \sum_{\{i,j\}\in E(G)} w_{i,j}(x_i - x_j)^2\)。
  2. 找一個可以執行 Python 程式碼的環境，試著跑看看文末的程式。把連結群聚的邊權重改得更小看看，觀察它的特徵值如何變化。
  3. 查找相關資訊，並描述 \(k\)-means clustering 演算法是如何運行的。
import numpy as np
import matplotlib.pyplot as plt
np.set_printoptions(precision=2, suppress=True)

# G has three clusters 0~4, 5~9, 10~14
# edges inside the same cluster are weighted 1
# the edge {4,5} and {9,10} are weighted 0.1
L = np.zeros((15,15), dtype=float)
L[:5, :5] = -1
L[5:10, 5:10] = -1
L[10:, 10:] = -1
L[[4,5,9,10], [5,4,10,9]] = -0.1
L[np.arange(15), np.arange(15)] -= L.sum(axis=1)
vals,vecs = np.linalg.eigh(L)
plt.plot(np.arange(1,16), vals)
print(vals)
print(vecs[:,:3])